Welcome to ShenZhenJia Knowledge Sharing Community for programmer and developer-Open, Learning and Share
menu search
person
Welcome To Ask or Share your Answers For Others

Categories

Please help me to understand why this code doesn't work. I know there is something very stupid wrong. This should be an implementation of the fourth order Runge kutta algorithm to solve Lorentz system of equations. I can't use any function because is an assignment.The errors are:

Warning (from warnings module): File "/Users/giuseppenegro/Desktop/loretnz.py", line 21 fy =x*ro-y-x*z RuntimeWarning: overflow encountered in double_scalars

Warning (from warnings module): File "/Users/giuseppenegro/Desktop/loretnz.py", line 22 fz =x*y-beta*z RuntimeWarning: overflow encountered in double_scalars

Warning (from warnings module): File "/Users/giuseppenegro/Desktop/loretnz.py", line 22 fz =x*y-beta*z RuntimeWarning: invalid value encountered in double_scalars

Warning (from warnings module): File "/Users/giuseppenegro/Desktop/loretnz.py", line 32 r+=(k1+2*k2+2*k3+k4)/6 RuntimeWarning: invalid value encountered in add

from numpy import array,arange,zeros
from pylab import plot,xlim,show
ro=28
beta=8/3
sigma=10
a=0
b=50
r=array([0.0,1.0,0.0],float)

N=10
h=(b-a)/N
tpoints=arange(a,b,h)
xpoints=[]
ypoints=[]
zpoints=[]
def f(r,t):
     x=r[0]
     y=r[1]
     z=r[2]
     fx =sigma*(y-x)
     fy =x*ro-y-x*z
     fz =x*y-beta*z
     return array([fx,fy,fz],float)
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,0.5*h+t)
    k3=h*f(r+0.5*k2+0.5,0.5*h+t)
    k4=h*f(r+k3,t+h)
    r+=(k1+2*k2+2*k3+k4)/6
See Question&Answers more detail:os

与恶龙缠斗过久,自身亦成为恶龙;凝视深渊过久,深渊将回以凝视…
thumb_up_alt 0 like thumb_down_alt 0 dislike
295 views
Welcome To Ask or Share your Answers For Others

1 Answer

The Runge-Kutta methods need to advance in small steps, 5 is not a small step. Try setting N to 1000.0 instead (the decimal is to make sure that (b-a)/N != 0).


与恶龙缠斗过久,自身亦成为恶龙;凝视深渊过久,深渊将回以凝视…
thumb_up_alt 0 like thumb_down_alt 0 dislike
Welcome to ShenZhenJia Knowledge Sharing Community for programmer and developer-Open, Learning and Share
...