c======================================================================= c obitnt: This is a routine for finding the periodic orbits of a c driven system. It uses a shooting method that integrates the c equations of motion of a system over a period nT (1<= n<=nmax) c combined with a Newton iteration that successively adjusts the c initial conditions until the initial and final conditions match. c c This program is currently set up to find the single periodic c orbit for the driven harmonic oscillator: c c H = p**2/(2*m) + m*w0**2*x**2/2 + V*x*cos(w*t+theta) c c To modify the routine to find periodic orbits for other systems, c one must rewrite the fcn and jac routines in 'fcn.f' and alter c the common block in 'fcn.inc'. Also, any common block c parameters must have their values set in this routine. c======================================================================= program orbitnt implicit none c----------------------------------------------------------------------- c Communication with routines in 'fcn.f' c----------------------------------------------------------------------- include 'fcn.inc' c----------------------------------------------------------------------- c Argument parsing variables c----------------------------------------------------------------------- integer iargc real*8 r8arg, r8_never parameter (r8_never = -1.0d-60) c----------------------------------------------------------------------- c Size of the system (neq = # of odes, npars = # of initial vals.) c----------------------------------------------------------------------- integer neq parameter (neq = 7) integer npars parameter (npars = 2) c----------------------------------------------------------------------- c Tolerances for lsoda (ltol) and Newton convergence (ntol) c----------------------------------------------------------------------- real*8 ltol, default_ltol parameter (default_ltol = 1.0d-8) real*8 ntol, default_ntol parameter (default_ntol = 1.0d-8) c----------------------------------------------------------------------- c ivs: array containing initial values (will be altered) c J: Jacobian used in Newton step c res: residual vector in Newton step c y: solution to odes c----------------------------------------------------------------------- real*8 ivs(npars), J(npars,npars) real*8 res(npars) real*8 y(neq) c----------------------------------------------------------------------- c S: action along classical path c dx,dp: finite differences used in calculation of J above c EPS: sqaure root of estimated machine precision c----------------------------------------------------------------------- real*8 S real*8 dx, dp, EPS c----------------------------------------------------------------------- c ipiv: storage needed by dgesv c info: return code from dgesv c ieq, ipar, iter: indices for eqns, parameters, iterations c mxiter: maximum number of Newton iterations at each period c nrhs: number of right-hand-sides sent to dgesv = 1 c----------------------------------------------------------------------- integer ipiv(neq) integer info integer ieq, ipar, iter integer mxiter, nrhs parameter (mxiter = 10, nrhs = 1) c----------------------------------------------------------------------- c T0: period of driving field c iN: index for orbit periods (period = iN*T) c maxN: maximum period to be searched c----------------------------------------------------------------------- real*8 T0, pi integer iN, maxN parameter (maxN = 10) c----------------------------------------------------------------------- c x0: permanent storage for inital value of x (not altered) c p0: ditto for p c R: residue of the periodic orbit c----------------------------------------------------------------------- real*8 x0, p0 real*8 R c----------------------------------------------------------------------- c norms used in monitoring convergence of Newton iterations c----------------------------------------------------------------------- real*8 nrm2res, nrm2upd, nrm2ivs real*8 drelabs, dv12norm c----------------------------------------------------------------------- c lsoda variables c----------------------------------------------------------------------- external fcn, jac real*8 tbgn, tend integer itol real*8 rtol, atol integer itask, istate, iopt integer ml, mu parameter (ml = 1, mu = 2) integer lrw parameter (lrw = 20 + 16*neq) real*8 rwork(lrw) integer liw parameter (liw = 20 + neq) integer iwork(liw) integer jt c----------------------------------------------------------------------- c Logical variables for tracing options c----------------------------------------------------------------------- logical ltrace parameter (ltrace = .true.) logical otrace parameter (otrace = .true.) c----------------------------------------------------------------------- c ofile: file into which orbit coords. will be written c----------------------------------------------------------------------- character*(*) ofile parameter (ofile = 'orbits') integer ufrom, getu c----------------------------------------------------------------------- c notpp: number of times orbit coords will be written per period c maxout: maximum number of output times c tout: array containing output times c----------------------------------------------------------------------- integer notpp parameter (notpp = 100) integer maxout, itout parameter (maxout = notpp*maxN) real*8 tout(maxout) c----------------------------------------------------------------------- c Parse command line arguments c----------------------------------------------------------------------- if (iargc() .lt. npars) go to 900 do ipar = 1, npars ivs(ipar) = r8arg(ipar, r8_never) if (ivs(ipar) .eq. r8_never) go to 900 end do ltol = r8arg(npars+1, default_ltol) ntol = r8arg(npars+2, default_ntol) if (ltol .le. 0.0d0) ltol = default_ltol if (ntol .le. 0.0d0) ntol = default_ntol pi = 4.0d0*atan(1.0d0) c----------------------------------------------------------------------- c Set values for common block parameters c w: frequency of driving field c theta: initial phase of driving field c V: strength of driving field c m: mass of particle c W: harmonic oscillator frequency c----------------------------------------------------------------------- w = 4.0d0*pi theta = 0.0d0 m = 1.0d0 w0 = 4.0d0 V = 1.0d-1 c----------------------------------------------------------------------- c Set values of x0, p0, T, EPS c----------------------------------------------------------------------- x0 = ivs(1) p0 = ivs(2) T0 = 2.0d0*pi/w EPS = sqrt(ltol) c----------------------------------------------------------------------- c Tracing of parameters c----------------------------------------------------------------------- if (ltrace) then write(0,*) 'w = ',w,', theta = ',theta,', m = ',m write(0,*) 'w0 = ',w0,', V = ',V,', T0 = ',T0 write(0,*) 'x0 = ',x0,', p0 = ',p0,', EPS = ',EPS write(0,*) 'ltol = ',ltol,', ntol = ',ntol end if c----------------------------------------------------------------------- c Set lsoda parameters c----------------------------------------------------------------------- itol = 1 rtol = ltol atol = ltol itask = 1 iopt = 0 jt = 4 iwork(1) = ml iwork(2) = mu c----------------------------------------------------------------------- c Open file ofile c----------------------------------------------------------------------- if (otrace) then ufrom = getu() open(ufrom,file=ofile,status='new') end if write(*,*) ' N R S' c----------------------------------------------------------------------- c Begin loop over periods NT c----------------------------------------------------------------------- do iN = 1, maxN write(0,*) 'For orbits N = ',iN write(0,*) write(0,*) ' Iter x0 p0 '// & ' log10(upd) log10(res)' write(0,*) c----------------------------------------------------------------------- c Reset initial values to original values c----------------------------------------------------------------------- ivs(1) = x0 ivs(2) = p0 c----------------------------------------------------------------------- c Begin Newton loop c----------------------------------------------------------------------- do iter = 1, mxiter dx = max(EPS*abs(ivs(1)),1d-8) dp = max(EPS*abs(ivs(2)),1d-8) y(1) = 0.0d0 y(2) = ivs(1) y(3) = ivs(2) y(4) = ivs(1) y(5) = ivs(2) + dp c----------------------------------------------------------------------- c Trick to reduce finite precision error (Num. Rec.) c----------------------------------------------------------------------- dp = y(5) - ivs(2) y(6) = ivs(1) + dx dx = y(6) - ivs(1) y(7) = ivs(2) c----------------------------------------------------------------------- c Define time interval and call lsoda c----------------------------------------------------------------------- tbgn = 0.0d0 tend = iN*T0 istate = 1 call lsoda(fcn,neq,y,tbgn,tend,itol,rtol,atol,itask, & istate,iopt,rwork,lrw,iwork,liw,jac,jt) c----------------------------------------------------------------------- c Error message if lsoda fails c----------------------------------------------------------------------- if (istate .ge. -1) go to 100 write(0,*) 'orbitnt: Error return ',istate, & ' from integrator lsoda.' if (otrace) then close(ufrom) end if stop 100 continue c----------------------------------------------------------------------- c If lsoda is successful, calculate residual c----------------------------------------------------------------------- res(1) = y(2) - ivs(1) res(2) = y(3) - ivs(2) nrm2res = dv12norm(res,npars) c----------------------------------------------------------------------- c Define Jacobian for this Newton step c----------------------------------------------------------------------- J(1,1) = (y(6)-y(2))/dx-1.0d0 J(1,2) = (y(4)-y(2))/dp J(2,1) = (y(7)-y(3))/dx J(2,2) = (y(5)-y(3))/dp-1.0d0 c----------------------------------------------------------------------- c call dgesv to solve Jdx=res c----------------------------------------------------------------------- call dgesv(npars,nrhs,J,npars,ipiv,res,npars,info) c----------------------------------------------------------------------- c Error message if dgesv fails c----------------------------------------------------------------------- if (info .ne. 0) then write(0,*) 'orbitnt: dgesv failed' if (otrace) then close(ufrom) end if stop else c----------------------------------------------------------------------- c if dgesv succeeds, calculate norms and update ivs c----------------------------------------------------------------------- nrm2ivs = dv12norm(ivs,npars) nrm2upd = dv12norm(res,npars) do ipar = 1, npars ivs(ipar) = ivs(ipar) - res(ipar) end do c----------------------------------------------------------------------- c Tracing output to monitor convergence c----------------------------------------------------------------------- write(0,200) iter, ivs(1), ivs(2), & log10(max(nrm2upd,1.0d-60)), & log10(max(nrm2res,1.0d-60)) 200 format(i2,1p,2e24.16,0p,2f8.2) c----------------------------------------------------------------------- c Check for convergence c----------------------------------------------------------------------- if (drelabs(nrm2upd,nrm2ivs,1.0d-6) .le. ntol) & go to 300 end if end do c----------------------------------------------------------------------- c No-convergence exit c----------------------------------------------------------------------- write(0,*) 'No convergence after ',mxiter, ' iterations'// & ' for orbit N = ',iN if (otrace) then close(ufrom) end if stop 300 continue c----------------------------------------------------------------------- c Normal exit, set up values of y for final call to lsoda c Must call lsoda once more to calculate R and S c----------------------------------------------------------------------- dx = max(EPS*abs(ivs(1)),1d-8) dp = max(EPS*abs(ivs(2)),1d-8) y(1) = 0.0d0 y(2) = ivs(1) y(3) = ivs(2) y(4) = ivs(1) y(5) = ivs(2) + dp dp = y(5) - ivs(2) y(6) = ivs(1) + dx dx = y(6) - ivs(1) y(7) = ivs(2) c----------------------------------------------------------------------- c If we are tracing the periodic orbit we will output the c coords. at several times c----------------------------------------------------------------------- if (otrace) then do itout = 1, notpp*iN tout(itout) = 0.0d0 + (itout-1)*T0/notpp end do tout(notpp*iN) = iN*T0 write(ufrom,*) ivs(1), ivs(2) do itout = 2, notpp*iN istate = 1 tbgn = tout(itout-1) tend = tout(itout) call lsoda(fcn,neq,y,tbgn,tend,itol,rtol,atol,itask, & istate,iopt,rwork,lrw,iwork,liw,jac,jt) if (istate .lt. -1) then write(0,*) 'orbitnt: Error return ',istate, & ' from lsoda on last call for N = ',iN close(ufrom) stop else write(ufrom,*) y(2), y(3) end if end do close(ufrom) else c----------------------------------------------------------------------- c If we are not tracing the orbit, we can integrate in one step c----------------------------------------------------------------------- istate = 1 tbgn = 0.0d0 tend = iN*T0 call lsoda(fcn,neq,tbgn,tend,itol,rtol,atol,itask,istate, & iopt,rwork,lrw,iwork,liw,jac,jt) if (istate .lt. -1) then write(0,*) 'orbitnt: Error return ',istate, & ' from lsoda on last call for N = ',iN stop end if end if c----------------------------------------------------------------------- c Compute R and S and write to standard out c----------------------------------------------------------------------- R = (2.0d0 - (y(6)-y(2))/dx - (y(5)-y(3))/dp)/4.0d0 S = y(1) write(*,*) iN, R, S end do stop c----------------------------------------------------------------------- c usage message c----------------------------------------------------------------------- 900 continue write(0,*) 'usage: orbitnt [ ]' write(0,*) write(0,*) 'x0,p0: initial values in orbit search' write(0,*) 'ltol: optional tolerance for lsoda' write(0,*) 'ntol: optional tolerance for Newton convergence' stop end c======================================================================= c dv12norm: Returns 12-norm of double precision vector c (by Matt Choptuik) c======================================================================= real*8 function dv12norm(v,n) implicit none integer n real*8 v(n) integer i dv12norm = 0.0d0 do i = 1, n dv12norm = dv12norm + v(i)*v(i) end do if (n .gt. 0) then dv12norm = sqrt(dv12norm/n) end if return end c======================================================================= c drelavs: Function useful for 'relativizing' quantity being c monitored for convergence (by Matt Choptuik) c======================================================================= real*8 function drelabs(dx,x,xfloor) implicit none real*8 dx, x, xfloor if (abs(x) .lt. abs(xfloor)) then drelabs = abs(dx) else drelabs = abs(dx/x) end if return end