2

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:

enter image description here

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.

martineau
  • 119,623
  • 25
  • 170
  • 301
Miguel Godinho
  • 125
  • 8
  • Do you understand what those warnings mean? Also, if you put a breakpoint right before you get each warning, what are all the values in those equations? Are they what you'd expect? – Random Davis Dec 09 '20 at 19:20
  • 1- I think my problem is that I 'm clear what these warnings means. I even thought it was errors. 2-what do u mean with putting "a breakpoint right before you get each warning" 3-"what are all the values in those equations?" didn't understand this. 4- If u mean I get what i want, with N=> 10000 yes, with N<= 1000, no. – Miguel Godinho Dec 09 '20 at 19:52
  • You're programming and you don't know about breakpoints? Learning to use a debugger is like critical day 1 info. – Random Davis Dec 09 '20 at 20:00
  • @RandomDavis : If using jupyter there are no breakpoints, one would need a proper IDE like pyzo or spyder etc. to have the usual debugging facilities available. In jupyter one could insert print statements inside the loops to see a list up to the last computed values. – Lutz Lehmann Dec 09 '20 at 20:06
  • @LutzLehmann `import pdb; pdb.set_trace()` works in any interpreter – Paul H Dec 09 '20 at 20:10
  • @PaulH : Thanks, I'll remember that. – Lutz Lehmann Dec 09 '20 at 20:14
  • I use jupyter, but this is math "problem" anyway, thx. – Miguel Godinho Dec 09 '20 at 20:21

1 Answers1

3

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.

Lutz Lehmann
  • 25,219
  • 2
  • 22
  • 51