I am trying to do Lorenz's equations in python (I am following exercise 8.3 from Mark Newman - Computational Physics (2012, CreateSpace Independent Publishing Platform)) I already got the graphics and everything looks "correct". This is probably a math problem and not really a progamming problem but I'm posting here to make sure. First, this is my code:
from numpy import arange,array
import pylab
def f(v,t):
s=10
r=28
b=8/3
x= v[0]
y= v[1]
z= v[2]
fx= s*(y - x)
fy= r*x - y - x*z
fz= x*y - b*z
return array([fx,fy,fz],float)
def d(N):
a=0.0
b=50.0
h=(b-a)/N
r=array([0.0,1.0,0.0],float)
tpoints=arange(a,b,h)
xpoints= []
ypoints= []
zpoints= []
for t in tpoints:
xpoints.append(r[0])
ypoints.append(r[1])
zpoints.append(r[2])
k1 = h*f(r,t)
k2 = h*f(r+0.5*k1,t+0.5*h)
k3 = h*f(r+0.5*k2,t+0.5*h)
k4 = h*f(r+k3,t+h)
r += (k1+2*k2+2*k3+k4)*(1/6)
return tpoints,xpoints,ypoints,zpoints
for i in range (1,6):
N=10**i
pylab.plot(d(N)[0],d(N)[1], label=N)
pylab.xlabel("t")
pylab.ylabel("x(t)")
pylab.title("Gráficos x em função de t")
pylab.legend()
pylab.show()
pylab.plot(d(N)[0],d(N)[2], label=N)
pylab.xlabel("t")
pylab.ylabel("y(t)")
pylab.title("Gráficos y em função de t")
pylab.legend()
pylab.show()
pylab.plot(d(N)[0],d(N)[3], label=N)
pylab.xlabel("t")
pylab.ylabel("z(t)")
pylab.title("Gráficos z em função de t")
pylab.legend()
pylab.show()
pylab.plot(d(N)[1],d(N)[3], label=N)
pylab.xlabel("x")
pylab.ylabel("z(x)")
pylab.title("Gráficos z em função de x")
pylab.legend()
pylab.show()
This give me the graphs to solve the problem and I think this is correct.
When I go from i=1
to i=3
it give me this error:
I think this is related with math problem, but when I search for the error it got me with something with arrays. So I'm checking it.
For i
equal or more then 4, it gave me no problem.
The Lorenz system with RK4 requires a step size of about 0.05
or less, that is, N=10**4
or larger for your construction. See the near duplicate Lorenz attractor with Runge-Kutta python.
For larger step sizes, that is, the cases that you get the error in, the method will return chaotic results that are mostly unrelated to the exact solution of the system and any bounds related to it. So divergence with floating point overflow is possible due to the quadratic, super-linear terms.