Using python built-in functions for coupled ODEs

Question

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 are given below.

# n-dimensional change in omega
def d_theta(omega):
    return omega

# n-dimensional change in omega
def d_omega(K,A,P,alpha,mask,n):
    def layer1(theta,omega):
        T = theta[:,None] - theta
        A[mask] = K[mask] * np.sin(T[mask])
        return - alpha*omega + P - A.sum(1)
    return layer1

These equations return vectors.

QUESTION 1

I know how to use odeint for two dimensions, (y,t). for my research I want to use a built-in Python function that works for higher dimensions.

QUESTION 2

I do not necessarily want to stop after a predetermined amount of time. I have other stopping conditions in the code below that will indicate whether the system of equations converges to the steady state. How do I incorporate these into a built-in Python solver?

WHAT I CURRENTLY HAVE

This is the code I am currently using to solve the system. I just implemented RK4 with constant time stepping in a loop.

# This function randomly samples initial values in the domain and returns whether the solution converged

# Inputs:
# f             change in theta (d_theta)
# g             change in omega (d_omega)
# tol           when step size is lower than tolerance, the solution is said to converge
# h             size of the time step
# max_iter      maximum number of steps Runge-Kutta will perform before giving up
# max_laps      maximum number of laps the solution can do before giving up
# fixed_t       vector of fixed points of theta
# fixed_o       vector of fixed points of omega
# n             number of dimensions

# theta         initial theta vector
# omega         initial omega vector

# Outputs:
# converges     true if it nodes restabilizes, false otherwise

def kuramoto_rk4_wss(f,g,tol_ss,tol_step,h,max_iter,max_laps,fixed_o,fixed_t,n):
    def layer1(theta,omega):
        lap = np.zeros(n, dtype = int)
        converges = False
        i = 0
        tau = 2 * np.pi

        while(i < max_iter): # perform RK4 with constant time step
            p_omega = omega
            p_theta = theta
            T1 = h*f(omega)
            O1 = h*g(theta,omega)
            T2 = h*f(omega + O1/2)
            O2 = h*g(theta + T1/2,omega + O1/2)
            T3 = h*f(omega + O2/2)
            O3 = h*g(theta + T2/2,omega + O2/2)
            T4 = h*f(omega + O3)
            O4 = h*g(theta + T3,omega + O3)

            theta = theta + (T1 + 2*T2 + 2*T3 + T4)/6 # take theta time step

            mask2 = np.array(np.where(np.logical_or(theta > tau, theta < 0))) # find which nodes left [0, 2pi]
            lap[mask2] = lap[mask2] + 1 # increment the mask
            theta[mask2] = np.mod(theta[mask2], tau) # take the modulus

            omega = omega + (O1 + 2*O2 + 2*O3 + O4)/6

            if(max_laps in lap): # if any generator rotates this many times it probably won't converge
                break
            elif(np.any(omega > 12)): # if any of the generators is rotating this fast, it probably won't converge
                break
            elif(np.linalg.norm(omega) < tol_ss and # assert the nodes are sufficiently close to the equilibrium
               np.linalg.norm(omega - p_omega) < tol_step and # assert change in omega is small
               np.linalg.norm(theta - p_theta) < tol_step): # assert change in theta is small
                converges = True
                break
            i = i + 1
        return converges
    return layer1

Thanks for your help!

For the lap counting, you should use the `np.floor_divide` function complimentary to the `np.mod`. Except if you are sure that the rotation is always forward. — Lutz Lehmann, Feb 24 '20 at 11:25
`np.floor_divide` will output negative numbers for negative inputs — Damien Beecroft, Feb 24 '20 at 16:32
This is bad why? Currently if you have an oscillation around the cutoff, say from 0.1 to -0.1, you will recognize a lap and change to 2*pi-0.1. Now if in the next instance the oscillation goes up by 0.2, you recognize another lap and change back to 0.1. Thus a small oscillation can quickly end the computation without showing what you are looking for. Alternatively, you have to check if the absolute value of the lap count is over the threshold. — Lutz Lehmann, Feb 24 '20 at 16:46

Lutz Lehmann · Accepted Answer · 2020-02-24T17:01:32.350

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

def odesys(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.

Can you explain your code a bit? What is the point of the second line and why did you write the `*` in front of the return function? Also, is N the length of the arrays? — Damien Beecroft, Feb 24 '20 at 16:42
Yes, `N=n`, I was writing without checking. `*a` where a is a list-like object unpacks the elements of the list, so `[*a,*b]` is the combined list of the elements of `a` and then `b`. The second line is just for reading convenience and code maintainability, there is no other reason to not use the slices directly in the function calls. — Lutz Lehmann, Feb 24 '20 at 16:50

Using python built-in functions for coupled ODEs

1 Answers1

Linked