Using Python Built-in Functions For Coupled Odes

THIS PART IS JUST BACKGROUND IF YOU NEED IT I am developing a numerical solver for the Second-Order Kuramoto Model. The functions I use to find the derivatives of theta and omega a

Solution 1:

You can wrap your existing functions into a function accepted by odeint (option tfirst=True) and solve_ivp as

defodesys(t,u):
    theta,omega = u[:n],u[n:]; # or = u.reshape(2,-1);return [ *f(omega), *g(theta,omega) ]; # or np.concatenate([f(omega), g(theta,omega)])

u0 = [*theta0, *omega0]
t = linspan(t0, tf, timesteps+1);
u = odeint(odesys, u0, t, tfirst=True);

#or

res = solve_ivp(odesys, [t0,tf], u0, t_eval=t)

The scipy methods pass numpy arrays and convert the return value into same, so that you do not have to care in the ODE function. The variant in comments is using explicit numpy functions.

While solve_ivp does have event handling, using it for a systematic collection of events is rather cumbersome. It would be easier to advance some fixed step, do the normalization and termination detection, and then repeat this.

If you want to later increase efficiency somewhat, use directly the stepper classes behind solve_ivp.