binaries/src/disembl/Tisean_3.0.1/source_f/upo.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   locate unstable periodic points
  23 c   author T. Schreiber (1998)
  24 c===========================================================================
  25       parameter(nx=1000000,mper=20)
  26       dimension x(nx)
  27       character*72 file, fout
  28       common /period/ x, nmax, m, eps
  29       data frac/0./, iper/1/, teq/-1./, tdis/-1./, tacc/-1./, h/-1./
  30       data iverb/1/
  31
  32       call whatido("locate unstable periodic points",iverb)
  33       m=max(imust("m"),1)
  34       eps=fcan("r",0.)
  35       frac=fcan("v",frac)
  36       teq=fcan("w",teq)
  37       tdis=fcan("W",tdis)
  38       h=fcan("s",h)
  39       tacc=fcan("a",tacc)
  40       iper=ican("p",iper)
  41       nmaxx=ican("l",nx)
  42       nexcl=ican("x",0)
  43       jcol=ican("c",0)
  44       icen=ican("n",nmaxx)
  45       isout=igetout(fout,iverb)
  46       if(eps.eq.0.and.frac.eq.0.) call usage()
  47
  48       do 10 ifi=1,nstrings()
  49          call nthstring(ifi,file)
  50          nmax=nmaxx
  51          call readfile(nmax,x,nexcl,jcol,file,iverb)
  52          if(file.eq."-") file="stdin"
  53          if(isout.eq.1) call addsuff(fout,file,"_upo_")
  54          call rms(nmax,x,sc,sd)
  55          if(frac.gt.0) eps=sd*frac
  56          if(teq.lt.0.) teq=eps
  57          if(tdis.lt.0.) tdis=eps
  58          if(tacc.lt.0.) tacc=eps
  59          if(h.lt.0.) h=eps
  60          if(isout.eq.1)
  61      .      write(fout(index(fout,"_upo_")+5:72),'(i2.2)') iper
  62          call outfile(fout,iunit,iverb)
  63          call findupo(iper,icen,teq,tdis,tacc,h,iunit,iverb)
  64  10      if(iunit.ne.istdout()) close(iunit)
  65       end
  66
  67       subroutine usage()
  68 c usage message
  69       call whatineed(
  70      .   "-m# [-r# | -v#] [-p# -w# -W# -a# -s# -n#"//
  71      .   " -o outfile -l# -x# -c# -V# -h] file(s)")
  72       call ptext("either -r or -v must be present")
  73       call popt("m","embedding dimension")
  74       call popt("r","absolute kernel bandwidth")
  75       call popt("v","same as fraction of standard deviation")
  76       call popt("p","period of orbit (1)")
  77       call popt("w","minimal separation of trial points (e)")
  78       call popt("W","minimal separation of distinct orbits (e)")
  79       call popt("a",
  80      .   "maximal error of orbit to be plotted (all plotted)")
  81       call popt("s","initial separation for stability (e)")
  82       call popt("n","number of trials (all points)")
  83       call popt("l","number of values to be read (all)")
  84       call popt("x","number of values to be skipped (0)")
  85       call popt("c","column to be read (1 or file,#)")
  86       call pout("file_upo_pp")
  87       call pall()
  88       call ptext("Verbosity levels (add what you want):")
  89       call ptext("          1 = input/output" )
  90       call ptext("          2 = print orbits found")
  91       call ptext("          4 = status after 1000 points")
  92       call ptext("          8 = status after 100 points")
  93       call ptext("         16 = status after 10 points")
  94       stop
  95       end
  96
  97       subroutine findupo(iper,icen,teq,tdis,tacc,h,iunit,iverb)
  98       parameter(nx=1000000,mper=20)
  99       external peri
 100       dimension x(nx),xp(mper),fvec(mper),xor(mper,nx),
 101      .   iw(mper),w0(mper,mper),w1(mper),w2(mper),w3(mper),
 102      .   w4(mper),w5(mper),w6(mper)
 103       common /period/ x, nmax, m, eps
 104
 105       if(iper.gt.mper) stop "findupo: make mper larger."
 106       tol=sqrt(r1mach(4))
 107       itry=0
 108       ior=0
 109       do 10 n=iper,nmax
 110          if(iv_10(iverb).eq.1) then
 111             if(mod(n,10).eq.0) write(istderr(),'(i7)') n
 112          else if(iv_100(iverb).eq.1) then
 113             if(mod(n,100).eq.0) write(istderr(),'(i7)') n
 114          else if(iv_1000(iverb).eq.1) then
 115             if(mod(n,1000).eq.0) write(istderr(),'(i7)') n
 116          endif
 117          if(known(n,iper,teq).eq.1) goto 10
 118          itry=itry+1
 119          if(itry.gt.icen) return
 120          do 20 i=1,iper
 121  20         xp(i)=x(n-iper+i)
 122          call snls1(peri,1,iper,iper,xp,fvec,w0,mper,tol,tol,0.,
 123      .      20*(iper+1),0.,w1,1,100.,0,info,nfev,ndum,iw,w2,w3,w4,w5,w6)
 124          err=enorm(iper,fvec)
 125          if(info.eq.-1.or.info.eq.5.or.err.gt.tacc) goto 10   ! unsuccessfull
 126          if(isold(iper,xp,ior,xor,tdis).eq.1) goto 10         ! already found
 127          ior=ior+1                                            ! a new orbit
 128          do 30 i=1,iper
 129  30         xor(i,ior)=xp(i)
 130          ipor=iperiod(iper,xp,tdis)
 131          sor=ipor*stab(iper,xp,h)/real(iper)
 132          call print(iper,xp,ipor,sor,err,iunit,iverb)
 133  10      continue
 134       end
 135
 136       function known(n,iper,tol)
 137 c return 1 if equivalent starting point has been tried
 138       parameter(nx=1000000)
 139       dimension x(nx)
 140       common /period/ x, nmax, m, eps
 141
 142       known=1
 143       do 10 nn=iper,n-1
 144          dis=0
 145          do 20 i=1,iper
 146  20         dis=dis+(x(n-iper+i)-x(nn-iper+i))**2
 147  10      if(sqrt(dis).lt.tol) return
 148       known=0
 149       end
 150
 151       function isold(iper,xp,ior,xor,toler)
 152 c determine if orbit is in data base
 153       parameter(mper=20)
 154       dimension xp(iper), xor(mper,*)
 155
 156       isold=1
 157       do 10 ip=1,iper
 158          do 20 io=1,ior
 159             dor=0
 160             do 30 i=1,iper
 161  30            dor=dor+(xp(i)-xor(i,io))**2
 162  20            if(sqrt(dor).le.toler) return
 163  10      call oshift(iper,xp)
 164       isold=0
 165       end
 166
 167       subroutine oshift(iper,xp)
 168 c leftshift orbit circularly by one position
 169       dimension xp(*)
 170
 171       h=xp(1)
 172       do 10 i=1,iper-1
 173  10      xp(i)=xp(i+1)
 174       xp(iper)=h
 175       end
 176
 177       function iperiod(iper,xp,tol)
 178 c determine shortest subperiod
 179       dimension xp(*)
 180
 181       do 10 iperiod=1,iper
 182          dis=0
 183          do 20 i=1,iper
 184             il=i-iperiod
 185             if(il.le.0) il=il+iper
 186  20         dis=dis+(xp(i)-xp(il))**2
 187  10      if(sqrt(dis).le.tol) return
 188       end
 189
 190       subroutine peri(iflag,mf,iper,xp,fvec,fjac,ldfjac)
 191 c built discrepancy vector (as called by snls1)
 192       dimension xp(*),fvec(*)
 193
 194       do 10 ip=1,iper
 195          fvec(ip)=xp(1)-fc(iper,xp,iflag)
 196  10      call oshift(iper,xp)
 197       end
 198
 199       function fc(iper,xp,iflag)
 200 c predict (cyclic) point 1, using iper,iper-1...
 201       parameter(nx=1000000)
 202       dimension  xp(*), x(nx)
 203       common /period/ x, nmax, m, eps
 204       data cut/20/
 205
 206       eps2=1./(2*eps*eps)
 207       ft=0
 208       sw=0
 209       fc=0
 210       do 10 n=m+1,nmax
 211          dis=0
 212          do 20 i=1,m
 213  20         dis=dis+(x(n-i)-xp(mod(m*iper-i,iper)+1))**2
 214          ddis=dis*eps2
 215          w=0
 216          if(ddis.lt.cut) w=exp(-ddis)
 217          ft=ft+w*x(n)
 218  10      sw=sw+w
 219       iflag=-1
 220       if(sw.eq.0) return   ! fc undefined, stop minimising
 221       fc=ft/sw
 222       iflag=1
 223       end
 224
 225       function stab(ilen,xp,h)
 226 c compute cycle stability by iteration of a tiny perturbation
 227       parameter(nx=1000000,mper=20,maxit=1000)
 228       dimension xp(*), x(nx), xcop(mper)
 229       common /period/ x, nmax, m, eps
 230
 231       if(mper.lt.ilen) stop "stability: make mper larger."
 232       iflag=1
 233       stab=0
 234       do 10 i=2,m
 235  10      xcop(i)=xp(mod(i-1,ilen)+1)
 236       xcop(1)=xp(1)+h
 237       do 20 it=1,maxit
 238          do 30 itt=1,ilen
 239             xx=fc(m,xcop,iflag)
 240             if(iflag.eq.-1) goto 1
 241             call oshift(m,xcop)
 242  30         xcop(m)=xx
 243          dis=0
 244          do 40 i=1,m
 245  40         dis=dis+(xcop(i)-xp(mod(i-1,ilen)+1))**2
 246          dis=sqrt(dis)
 247          stab=stab+log(dis/h)
 248          do 20 i=1,m
 249  20         xcop(i)=xp(mod(i-1,ilen)+1)*(1-h/dis) + xcop(i)*h/dis
 250  1    stab=stab/max(it-1,1)
 251       end
 252
 253       subroutine print(iper,xp,ipor,sor,err,iunit,iverb)
 254 c write orbit to iunit and to stdout
 255       dimension xp(*)
 256
 257       write(iunit,*)
 258       write(iunit,*) "period / accuracy / stability"
 259       write(iunit,*) ipor, err, exp(sor)
 260       do 10 i=1,ipor
 261  10      write(iunit,*) i, xp(i)
 262       if(iv_upo(iverb).eq.0) return
 263       write(istderr(),*)
 264       write(istderr(),*) "period / accuracy / stability"
 265       write(istderr(),*) ipor, err, exp(sor)
 266       do 20 i=1,ipor
 267  20      write(istderr(),*) i, xp(i)
 268       end