'Program LEVI.BAS (c) 2008 by J. C. Sprott 'Finds initial conditions for a system of ODEs that give bounded chaotic solutions 'using Levi-Civita regularization 'Use PowerBasic Console Compiler '#DEBUG ERROR ON DEFEXT a-z FUNCTION PBMAIN() ne&=4 'Number of equations dr=1e-5## 'Separation of orbits h=.001## 'Iteration step size nm&=1e4 'Number of iterates per test lemin=0.001## 'Minimum LE to qualify as chaotic prec=.5## 'Relative precision in calculated LE bail&=2e9 'Bailout condition (maximum # of iterations) xbound=1e6## 'Upper bound on attractor size ybound=0## 'LOWER BOUND ON y dmax&=2 'Maximum number of digits of parameter resolution eps=1e-9## 'Size of search space nudge=0## 'Amount to nudge parameters back toward zero dr2=dr*dr RANDOMIZE TIMER DIM x(ne&),dxdt(ne&),xo(ne&),xobest(ne&) DIM pi AS GLOBAL EXT pi=3.141592653589794## DIM energy AS GLOBAL EXT CONSOLE NAME EXE.NAME\$ CLS CONSOLE SET FOCUS DO REDIM xe(ne&),dx(ne&),lehist(nm&-1) xo(1)=0.3##: xo(2)=0.5##: xo(3)=0##: xo(4)=0.4## FOR j&=1 TO ne& 'xo(j&)=ROUND(xobest(j&)+eps*MAX(1,ABS(xobest(j&)))*(gauss-nudge*SGN(CINT(xobest(j&)))),dmax&) x(j&)=xo(j&) xe(j&)=x(j&)+dr/SQR(ne&) 'Perturb nearby initial condition NEXT j& m=1.00025##: q=0.125## energy=m*xo(3)*xo(3)+xo(4)*xo(4)/(xo(2)*xo(2))+2##*q/(xo(2)*xo(2))-2##/SQR(m*m*xo(1)*xo(1)+xo(2)^4) n&=0: ltot=0##: tr=0##: ntot&=0: ymin=10 DO INCR n& CALL rk5(x(),dxdt(),ne&,h,x()) CALL rk5(xe(),dxdt(),ne&,h,xe()) rs=0## FOR i&=1 TO ne& IF ABS(x(i&))>xbound THEN PRINT"Orbit exceeded bound: x("+TRIM\$(STR\$(i&))+") =";x(i&): EXIT LOOP dx(i&)=xe(i&)-x(i&) 'IF ABS(dx(i&))>6## THEN dx(i&)=dx(i&)-2##*pi*SGN(dx(i&)) '?FIX FOR MOD X rs=rs+dx(i&)*dx(i&) NEXT i& ymin=MIN(ymin,x(2)) 'IF ymin0.1 THEN PRINT"Orbits separating too much: rs =";rs: EXIT LOOP rs=SQR(rs/dr2) IF rs>0## THEN FOR i&=1 TO ne& xe(i&)=ROUND(x(i&)-dx(i&)/rs,99) NEXT i& IF n&>=nm& THEN ltot=ltot+LOG(rs) INCR ntot& FOR i&=1 TO ne& 'Calculate trace of Jacobian xsave=x(i&) x(i&)=x(i&)-0.5##*dr CALL derivs(x(),dxdt(),ne&) dxdt=dxdt(i&) x(i&)=x(i&)+dr CALL derivs(x(),dxdt(),ne&) tr=tr+(dxdt(i&)-dxdt)/dr x(i&)=xsave NEXT i& le=ltot/(h*ntot&) 'Largest Lyapunov exponent lesum=tr/ntot& 'Sum of Lyapunov exponents lehist(ntot& MOD nm&)=le IF ntot&>=nm& AND (ntot& MOD nm&)=nm&-1 THEN'See if LE has converged CALL lyapext(lehist(),le,dle) IF le+dle/prec0.01## THEN PRINT"Energy not conserved: error =";de;"%": WAITKEY\$: EXIT LOOP END IF END IF ELSE n&=0: PRINT"Orbits collapsed together: rs=";rs: EXIT LOOP END IF IF n&>bail& THEN 'Stop before n& overflows CALL best(ebest&,xo(),le,dle,lesum,ntot&,z\$,ymin) PRINT"Done!" WAITKEY\$ EXIT FUNCTION END IF LOOP LOOP END FUNCTION SUB DERIVS (x(), dxdt(), n&) 'Returns the time derivatives dxdt(i&) of x(i&) m=1.00025##: q=0.125## denom=m*m*x(1)*x(1)+x(2)^4 dxdt(1)=2##*x(2)*x(2)*x(3) dxdt(2)=x(4) dxdt(3)=-2##*m*x(1)*x(2)*x(2)/denom^1.5## dxdt(4)=energy*x(2)-m*x(3)*x(3)*x(2)+2##*m*m*x(1)*x(1)*x(2)/denom^1.5## END SUB SUB RK5 (x(), dxdt(), n&, h, xnew()) 'Calls RK4 (4th order Runge-Kutta) repeatedly with adjustable step size eallow=1e-12 'This is the allowed absolute iteration error DIM xo(n&),xh(n&) 'Temporary required arrays FOR i&=1 TO n&: xo(i&)=x(i&): NEXT i& 'Save initial conditions htry=h CALL RK4(xo(),dxdt(),n&,htry,xnew()) 'Take a full step htry=h/2## CALL RK4(xo(),dxdt(),n&,htry,xh()) 'Take 2 half-steps CALL RK4(xh(),dxdt(),n&,htry,xh()) e=0 FOR i&=1 TO n& 'Calculate maximum error in the difference e=MAX(e,ABS(xnew(i&)-xh(i&))) NEXT i& steps&=CEIL(1.1+(e/eallow)^.25) 'Number of steps required to achieve eallow steps&=MIN(steps&,1e4) 'Bailout after 1e4 steps IF steps&>2 THEN 'The two half steps were not enough htry=h/steps& 'so try more CALL RK4(xo(),dxdt(),n&,htry,xh()) FOR i&=2 TO steps& CALL RK4(xh(),dxdt(),n&,htry,xh()) NEXT i& END IF fcor=1##/(steps&^4-1##) 'Fifth order correction in h FOR i&=1 TO n& xnew(i&)=xh(i&)+fcor*(xnew(i&)-xh(i&)) IF ABS(xnew(i&))<1e-999 THEN xnew(i&)=xnew(i&)+1e-999*(RND-.5) 'Perturb off the origin NEXT i& 'IF xnew(1)>pi THEN xnew(1)=xnew(1)-2##*pi '?X IS PERIODIC 'IF xnew(1)<-pi THEN xnew(1)=xnew(1)+2##*pi 'IF xnew(n&)<0 THEN xnew(n&)=xnew(n&)+2##*pi '?Z IS PERIODIC 'IF xnew(n&)>2##*pi THEN xnew(n&)=xnew(n&)-2##*pi ERASE xh,xo END SUB SUB SI4 (X(), DXDT(), N&, H, XNEW()) 'Fourth-order symplectic integrator (Forest90) 'Assumes derivatives are ordered x, vx, y, vy, ... c1 = 0.675603595979828817##*h c2 = -0.175603595979828817##*h d1 = 1.351207191959657630##*h d2 = -1.702414383919315270##*h DIM xt(n&) FOR i&=1 TO n& xt(i&)=x(i&) NEXT i& CALL derivs(xt(), dxdt(), n&) FOR i&=2 TO n& STEP 2 xt(i&)=xt(i&)+c1*dxdt(i&) NEXT i& CALL derivs(xt(), dxdt(), n&) FOR i&=1 TO n& STEP 2 xt(i&)=xt(i&)+d1*dxdt(i&) NEXT i& CALL derivs(xt(), dxdt(), n&) FOR i&=2 TO n& STEP 2 xt(i&)=xt(i&)+c2*dxdt(i&) NEXT i& CALL derivs(xt(), dxdt(), n&) FOR i&=1 TO n& STEP 2 xt(i&)=xt(i&)+d2*dxdt(i&) NEXT i& CALL derivs(xt(), dxdt(), n&) FOR i&=2 TO n& STEP 2 xt(i&)=xt(i&)+c2*dxdt(i&) NEXT i& CALL derivs(xt(), dxdt(), n&) FOR i&=1 TO n& STEP 2 xt(i&)=xt(i&)+d1*dxdt(i&) NEXT i& CALL derivs(xt(), dxdt(), n&) FOR i&=1 TO n& STEP 2 xnew(i&)=xt(i&) xnew(i&+1)=xt(i&+1)+c1*dxdt(i&+1) NEXT i& ERASE xt END SUB SUB RK4 (X(), DXDT(), N&, H, XNEW()) 'Fourth-order Runge-Kutta integrator DIM XT(N&), DXT(N&), DXM(N&) HH = H * .5## H6 = H / 6## CALL DERIVS(X(), DXDT(), N&) FOR I& = 1 TO N& XT(I&) = X(I&) + HH * DXDT(I&) NEXT I& CALL DERIVS(XT(), DXT(), N&) FOR I& = 1 TO N& XT(I&) = X(I&) + HH * DXT(I&) NEXT I& CALL DERIVS(XT(), DXM(), N&) FOR I& = 1 TO N& XT(I&) = X(I&) + H * DXM(I&) DXM(I&) = DXT(I&) + DXM(I&) NEXT I& CALL DERIVS(XT(), DXT(), N&) FOR I& = 1 TO N& XNEW(I&) = X(I&) + H6 * (DXDT(I&) + DXT(I&) + 2## * DXM(I&)) NEXT I& ERASE DXM, DXT, XT END SUB SUB lyapext(lehist(), le, dle) 'Estimates the asymptotic Lyapunov exponent (le) 'and its error (dle) from the recent history (lehist) n&=UBOUND(lehist): n3&=n&/3 'See how many points were stored l1=0##: l2=0##: l3=0## FOR i&=0 TO n3& 'Average the LE over each third l1=l1+lehist(i&) NEXT i& FOR i&=n3& TO 2*n3& l2=l2+lehist(i&) NEXT i& FOR i&=2*n3& TO n& l3=l3+lehist(i&) NEXT i& l1=l1/(n3&+1): l2=l2/(n3&+1): l3=l3/(n3&+1) IF ABS(l3-l2)= 1## OR R = 0## GAUSS = V2 * SQR(-2## * LOG(R) / R) END FUNCTION FUNCTION tanh (x) 'Returns hyperbolic tangent of x IF ABS(x)>22 THEN tanh = SGN(x) ELSE tanh = 1## - 2## / (EXP(2## * x) + 1##) END IF END FUNCTION FUNCTION sinh (x) 'Returns hyperbolic sine of x sinh = (EXP(x) - EXP(-x)) / 2## END FUNCTION FUNCTION cosh (x) 'Returns hyperbolic cosine of x cosh = (EXP(x) + EXP(-x)) / 2## END FUNCTION