+++ /dev/null
-c===========================================================================
-c
-c This file is part of TISEAN
-c
-c Copyright (c) 1998-2007 Rainer Hegger, Holger Kantz, Thomas Schreiber
-c
-c TISEAN is free software; you can redistribute it and/or modify
-c it under the terms of the GNU General Public License as published by
-c the Free Software Foundation; either version 2 of the License, or
-c (at your option) any later version.
-c
-c TISEAN is distributed in the hope that it will be useful,
-c but WITHOUT ANY WARRANTY; without even the implied warranty of
-c MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-c GNU General Public License for more details.
-c
-c You should have received a copy of the GNU General Public License
-c along with TISEAN; if not, write to the Free Software
-c Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-c
-c===========================================================================
-c lorenz.f
-c integrates the Lorenz system with Runge Kutta fourth order
-c author: H. Kantz (2007) based on earlier versions
-c with optional noise
-c===========================================================================
-c
-
- real*8 x(3),u(3,3),sliap(3),bb,ss,rr,r1,r2,dh,s
- character*72 fout
- data iverb/1/
-
- iverb=ican('V',iverb)
- call whatido("integration of the Lorenz system",iverb)
- irun=imust('l')
- itrans=ican('x',100)
- rr=fcan('R',28.0)
- ss=fcan('S',10.0)
- bb=fcan('B',2.666666667)
- isamp=ican('f',100)
- sn=fcan('r',0.)
-c ilyap=lopt('L',1)
-
- isout=igetout(fout,iverb)
-
- if(isout.eq.1) fout="lorenz.dat"
- call outfile(fout,iunit,iverb)
-
-cc intermittency parameters
-c ss=10.d0
-c rr=166.11d0
-c bb=8.d0/3.d0
-
- iseed1=6456423
- iseed2=7243431
-
- xav=0.
- xsq=0.
- rsq=0.
-
-c step width of Runge Kutta integration dh:
- dh=.0005d0
-c time intervals between re-orthogonalization of tangent space
-c vectors: 0.01 units of time.
- ireno=.01d0/dh
-c length of transient in iteration steps:
- itrans=real(itrans)/dh
- totaltime=real(irun)/real(isamp)
- istep=1.d0/dh/isamp
-
- if (iverb.eq.1)
- . write(istderr(),*)'Lorenz trajectory covering',totaltime,
- . ' time units'
-
-c x(1)=sqrt(s*(r+1.d0))+2.
-c x(2)=x(1)-1.d0
-c x(3)=r
-
- x(1)=5.
- x(2)=-10.
- x(3)=3.
-
- do 1 i=1,3
- sliap(i)=0.d0
- do j=1,3
- u(i,j)=0.d0
- enddo
- u(i,i)=1.d0
-1 continue
-
- do 10 i2=1,itrans
-
- call RUKU(3,x,u,dh,bb,ss,rr)
-
- if (mod(i2,ireno).eq.0) then
- call norm(u,1,s)
- do i=2,3
- do j=1,i-1
- call proj(u,i,j)
- enddo
- call NORM(u,i,s)
- enddo
- endif
-
-10 continue
-
- write(iunit,101)x(1),x(2),x(3)
-
- 100 continue
- do 20 i2=1,irun*istep
-c add noise
- if (sn.gt.0.0) then
- call gauss(r1,r2,iseed1,iseed2)
- x(1)=x(1)+r1*sn
- x(2)=x(2)+r2*sn
- call gauss(r1,r2,iseed1,iseed2)
- x(3)=x(3)+r1*sn
- xav=xav+x(1)
- xsq=xsq+x(1)**2
- rsq=rsq+r1*r1
- endif
- call RUKU(3,x,u,dh,bb,ss,rr)
- if (mod(i2,istep).eq.0) write(iunit,101)x(1),x(2),x(3)
- if (mod(i2,ireno).eq.0) then
-c Gram Schmidt Orthonormierung
- call norm(u,1,s)
- sliap(1)=sliap(1)+log(s)
- do i=2,3
- do j=1,i-1
- call proj(u,i,j)
- enddo
- call NORM(u,i,s)
- sliap(i)=sliap(i)+log(s)
- enddo
- endif
-
- 20 continue
-
- if (sn.gt.0.0) then
- xav=xav/irun/istep
- xsq=xsq/irun/istep
- rsq=rsq/irun/istep
- rlevel=sqrt(rsq)/sqrt(xsq-xav*xav)*100.
- if (iverb.eq.1)
- . print*,'noise level in percent of x-coordinate',rlevel
- endif
- if (iverb.eq.1) then
- write(istderr(),*)
- write(istderr(),*)'Lyapunov exponents [1/unit time]'
- do i=1,3
- write(istderr(),*)real(sliap(i)/totaltime)
- enddo
- endif
-
- 101 format(2x,3f10.3)
-
- stop
- end
-
- subroutine FORCE(x,ff,dh,bb,ss,rr)
- real*8 x(3),ff(3),dh,bb,ss,rr
-
- ff(1)=dh*ss*(x(2)-x(1))
- ff(2)=dh*(x(1)*(-x(3)+rr)-x(2))
- ff(3)=dh*(x(1)*x(2)-bb*x(3))
-
- return
- end
-
- subroutine LFORCE(x,u,fl,dh,bb,ss,rr)
- real*8 x(3),u(3,3),dh,fl(3,3),bb,ss,rr
-
- do j=1,3
- fl(j,1)=dh*ss*(u(j,2)-u(j,1))
- fl(j,2)=dh*(u(j,1)*(rr-x(3))-x(1)*u(j,3)-u(j,2))
- fl(j,3)=dh*(u(j,1)*x(2)+x(1)*u(j,2)-bb*u(j,3))
- enddo
- return
- end
-
- subroutine RUKU(n,x,u,dh,bb,ss,rr)
-c 4th-order Runge Kutta
-c initial point x
-c final point y
-c stepsize dh
-c add subroutine force
-
- implicit real*8 (a-h,o-z)
- real*8 x(3),ff1(3),ff2(3),ff3(3),ff4(3),dummy(3)
- real*8 u(3,3),fl1(3,3),fl2(3,3),fl3(3,3),fl4(3,3)
- real*8 dl(3,3)
-
- call force(x,ff1,dh,bb,ss,rr)
- call LFORCE(x,u,fl1,dh,bb,ss,rr)
-
- do i=1,n
- dummy(i)=ff1(i)*.5d0+x(i)
- do j=1,3
- dl(i,j)=fl1(i,j)*.5d0+u(i,j)
- enddo
- enddo
-
- call force(dummy,ff2,dh,bb,ss,rr)
- call LFORCE(dummy,dl,fl2,dh,bb,ss,rr)
-
- do i=1,n
- dummy(i)=ff2(i)*.5d0+x(i)
- do j=1,3
- dl(i,j)=fl2(i,j)*.5d0+u(i,j)
- enddo
- enddo
-
- call force(dummy,ff3,dh,bb,ss,rr)
- call LFORCE(dummy,dl,fl3,dh,bb,ss,rr)
-
- do i=1,n
- dummy(i)=ff3(i)+x(i)
- do j=1,3
- dl(i,j)=fl3(i,j)+u(i,j)
- enddo
- enddo
-
- call force(dummy,ff4,dh,bb,ss,rr)
- call LFORCE(dummy,dl,fl4,dh,bb,ss,rr)
-
- do i=1,n
- x(i)=x(i)+ff1(i)/6.d0+ff2(i)/3.d0+ff3(i)/3.d0+ff4(i)/6.d0
- do j=1,3
- u(i,j)=u(i,j)+fl1(i,j)/6.d0+fl2(i,j)/3.d0+fl3(i,j)/3.d0
- + +fl4(i,j)/6.d0
- enddo
- enddo
-
- return
- end
-
- subroutine NORM(u,i,s)
- real*8 u(3,3),s
-
- s=0.d0
- do 10 j=1,3
-10 s=s+u(i,j)**2
- s=sqrt(s)
- si=1.d0/s
- do 20 j=1,3
-20 u(i,j)=u(i,j)*si
- return
- end
-
- subroutine PROJ(u,i,j)
- real*8 u(3,3),s
- s=0.d0
- do 10 k=1,3
-10 s=s+u(i,k)*u(j,k)
- do 20 k=1,3
-20 u(i,k)=u(i,k)-s*u(j,k)
- return
- end
-
-c>-------------------------------------------------------
- subroutine gauss(r1,r2,iseed1,iseed2)
-
- real*8 r1,r2,p,phi,r
- pii=8.d0*atan(1.d0)
-
- call RANDOM1(p,iseed1)
- call RANDOM1(phi,iseed2)
-
- phi=phi*pii
- r=sqrt(-log(1.d0-p)*2.d0)
-
- r1=r*sin(phi)
- r2=r*cos(phi)
- return
- end
-c>-------------------------------------------------------
- subroutine RANDOM1(r,iseed)
-c
-c random number generator of Park & Miller
-c random numbers in [0,1] !!!
- real*8 r
- integer*8 ia,im,ix
- ia=7**5
- im=2147483647
- ix=iseed
- ix=mod(ia*ix,im)
- r=dfloat(ix)/dfloat(im)
- iseed=ix
- return
- end
-c>------------------------------------------------------------------
- subroutine usage()
-c usage message
-
- call whatineed(
- . "-l# [-f# -r# -R# -S# -B# -o outfile -x# -V# -h]")
- call popt("l","length of trajectory x,y,z")
- call popt("f","sample points per unit time [100]")
- call popt("r","absolute noise amplitute [0]")
- call popt("R","parameter r [28]")
- call popt("S","parameter sigma [10]")
- call popt("B","parameter b [8/3]")
- call popt("x","transient discarded [100 units of time]")
-c call popt("L","if present: compute Lyapunov exponents")
- call pout("lorenz.dat")
- call pall()
- stop
- end
-
-
-
git://source.jalview.org / jabaws.git / blobdiff