website/archive/binaries/mac/src/disembl/Tisean_3.0.1/source_f/lorenz.f

   1 c===========================================================================
   2 c
   3 c   This file is part of TISEAN
   4 c
   5 c   Copyright (c) 1998-2007 Rainer Hegger, Holger Kantz, Thomas Schreiber
   6 c
   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.
  11 c
  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.
  16 c
  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
  20 c
  21 c===========================================================================
  22 c   lorenz.f
  23 c   integrates the Lorenz system with Runge Kutta fourth order
  24 c   author: H. Kantz (2007) based on earlier versions
  25 c   with optional noise
  26 c===========================================================================
  27 c
  28
  29       real*8 x(3),u(3,3),sliap(3),bb,ss,rr,r1,r2,dh,s
  30       character*72 fout
  31       data iverb/1/
  32
  33       iverb=ican('V',iverb)
  34       call whatido("integration of the Lorenz system",iverb)
  35       irun=imust('l')
  36       itrans=ican('x',100)
  37       rr=fcan('R',28.0)
  38       ss=fcan('S',10.0)
  39       bb=fcan('B',2.666666667)
  40       isamp=ican('f',100)
  41       sn=fcan('r',0.)
  42 c      ilyap=lopt('L',1)
  43
  44       isout=igetout(fout,iverb)
  45
  46       if(isout.eq.1) fout="lorenz.dat"
  47       call outfile(fout,iunit,iverb)
  48
  49 cc intermittency parameters
  50 c       ss=10.d0
  51 c       rr=166.11d0
  52 c       bb=8.d0/3.d0
  53
  54       iseed1=6456423
  55       iseed2=7243431
  56
  57       xav=0.
  58       xsq=0.
  59       rsq=0.
  60
  61 c     step width of Runge Kutta integration dh:
  62       dh=.0005d0
  63 c     time intervals between re-orthogonalization of tangent space
  64 c            vectors: 0.01 units of time.
  65       ireno=.01d0/dh
  66 c     length of transient in iteration steps:
  67       itrans=real(itrans)/dh
  68       totaltime=real(irun)/real(isamp)
  69       istep=1.d0/dh/isamp
  70
  71       if (iverb.eq.1)
  72      . write(istderr(),*)'Lorenz trajectory covering',totaltime,
  73      .                  ' time units'
  74
  75 c      x(1)=sqrt(s*(r+1.d0))+2.
  76 c      x(2)=x(1)-1.d0
  77 c      x(3)=r
  78
  79       x(1)=5.
  80       x(2)=-10.
  81       x(3)=3.
  82
  83       do 1 i=1,3
  84        sliap(i)=0.d0
  85        do j=1,3
  86         u(i,j)=0.d0
  87        enddo
  88        u(i,i)=1.d0
  89 1     continue
  90
  91       do 10 i2=1,itrans
  92
  93         call RUKU(3,x,u,dh,bb,ss,rr)
  94
  95         if (mod(i2,ireno).eq.0) then
  96           call norm(u,1,s)
  97           do i=2,3
  98             do j=1,i-1
  99               call proj(u,i,j)
 100             enddo
 101            call NORM(u,i,s)
 102           enddo
 103         endif
 104
 105 10    continue
 106
 107       write(iunit,101)x(1),x(2),x(3)
 108
 109  100  continue
 110       do 20 i2=1,irun*istep
 111 c       add noise
 112         if (sn.gt.0.0) then
 113           call gauss(r1,r2,iseed1,iseed2)
 114           x(1)=x(1)+r1*sn
 115           x(2)=x(2)+r2*sn
 116           call gauss(r1,r2,iseed1,iseed2)
 117           x(3)=x(3)+r1*sn
 118           xav=xav+x(1)
 119           xsq=xsq+x(1)**2
 120           rsq=rsq+r1*r1
 121         endif
 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
 126           call norm(u,1,s)
 127           sliap(1)=sliap(1)+log(s)
 128           do i=2,3
 129             do j=1,i-1
 130              call proj(u,i,j)
 131             enddo
 132             call NORM(u,i,s)
 133             sliap(i)=sliap(i)+log(s)
 134           enddo
 135         endif
 136
 137  20   continue
 138
 139       if (sn.gt.0.0) then
 140         xav=xav/irun/istep
 141         xsq=xsq/irun/istep
 142         rsq=rsq/irun/istep
 143         rlevel=sqrt(rsq)/sqrt(xsq-xav*xav)*100.
 144         if (iverb.eq.1)
 145      .   print*,'noise level in percent of x-coordinate',rlevel
 146       endif
 147       if (iverb.eq.1) then
 148        write(istderr(),*)
 149        write(istderr(),*)'Lyapunov exponents [1/unit time]'
 150        do i=1,3
 151         write(istderr(),*)real(sliap(i)/totaltime)
 152        enddo
 153       endif
 154
 155  101  format(2x,3f10.3)
 156
 157       stop
 158       end
 159
 160       subroutine FORCE(x,ff,dh,bb,ss,rr)
 161       real*8 x(3),ff(3),dh,bb,ss,rr
 162
 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))
 166
 167       return
 168       end
 169
 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
 172
 173        do j=1,3
 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))
 177        enddo
 178       return
 179       end
 180
 181       subroutine RUKU(n,x,u,dh,bb,ss,rr)
 182 c     4th-order Runge Kutta
 183 c     initial point x
 184 c     final point y
 185 c     stepsize dh
 186 c     add subroutine force
 187
 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)
 191       real*8 dl(3,3)
 192
 193       call force(x,ff1,dh,bb,ss,rr)
 194       call LFORCE(x,u,fl1,dh,bb,ss,rr)
 195
 196       do i=1,n
 197       dummy(i)=ff1(i)*.5d0+x(i)
 198         do j=1,3
 199         dl(i,j)=fl1(i,j)*.5d0+u(i,j)
 200         enddo
 201       enddo
 202
 203       call force(dummy,ff2,dh,bb,ss,rr)
 204       call LFORCE(dummy,dl,fl2,dh,bb,ss,rr)
 205
 206       do i=1,n
 207       dummy(i)=ff2(i)*.5d0+x(i)
 208         do j=1,3
 209         dl(i,j)=fl2(i,j)*.5d0+u(i,j)
 210         enddo
 211       enddo
 212
 213       call force(dummy,ff3,dh,bb,ss,rr)
 214       call LFORCE(dummy,dl,fl3,dh,bb,ss,rr)
 215
 216       do i=1,n
 217       dummy(i)=ff3(i)+x(i)
 218         do j=1,3
 219         dl(i,j)=fl3(i,j)+u(i,j)
 220         enddo
 221       enddo
 222
 223       call force(dummy,ff4,dh,bb,ss,rr)
 224       call LFORCE(dummy,dl,fl4,dh,bb,ss,rr)
 225
 226       do i=1,n
 227       x(i)=x(i)+ff1(i)/6.d0+ff2(i)/3.d0+ff3(i)/3.d0+ff4(i)/6.d0
 228         do j=1,3
 229         u(i,j)=u(i,j)+fl1(i,j)/6.d0+fl2(i,j)/3.d0+fl3(i,j)/3.d0
 230      +               +fl4(i,j)/6.d0
 231         enddo
 232       enddo
 233
 234       return
 235       end
 236
 237       subroutine NORM(u,i,s)
 238       real*8 u(3,3),s
 239
 240       s=0.d0
 241       do 10 j=1,3
 242 10    s=s+u(i,j)**2
 243       s=sqrt(s)
 244       si=1.d0/s
 245       do 20 j=1,3
 246 20    u(i,j)=u(i,j)*si
 247       return
 248       end
 249
 250       subroutine PROJ(u,i,j)
 251       real*8 u(3,3),s
 252       s=0.d0
 253       do 10 k=1,3
 254 10      s=s+u(i,k)*u(j,k)
 255       do 20 k=1,3
 256 20      u(i,k)=u(i,k)-s*u(j,k)
 257       return
 258       end
 259
 260 c>-------------------------------------------------------
 261       subroutine gauss(r1,r2,iseed1,iseed2)
 262
 263       real*8 r1,r2,p,phi,r
 264       pii=8.d0*atan(1.d0)
 265
 266       call RANDOM1(p,iseed1)
 267       call RANDOM1(phi,iseed2)
 268
 269       phi=phi*pii
 270       r=sqrt(-log(1.d0-p)*2.d0)
 271
 272       r1=r*sin(phi)
 273       r2=r*cos(phi)
 274       return
 275       end
 276 c>-------------------------------------------------------
 277       subroutine RANDOM1(r,iseed)
 278 c
 279 c     random number generator of Park & Miller
 280 c     random numbers in [0,1] !!!
 281       real*8 r
 282       integer*8 ia,im,ix
 283       ia=7**5
 284       im=2147483647
 285       ix=iseed
 286       ix=mod(ia*ix,im)
 287       r=dfloat(ix)/dfloat(im)
 288       iseed=ix
 289       return
 290       end
 291 c>------------------------------------------------------------------
 292       subroutine usage()
 293 c usage message
 294
 295       call whatineed(
 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")
 306       call pall()
 307       stop
 308       end
 309
 310
 311