'Program INITCOND.BAS (c) 2008 by J. C. Sprott 'Finds initial conditions for a system of ODEs that give bounded chaotic solutions 'Use PowerBasic Console Compiler #DEBUG ERROR OFF DEFEXT a-z FUNCTION PBMAIN() ne&=4 'Number of equations dr=1e-8## '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.01 'LOWER BOUND ON y dmax&=1 'Maximum number of digits of parameter resolution eps=1## 'Size of search space nudge=0.5## '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## CONSOLE NAME EXE.NAME\$ CLS CONSOLE SET FOCUS DO REDIM xe(ne&),dx(ne&),lehist(nm&-1) 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& 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(3)) IF ymindr 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/prec2 THEN 'The two half steps were not enough htry=h/steps& 'so try more CALL SI4(xo(),dxdt(),n&,htry,xh()) FOR i&=2 TO steps& CALL SI4(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&)) '?DON'T USE THIS WITH SI4 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