Solving The Lorentz Model Using Runge Kutta 4th Order In Python Without A Package

I wish to solve the Lorentz model in Python without the help of a package and my codes seems not to work to my expectation. I do not know why I am not getting the expected results

Solution 1:

I changed the integration step (btw., classical 4th order Runge-Kutta, not the adaptive RK54) to use the python core concept of lists and list operations extensively to reduce the number of places where the computation is defined. There were no errors there to correct, but I think the algorithm itself is more concentrated.

You had an error in the system that changed it into a system that rapidly diverges. You had fx = sigma*(y-z) while the Lorentz system has fx = sigma*(y-x).

Next your main loop has some strange assignments. In every loop you first set the previous coordinates to the initial conditions and then replace the full arrays with the RK step applied to the full arrays. I replaced that completely, there are no small steps to a correct solution.

import numpy as np
import matplotlib.pyplot as plt
#from scipy.integrate import odeintdeffx(x,y,z,t): return sigma*(y-x)
deffy(x,y,z,t): return x*(rho-z)-y
deffz(x,y,z,t): return x*y-beta*z

#a) Defining the classical Runge-Kutta 4th order methoddefRungeKutta45(x,y,z,fx,fy,fz,t,h):
    k1x,k1y,k1z = ( h*f(x,y,z,t) for f in (fx,fy,fz) )
    xs, ys,zs,ts = ( r+0.5*kr for r,kr inzip((x,y,z,t),(k1x,k1y,k1z,h)) )
    k2x,k2y,k2z = ( h*f(xs,ys,zs,ts) for f in (fx,fy,fz) )
    xs, ys,zs,ts = ( r+0.5*kr for r,kr inzip((x,y,z,t),(k2x,k2y,k2z,h)) )
    k3x,k3y,k3z = ( h*f(xs,ys,zs,ts) for f in (fx,fy,fz) )
    xs, ys,zs,ts = ( r+kr for r,kr inzip((x,y,z,t),(k3x,k3y,k3z,h)) )
    k4x,k4y,k4z  =( h*f(xs,ys,zs,ts) for f in (fx,fy,fz) )
    return (r+(k1r+2*k2r+2*k3r+k4r)/6for r,k1r,k2r,k3r,k4r inzip((x,y,z),(k1x,k1y,k1z),(k2x,k2y,k2z),(k3x,k3y,k3z),(k4x,k4y,k4z)))



sigma=10.
beta=8./3.
rho=28.
tIn=0.
tFin=10.
h=0.01
totalSteps=int(np.floor((tFin-tIn)/h))

t = totalSteps * [0.0]
x = totalSteps * [0.0]
y = totalSteps * [0.0]
z = totalSteps * [0.0]

x[0],y[0],z[0],t[0] = 1., 1., 1., 0.#Initial conditionfor i inrange(1, totalSteps):
    x[i],y[i],z[i] = RungeKutta45(x[i-1],y[i-1],z[i-1], fx,fy,fz, t[i-1], h)

Using tFin = 40 and h=0.01 I get the image

looking like the typical image of the Lorentz attractor.