Hi all,
I am trying to solve the famous Lorenz 3x3 system of ODEs using rk4.f from NR.

I studied an example from the NR example book and then I tried to code my program.
Here is the problem
C................................................. .....
C
C NUMERICAL SOLUTION OF THE FOLLOWING 3x3 SYSTEM
C OF ODEs
C LORENZ
C
C DX
C ---- = 10 * ( Y - X )
C DT
C
C DY
C ---- = X * (28 -Z) - Y
C DT
C
C DZ
C ---- = X * Y - (8/3) * Z
C DT
C
C WITH INITIAL CONDITIONS
C X = 3
C 0
C
C Y = 15
C 0
C
C Z = 1
C 0
C................................................. .....

But when I run the following program the solution is not satisfcatory. Could you please give me a hint of what is goint wrong?

INTEGER N
PARAMETER(N=3)
INTEGER i,j
REAL h,x,y(N),dydx(N),yout(N)
EXTERNAL derivs
x=0.0
y(1)=3
y(2)=15
y(3)=1

call derivs(x,y,dydx)
WRITE(*,*) ' time X Y Z '
WRITE(*,*) '================================================= ===='
do 11 i=1,100
h=0.0001*i
call rk4(y,dydx,N,x,h,yout,derivs)

write(*,*) h , yout(1), yout(2), yout(3)


11 continue
END

SUBROUTINE derivs(x,y,dydx)
REAL x,y(*),dydx(*)
dydx(1)= 10.0 * ( y(2) - y(1) )
dydx(2)= y(1) *( 28.0 - y(3) ) - y(2)
dydx(3)= y(1) * y(2) - (2.666666666*y(3))
return
END


Thanks

Makris

I did just that using Gill's method, a variation on 4th order Runge-Kutta, but using a macro in Excel. So you should be able to do the same with the NR subroutine. Your problem seems to be that your step size is too small and so is your number of steps. Try h = 0.01 and maybe 2100 steps. You might also want to use double precision (REAL*8), but do try single precision first. Even zipped, the spreadsheet is too big for this Forum, so send me your email address if you would like a copy (TO: turtle_hollow@sigmaxi.net). NB: Underscore, not space, between e_h.
Ciao, Foster Morrison