'Calculate theoretical values of S(j) DEFEXT a-z DIM file\$(15),x(3),dxdt(3) nmax&=512 'Number of data points to save to file (plus 18) CLS PRINT"A: Henon map": file\$(1)="henon.dat" PRINT"B: Chaostsa preface": file\$(2)="preface.dat" PRINT"C: Time-delayed Henon map": file\$(3)="delayhen.dat" PRINT"D: Logistic map": file\$(4)="logmap.dat" PRINT"E: Goutte map": file\$(5)="goutte.dat" PRINT"F: Cantor set": file\$(6)="cantor.dat" PRINT"G: Lorenz 3-d map": file\$(7)="lor3dmap.dat" PRINT"H: Jerk system": file\$(8)="quadjerk.dat" PRINT"I: Lorenz attractor": file\$(9)="lorenz.dat" PRINT"J: Rossler attractor": file\$(10)="rosslery.dat" PRINT"K: Henon map x+y": file\$(11)="henonxy.dat" PRINT"L: Logistic map x+y": file\$(12)="logmapxy.dat" PRINT"M: Quadratic map x+y": file\$(13)="quadxy.dat" PRINT"N: Binary shift map": file\$(14)="shiftmap.dat" PRINT"O: Hyperchaotic map": file\$(15)="hypchaos.dat" PRINT PRINT"Case (A-O)?"; WHILE q\$="": q\$=UCASE\$(INKEY\$): WEND SELECT CASE q\$ 'Set initial conditions CASE "A": x1=0.1 CASE "B": x1=0.1 CASE "C": x1=0.1 CASE "D": x1=0.1 CASE "E": x1=0.1 CASE "F": x1=0.1 CASE "G": x1=0.1: x2=0.5: x3=-1 CASE "H": h=0.1##: x1=0: x2=-0.1##*h: x3=-0.2##*h CASE "I": x(1)=0.1 CASE "J": x(1)=-9 CASE "K": x1=0.1 CASE "L": x1=0.1 CASE "M": x1=0.2 CASE "N": x1=0.1 CASE "O": x1=1: x2=0.1: x3=0 END SELECT CLS n&=0 WHILE INKEY\$<>CHR\$(27) IF i&=7000 THEN OPEN file\$(ASC(q\$)-64) FOR OUTPUT AS #1 INCR i& SELECT CASE q\$ CASE"A" 'Henon map x=1##-1.4##*x1*x1+0.3##*x2 x2=x1 x1=x s0=s0+1## s1=s1+ABS(2.8##*x) s2=s2+ABS(0.3##) IF (i& MOD 1e5)=0 THEN PRINT USING"####.#####";x, s0/i&, s1/i&, s2/i& CASE"B" 'Chaostsa preface x=x1*x1-0.2##*x1-0.9##*x2+0.6##*x3 x3=x2 x2=x1 x1=x s0=s0+0 s1=s1+ABS(2##*x-0.2##) s2=s2+ABS(-0.9##) s3=s3+ABS(0.6##) IF (i& MOD 1e5)=0 THEN PRINT USING"####.#####";x, s0/i&, s1/i&, s2/i&, s3/i& CASE"C" 'Time-delayed Henon map x=1##-1.6##*x1*x1+0.1##*x4 x4=x3 x3=x2 x2=x1 x1=x s0=s0+1 s1=s1+ABS(-3.2##*x) s4=s4+ABS(0.1##) IF (i& MOD 1e5)=0 THEN PRINT USING"####.#####";x, s0/i&, s1/i&, s2/i&, s3/i&, s4/i& CASE "D" 'Logistic map x=4##*x1*(1##-x1) x1=x s1=s1+4##*ABS(1##-2##*x) IF (i& MOD 1e5)=0 THEN PRINT USING"####.#####";x, s0/i&, s1/i& CASE"E" 'Goutte map x=1##-1.4##*x2*x2+0.3##*x4 x4=x3 x3=x2 x2=x1 x1=x s0=s0+1 s2=s2+ABS(-2.8##*x) s4=s4+ABS(0.3##) IF (i& MOD 1e5)=0 THEN PRINT USING"####.#####";x, s0/i&, s1/i&, s2/i&, s3/i&, s4/i& CASE "F" 'Cantor set x=3.5699456718##*x1*(1##-x1) x1=x s1=s1+3.5699456718##*ABS(1##-2##*x) IF (i& MOD 1e5)=0 THEN PRINT USING"####.#####";x, s0/i&, s1/i& CASE"G" 'Lorenz 3-D map x=x1*x2-x3 x3=x2 x2=x1 x1=x s0=s0+0 s1=s1+ABS(x) s2=s2+ABS(x) s3=s3+ABS(-1) IF (i& MOD 1e5)=0 THEN PRINT USING"####.#####";x, s0/i&, s1/i&, s2/i&, s3/i& CASE"H" 'Jerk system (Case M03 discretized) x=3##*x1-3##*x2+x3 x=x+h*(x1-x2/2##) x=x-h*h*(x1/2##-x2+x3/6##) x=x+h*h*h*x1*(x1-1##) d=(1##+h/2##+h*h/3##) x=x/d x3=x2 x2=x1 x1=x s0=s0+0 s1=s1+ABS(3+h-h*h/2+2*h*h*h*x-h*h*h)/d s2=s2+ABS(-3-h/2+h*h)/d s3=s3+ABS(1-h*h/6)/d IF (i& MOD 1e5)=0 THEN PRINT USING"####.#####";x, s0/i&, s1/i&, s2/i&, s3/i& CASE"I" 'Lorenz attractor h=0.1## CALL RK5 (x(), dxdt(), 3, h, x()) x=x(1) s1=s1+ABS(dxdt(1)) IF (i& MOD 1e5)=0 THEN PRINT USING"####.#####";x, s0/i&, s1/i& CASE"J" 'Rossler attractor h=0.2## CALL RK5 (x(), dxdt(), 3, h, x()) x=x(2) s1=s1+ABS(dxdt(2)) IF (i& MOD 1e5)=0 THEN PRINT USING"####.#####";x, s0/i&, s1/i& CASE"K" 'Henon map x+y x=1##-1.4##*x1*x1+0.3##*x2 x2=x1 x1=x x=x1+x2 s0=0 s1=0 s2=0 IF (i& MOD 1e5)=0 THEN PRINT USING"####.#####";x, s0/i&, s1/i&, s2/i& CASE"L" 'Logistic map x+y x=4##*x1*(1##-x1) x2=x1 x1=x x=x1+x2 s0=0 s1=0 s2=0 IF (i& MOD 1e5)=0 THEN PRINT USING"####.#####";x, s0/i&, s1/i&, s2/i& CASE"M" 'Quadratic map x+y x=1##-1.5##*x1*x1 x2=x1 x1=x x=x1+x2 s0=0 s1=0 s2=0 IF (i& MOD 1e5)=0 THEN PRINT USING"####.#####";x, s0/i&, s1/i&, s2/i& CASE "N" 'Binary shift map x=1.99999999999999999##*x1 IF x>1## THEN x=x-1## x1=x s1=s1+1.999999999999999## IF (i& MOD 1e5)=0 THEN PRINT USING"####.#####";x, s0/i&, s1/i& CASE "O" 'Hyperchaotic map x=1.76##-x2*x2-0.1##*x3 x3=x2 x2=x1 x1=x s0=s0+1.76## s1=s1+0## s2=s2+ABS(2##*x) s3=s3+0.1## IF (i& MOD 1e5)=0 THEN PRINT USING"####.#####";x, s0/i&, s1/i&, s2/i&, s3/i& CASE ELSE PRINT"Not a legal case" END END SELECT IF i&>7000 AND i&<=7000+nmax&+18 THEN PRINT#1,x IF i&=7000+nmax&+18 THEN CLOSE #1 PRINT"File "+file\$(ASC(q\$)-64)+" written" END IF WEND END SUB DERIVS (x(), dxdt(), n&) 'Returns the time derivatives dxdt(i&) of x(i&) 'Chaotic flow SHARED q\$ SELECT CASE q\$ CASE"I" 'Lorenz dxdt(1)=10##*(x(2)-x(1)) dxdt(2)=-x(1)*x(3)+28##*x(1)-x(2) dxdt(3)=x(1)*x(2)-8##*x(3)/3## CASE"J" 'Rossler dxdt(1)=-x(2)-x(3) dxdt(2)=x(1)+0.2##*x(2) dxdt(3)=0.2##+x(3)*(x(1)-5.7##) END SELECT 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&,1000) 'Bailout after 100 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 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