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===========================================================================
22 c locate unstable periodic points
23 c author T. Schreiber (1998)
24 c===========================================================================
25 parameter(nx=1000000,mper=20)
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./
32 call whatido("locate unstable periodic points",iverb)
45 isout=igetout(fout,iverb)
46 if(eps.eq.0.and.frac.eq.0.) call usage()
48 do 10 ifi=1,nstrings()
49 call nthstring(ifi,file)
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
57 if(tdis.lt.0.) tdis=eps
58 if(tacc.lt.0.) tacc=eps
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)
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)")
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")
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")
97 subroutine findupo(iper,icen,teq,tdis,tacc,h,iunit,iverb)
98 parameter(nx=1000000,mper=20)
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
105 if(iper.gt.mper) stop "findupo: make mper larger."
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
117 if(known(n,iper,teq).eq.1) goto 10
119 if(itry.gt.icen) return
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)
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
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)
136 function known(n,iper,tol)
137 c return 1 if equivalent starting point has been tried
138 parameter(nx=1000000)
140 common /period/ x, nmax, m, eps
146 20 dis=dis+(x(n-iper+i)-x(nn-iper+i))**2
147 10 if(sqrt(dis).lt.tol) return
151 function isold(iper,xp,ior,xor,toler)
152 c determine if orbit is in data base
154 dimension xp(iper), xor(mper,*)
161 30 dor=dor+(xp(i)-xor(i,io))**2
162 20 if(sqrt(dor).le.toler) return
163 10 call oshift(iper,xp)
167 subroutine oshift(iper,xp)
168 c leftshift orbit circularly by one position
177 function iperiod(iper,xp,tol)
178 c determine shortest subperiod
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
190 subroutine peri(iflag,mf,iper,xp,fvec,fjac,ldfjac)
191 c built discrepancy vector (as called by snls1)
192 dimension xp(*),fvec(*)
195 fvec(ip)=xp(1)-fc(iper,xp,iflag)
196 10 call oshift(iper,xp)
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
213 20 dis=dis+(x(n-i)-xp(mod(m*iper-i,iper)+1))**2
216 if(ddis.lt.cut) w=exp(-ddis)
220 if(sw.eq.0) return ! fc undefined, stop minimising
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
231 if(mper.lt.ilen) stop "stability: make mper larger."
235 10 xcop(i)=xp(mod(i-1,ilen)+1)
240 if(iflag.eq.-1) goto 1
245 40 dis=dis+(xcop(i)-xp(mod(i-1,ilen)+1))**2
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)
253 subroutine print(iper,xp,ipor,sor,err,iunit,iverb)
254 c write orbit to iunit and to stdout
258 write(iunit,*) "period / accuracy / stability"
259 write(iunit,*) ipor, err, exp(sor)
261 10 write(iunit,*) i, xp(i)
262 if(iv_upo(iverb).eq.0) return
264 write(istderr(),*) "period / accuracy / stability"
265 write(istderr(),*) ipor, err, exp(sor)
267 20 write(istderr(),*) i, xp(i)