1 c===========================================================================
3 c This file is part of TISEAN
5 c Copyright (c) 1998-2007 Rainer Hegger, Holger Kantz, Thomas Schreiber
7 c TISEAN is free software; you can redistribute it and/or modify
8 c it under the terms of the GNU General Public License as published by
9 c the Free Software Foundation; either version 2 of the License, or
10 c (at your option) any later version.
12 c TISEAN is distributed in the hope that it will be useful,
13 c but WITHOUT ANY WARRANTY; without even the implied warranty of
14 c MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 c GNU General Public License for more details.
17 c You should have received a copy of the GNU General Public License
18 c along with TISEAN; if not, write to the Free Software
19 c Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
21 c===========================================================================
23 c integrates the Lorenz system with Runge Kutta fourth order
24 c author: H. Kantz (2007) based on earlier versions
26 c===========================================================================
29 real*8 x(3),u(3,3),sliap(3),bb,ss,rr,r1,r2,dh,s
34 call whatido("integration of the Lorenz system",iverb)
39 bb=fcan('B',2.666666667)
44 isout=igetout(fout,iverb)
46 if(isout.eq.1) fout="lorenz.dat"
47 call outfile(fout,iunit,iverb)
49 cc intermittency parameters
61 c step width of Runge Kutta integration dh:
63 c time intervals between re-orthogonalization of tangent space
64 c vectors: 0.01 units of time.
66 c length of transient in iteration steps:
67 itrans=real(itrans)/dh
68 totaltime=real(irun)/real(isamp)
72 . write(istderr(),*)'Lorenz trajectory covering',totaltime,
75 c x(1)=sqrt(s*(r+1.d0))+2.
93 call RUKU(3,x,u,dh,bb,ss,rr)
95 if (mod(i2,ireno).eq.0) then
107 write(iunit,101)x(1),x(2),x(3)
110 do 20 i2=1,irun*istep
113 call gauss(r1,r2,iseed1,iseed2)
116 call gauss(r1,r2,iseed1,iseed2)
122 call RUKU(3,x,u,dh,bb,ss,rr)
123 if (mod(i2,istep).eq.0) write(iunit,101)x(1),x(2),x(3)
124 if (mod(i2,ireno).eq.0) then
125 c Gram Schmidt Orthonormierung
127 sliap(1)=sliap(1)+log(s)
133 sliap(i)=sliap(i)+log(s)
143 rlevel=sqrt(rsq)/sqrt(xsq-xav*xav)*100.
145 . print*,'noise level in percent of x-coordinate',rlevel
149 write(istderr(),*)'Lyapunov exponents [1/unit time]'
151 write(istderr(),*)real(sliap(i)/totaltime)
155 101 format(2x,3f10.3)
160 subroutine FORCE(x,ff,dh,bb,ss,rr)
161 real*8 x(3),ff(3),dh,bb,ss,rr
163 ff(1)=dh*ss*(x(2)-x(1))
164 ff(2)=dh*(x(1)*(-x(3)+rr)-x(2))
165 ff(3)=dh*(x(1)*x(2)-bb*x(3))
170 subroutine LFORCE(x,u,fl,dh,bb,ss,rr)
171 real*8 x(3),u(3,3),dh,fl(3,3),bb,ss,rr
174 fl(j,1)=dh*ss*(u(j,2)-u(j,1))
175 fl(j,2)=dh*(u(j,1)*(rr-x(3))-x(1)*u(j,3)-u(j,2))
176 fl(j,3)=dh*(u(j,1)*x(2)+x(1)*u(j,2)-bb*u(j,3))
181 subroutine RUKU(n,x,u,dh,bb,ss,rr)
182 c 4th-order Runge Kutta
186 c add subroutine force
188 implicit real*8 (a-h,o-z)
189 real*8 x(3),ff1(3),ff2(3),ff3(3),ff4(3),dummy(3)
190 real*8 u(3,3),fl1(3,3),fl2(3,3),fl3(3,3),fl4(3,3)
193 call force(x,ff1,dh,bb,ss,rr)
194 call LFORCE(x,u,fl1,dh,bb,ss,rr)
197 dummy(i)=ff1(i)*.5d0+x(i)
199 dl(i,j)=fl1(i,j)*.5d0+u(i,j)
203 call force(dummy,ff2,dh,bb,ss,rr)
204 call LFORCE(dummy,dl,fl2,dh,bb,ss,rr)
207 dummy(i)=ff2(i)*.5d0+x(i)
209 dl(i,j)=fl2(i,j)*.5d0+u(i,j)
213 call force(dummy,ff3,dh,bb,ss,rr)
214 call LFORCE(dummy,dl,fl3,dh,bb,ss,rr)
219 dl(i,j)=fl3(i,j)+u(i,j)
223 call force(dummy,ff4,dh,bb,ss,rr)
224 call LFORCE(dummy,dl,fl4,dh,bb,ss,rr)
227 x(i)=x(i)+ff1(i)/6.d0+ff2(i)/3.d0+ff3(i)/3.d0+ff4(i)/6.d0
229 u(i,j)=u(i,j)+fl1(i,j)/6.d0+fl2(i,j)/3.d0+fl3(i,j)/3.d0
237 subroutine NORM(u,i,s)
250 subroutine PROJ(u,i,j)
256 20 u(i,k)=u(i,k)-s*u(j,k)
260 c>-------------------------------------------------------
261 subroutine gauss(r1,r2,iseed1,iseed2)
266 call RANDOM1(p,iseed1)
267 call RANDOM1(phi,iseed2)
270 r=sqrt(-log(1.d0-p)*2.d0)
276 c>-------------------------------------------------------
277 subroutine RANDOM1(r,iseed)
279 c random number generator of Park & Miller
280 c random numbers in [0,1] !!!
287 r=dfloat(ix)/dfloat(im)
291 c>------------------------------------------------------------------
296 . "-l# [-f# -r# -R# -S# -B# -o outfile -x# -V# -h]")
297 call popt("l","length of trajectory x,y,z")
298 call popt("f","sample points per unit time [100]")
299 call popt("r","absolute noise amplitute [0]")
300 call popt("R","parameter r [28]")
301 call popt("S","parameter sigma [10]")
302 call popt("B","parameter b [8/3]")
303 call popt("x","transient discarded [100 units of time]")
304 c call popt("L","if present: compute Lyapunov exponents")
305 call pout("lorenz.dat")