Lorenz Attractor With Runge-kutta Python

Hello I have to program a python function to solve Lorenz differential equations using Runge-Kutta 2cond grade sigma=10, r=28 and b=8/3 with initial conditions (x,y,z)=(0,1,0) t

Solution 1:

Solving differential equations numerically can be challenging. If you choose too high step sizes, the solution will accumulate high errors and can even become unstable, as in your case.

Either you should drastically reduce the step size (h) or just use the adaptive Runge Kutta method provided by scipy:

from numpy import array, linspace
from scipy.integrate import solve_ivp
import pylab
from mpl_toolkits import mplot3d

deffunc(t, r):
    x, y, z = r 
    fx = 10 * (y - x)
    fy = 28 * x - y - x * z
    fz = x * y - (8.0 / 3.0) * z
    return array([fx, fy, fz], float)

# and I run it
r0 = [0, 1, 0]
sol = solve_ivp(func, [0, 50], r0, t_eval=linspace(0, 50, 5000))

# and plot it
fig = pylab.figure()
ax = pylab.axes(projection="3d")
ax.plot3D(sol.y[0,:], sol.y[1,:], sol.y[2,:], 'blue')
pylab.show()

This solver uses 4th and 5th order Runge Kutta combination and controls the deviation between them by adapting the step size. See more usage documentation here: https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html

Solution 2:

You use a step size of h=2.5.

For RK4 the useful step sizes given a Lipschitz constant L are in the range L*h=1e-3 to 0.1, one might get somewhat right looking results up to L*h=2.5. Above that the method turns chaotic, any resemblance to the underlying ODE is lost.

The Lorenz system has a Lipschitz constant of about L=50, see Chaos and continuous dependency of ODE solution, so h<0.05 is absolutely required, h=0.002 is better and h=2e-5 gives the numerically best results for this numerical method.

Solution 3:

It can be related to a division by zero or when a limit of a type is exceeded (float type).

To figure out where and when it happens you can set numpy.seterr('raise') and it will raise an exception so you can debug and see what it's happening. It seems your algorithm is diverging.

Here you can se how to use numpy.seterr