I am currently trying to write some python code to solve an arbitrary system of first order ODEs, using a general explicit Runge-Kutta method defined by the values alpha, gamma (both vectors of dimension m) and beta (lower triangular matrix of dimension m x m) of the Butcher table which are passed in by the user. My code appears to work for single ODEs, having tested it on a few different examples, but I’m struggling to generalise my code to vector valued ODEs (i.e. systems).

In particular, I try to solve a Van der Pol oscillator ODE (reduced to a first order system) using Heun’s method defined by the Butcher Tableau values given in my code, but I receive the errors

- “RuntimeWarning: overflow encountered in double_scalars
`f = lambda t,u: np.array(... etc)`

” and- “RuntimeWarning: invalid value encountered in add
`kvec[i] = f(t+alpha[i]*h,y+h*sum)`

”

followed by my solution vector that is clearly blowing up. Note that the commented out code below is one of the examples of single ODEs that I tried and is solved correctly. Could anyone please help? Here is my code:

import numpy as np def rk(t,y,h,f,alpha,beta,gamma): '''Runga Kutta iteration''' return y + h*phi(t,y,h,f,alpha,beta,gamma) def phi(t,y,h,f,alpha,beta,gamma): '''Phi function for the Runga Kutta iteration''' m = len(alpha) count = np.zeros(len(f(t,y))) kvec = k(t,y,h,f,alpha,beta,gamma) for i in range(1,m+1): count = count + gamma[i-1]*kvec[i-1] return count def k(t,y,h,f,alpha,beta,gamma): '''returning a vector containing each step k_{i} in the m step Runga Kutta method''' m = len(alpha) kvec = np.zeros((m,len(f(t,y)))) kvec[0] = f(t,y) for i in range(1,m): sum = np.zeros(len(f(t,y))) for l in range(1,i+1): sum = sum + beta[i][l-1]*kvec[l-1] kvec[i] = f(t+alpha[i]*h,y+h*sum) return kvec def timeLoop(y0,N,f,alpha,beta,gamma,h,rk): '''function that loops through time using the RK method''' t = np.zeros([N+1]) y = np.zeros([N+1,len(y0)]) y[0] = y0 t[0] = 0 for i in range(1,N+1): y[i] = rk(t[i-1],y[i-1], h, f,alpha,beta,gamma) t[i] = t[i-1]+h return t,y ################################################################# '''f = lambda t,y: (c-y)**2 Y = lambda t: np.array([(1+t*c*(c-1))/(1+t*(c-1))]) h0 = 1 c = 1.5 T = 10 alpha = np.array([0,1]) gamma = np.array([0.5,0.5]) beta = np.array([[0,0],[1,0]]) eff_rk = compute(h0,Y(0),T,f,alpha,beta,gamma,rk, Y,11)''' #constants mu = 100 T = 1000 h = 0.01 N = int(T/h) #initial conditions y0 = 0.02 d0 = 0 init = np.array([y0,d0]) #Butcher Tableau for Heun's method alpha = np.array([0,1]) gamma = np.array([0.5,0.5]) beta = np.array([[0,0],[1,0]]) #rhs of the ode system f = lambda t,u: np.array([u[1],mu*(1-u[0]**2)*u[1]-u[0]]) #solving the system time, sol = timeLoop(init,N,f,alpha,beta,gamma,h,rk) print(sol)

Your step size is not small enough. The Van der Pol oscillator with `mu=100`

is a fast-slow system with very sharp turns at the switching of the modes, so rather stiff. With explicit methods this requires small step sizes, the smallest sensible step size is `1e-5`

to `1e-6`

. You get a solution on the limit cycle already for `h=0.001`

, with resulting velocities up to 150.

You can reduce some of that stiffness by using a different velocity/impulse variable. In the equation

x'' - mu*(1-x^2)*x' + x = 0

you can combine the first two terms into a derivative,

mu*v = x' - mu*(1-x^2/3)*x

so that

x' = mu*(v+(1-x^2/3)*x) v' = -x/mu

The second equation is now uniformly slow close to the limit cycle, while the first has long relatively straight jumps when `v`

leaves the cubic `v=x^3/3-x`

.

This integrates nicely with the original `h=0.01`

, keeping the solution inside the box `[-3,3]x[-2,2]`

, even if it shows some strange oscillations that are not present for smaller step sizes and the exact solution.

## Recent Comments