Ask Question

Name:
Title:
Your Question:

Answer Question

Name:
Your Answer:
User Submitted Source Code!


Description:
  rossler
Language: C/C++
Code:
c===========================================================================
c
c   This file is part of TISEAN

c   Copyright (c) 1998-2007 Rainer Hegger, Holger Kantz, Thomas Schreiber

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   rossler.f
c   integrates the rossler 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,cc,aa,r1,r2,dh,s
      character*72 fout
      data iverb/1/

      iverb=ican('V',iverb)
      call whatido("integration of the rossler system",iverb)
      irun=imust('l')
      itrans=ican('x',100)
      aa=fcan('A',0.2)
      cc=fcan('C',5.7)
      bb=fcan('B',0.2)
      isamp=ican('f',100)
      sn=fcan('r',0.)
c      ilyap=lopt('L',1)

      isout=igetout(fout,iverb)

      if(isout.eq.1) fout="rossler.dat"
      call outfile(fout,iunit,iverb)


      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(),*)'rossler 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)=0.1.
      x(2)=0.1.
      x(3)=0.

      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,cc,aa)

        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,cc,aa)
        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,cc,aa)
      real*8 x(3),ff(3),dh,bb,cc,aa

        ff(1)=dh*(x(2)-x(3))
        ff(2)=dh*(x(1)+aa*x(2))
        ff(3)=dh*(bb+x(3)*(x(1)-cc))

      return
      end

      subroutine LFORCE(x,u,fl,dh,bb,cc,aa)
      real*8 x(3),u(3,3),dh,fl(3,3),bb,cc,aa

       do j=1,3
         fl(j,1)=dh*cc*(u(j,2)-u(j,3))
         fl(j,2)=dh*(u(j,1)+aa*u(j,2))
         fl(j,3)=dh*(bb+u(j,3)*(u(j,1)-cc))
       enddo
      return
      end

      subroutine RUKU(n,x,u,dh,bb,cc,aa)
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,cc,aa)
      call LFORCE(x,u,fl1,dh,bb,cc,aa)

      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,cc,aa)
      call LFORCE(dummy,dl,fl2,dh,bb,cc,aa)

      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,cc,aa)
      call LFORCE(dummy,dl,fl3,dh,bb,cc,aa)

      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,cc,aa)
      call LFORCE(dummy,dl,fl4,dh,bb,cc,aa)

      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("rossler.dat")
      call pall()
      stop
      end



Comments: