'Program REGIONS.BAS (c) 2008 by J. C. Sprott 'Makes a scatter plot of different dynamic behaviors in 2-D parameter space 'Use PowerBasic Console Compiler '#DEBUG ERROR ON '#DEBUG DISPLAY ON DEFEXT a-z GLOBAL dots&,maxdots& GLOBAL fhndl& FUNCTION PBMAIN() ne&=3 'Number of equations np&=2 'Number of parameters dr=1e-8## 'Separation of orbits h=.1## 'Iteration step size nm&=2e3 'Number of iterates per test xbound=1e3## 'Upper bound on attractor size w&=512 'Width of plot h&=512 'Height of plot b&=40 'Border size dmax&=4 'Maximum number of digits of parameter resolution DIM x(ne&),dxdt(ne&),xo(ne&),xe(ne&),dx(ne&),amin(np&),amax(np&),lehist(nm&-1) DIM wr(ne&),wi(ne&),jacob(ne&,ne&) DIM xobest(ne&),abest(np&) DIM xeq(ne&),lehist(nm&-1),x1hist(nm&-1),d&(np&) DIM a(np&) AS GLOBAL EXT ax\$="a": amin(1)=0##: amax(1)=1## ay\$="b": amin(2)=0##: amax(2)=3## xticks&=ticks&(amin(1),amax(1),14) 'Number of ticks on x-axis yticks&=ticks&(amin(2),amax(2),14) 'Number of ticks on y-axis dr2=dr*dr FONT NEW "MS Sans Serif",14,1 TO fhndl& RANDOMIZE TIMER DO CONSOLE SET FOCUS COLOR 7 CLS PRINT,," REGIONS" PRINT,,"(c) 2012 by J. C. Sprott" PRINT PRINT,"Press D for dynamic regions" PRINT," R for hyperchaotic regions" PRINT," S for stability of equilibria" PRINT," H for homoclinic orbit" PRINT," L for maximum Lyapunov exponent" PRINT," K for maximum Kaplan-Yorke dimension" PRINT," A for maximum asymmetry" PRINT," to end program: "; q\$=UCASE\$(WAITKEY\$) SELECT CASE q\$ CASE \$ESC EXIT FUNCTION CASE"D" lemin=0.001## 'Minimum LE to qualify as chaotic prec=1## 'Relative precision in calculated LE bail&=1e5 'Bailout condition (maximum # of iterations) DESKTOP GET SIZE TO ddw&, dh& CONSOLE SET LOC w&+b&+6,0 CLS GRAPHIC WINDOW EXE.FULL\$, 0, 0, w&+b&, h&+30 TO hWin??? GRAPHIC ATTACH hWin???, 0, REDRAW GRAPHIC COLOR %BLACK, %WHITE GRAPHIC CLEAR GRAPHIC SET FONT fhndl& xm\$=TRIM\$(STR\$(amax(2))) IF LEFT\$(xm\$,1)="." THEN xm\$="0"+xm\$ REPLACE"-." WITH"-0." IN xm\$ GRAPHIC TEXT SIZE xm\$ TO xmw, xmh GRAPHIC SET POS(b&-5-xmw,0) GRAPHIC PRINT xm\$ xm\$=ay\$ 'a(2) - Vertical axis GRAPHIC TEXT SIZE xm\$ TO xmw, xmh GRAPHIC SET POS(b&-5-xmw,h&/2-13) GRAPHIC PRINT xm\$ xm\$=TRIM\$(STR\$(amin(2))) REPLACE"-." WITH"-0." IN xm\$ IF LEFT\$(xm\$,1)="." THEN xm\$="0"+xm\$ GRAPHIC TEXT SIZE xm\$ TO xmw, xmh GRAPHIC SET POS(b&-5-xmw,h&-xmh) GRAPHIC PRINT xm\$ xm\$=TRIM\$(STR\$(amin(1))) IF LEFT\$(xm\$,1)="." THEN xm\$="0"+xm\$ REPLACE"-." WITH"-0." IN xm\$ GRAPHIC TEXT SIZE xm\$ TO xmw, xmh GRAPHIC SET POS(b&,h&) GRAPHIC PRINT xm\$ xm\$=ax\$ 'a(1) - Horizontal axis GRAPHIC TEXT SIZE xm\$ TO xmw, xmh GRAPHIC SET POS(b&+w&/2-xmw/2,h&) GRAPHIC PRINT xm\$ xm\$=TRIM\$(STR\$(amax(1))) IF LEFT\$(xm\$,1)="." THEN xm\$="0"+xm\$ REPLACE"-." WITH"-0." IN xm\$ GRAPHIC TEXT SIZE xm\$ TO xmw, xmh GRAPHIC SET POS(b&+w&-xmw,h&) GRAPHIC PRINT xm\$ GRAPHIC BOX (b&, 0) - (b&+w&, h&+1) 'GRAPHIC RENDER exe.name\$+".bmp",(0,0)-(w&+b&-1,h&+30-1) GRAPHIC REDRAW CONSOLE SET FOCUS ddw&=1024: done&=0 DO IF ddw&>0 THEN ddw&=ddw&/2 ELSE done&=1 'Box size (gradually improve the resolution) FOR xp&=0 TO w&-ddw& STEP ddw& FOR yp&=0 TO h&-ddw& STEP ddw& IF done&=1 THEN GRAPHIC GET PIXEL(b&+xp&,yp&) TO c& IF c&<>%RGB_WHITE THEN ITERATE 'This point has already been colored GRAPHIC GET PIXEL(b&+xp&,yp&-1) TO cn& 'North point GRAPHIC GET PIXEL(b&+xp&,yp&+1) TO cs& 'South point GRAPHIC GET PIXEL(b&+xp&+1,yp&) TO ce& 'East point GRAPHIC GET PIXEL(b&+xp&-1,yp&) TO cw& 'West point IF cn&=c& AND cs&=c& AND ce&=c& AND cw&=c& THEN ITERATE 'This point is surrounded by white ELSE IF ddw&xbound OR rs>dr THEN 'Unbounded orbit dyn\$="U" EXIT LOOP 'EXIT,ITERATE END IF rs=SQR(rs/dr2) IF rs>0## THEN FOR i&=1 TO ne& xe(i&)=x(i&)-dx(i&)/rs NEXT i& END IF IF n&>=nm& THEN ltot=ltot+LOG(rs) INCR ntot& xtot=xtot+x(1): xabs=xabs+ABS(x(1)): asym=xtot/xabs 'Calculate asymmetry IF ne&>1 THEN ytot=ytot+x(2): yabs=yabs+ABS(x(2)): asym=asym+ytot/yabs 'IF ne&>2 THEN ztot=ztot+x(3): zabs=zabs+ABS(x(3)): asym=asym+ztot/zabs le=ltot/(h*ntot&) 'Largest Lyapunov exponent lehist(ntot& MOD nm&)=le IF ntot&>=nm& AND (ntot& MOD nm&)=nm&-1 THEN'See if LE has converged CALL mouseloc(w&,h&,b&,amin(),amax(),ax\$,ay\$) CALL lyapext(lehist(),le,dle) IF le>lemin+dle/prec THEN dyn\$="C": EXIT LOOP IF ABS(le)bail& THEN 'It's taking too long SELECT CASE le CASE >lemin: dyn\$="C" CASE -lemin TO lemin: dyn\$="P" CASE <-lemin: dyn\$="S" END SELECT EXIT LOOP END IF END IF END IF LOOP 'IF dyn\$="C" AND ABS(asym)<0.5## THEN dyn\$="D" 'Symmetric strange attractor 'IF dyn\$="P" AND ABS(asym)<0.5## THEN dyn\$="Q" 'Symmetric limit cycle SELECT CASE dyn\$ 'Select colors for different dynamics CASE "U": COLOR 7: GRAPHIC COLOR %RGB_WHITE CASE "C": COLOR 15: GRAPHIC COLOR %RGB_BLACK CASE "D": COLOR 13: GRAPHIC COLOR %RGB_MAGENTA CASE "P": COLOR 11: GRAPHIC COLOR %RGB_CYAN CASE "Q": COLOR 7: GRAPHIC COLOR %RGB_GOLD CASE "S": COLOR 12: GRAPHIC COLOR %RGB_RED CASE ELSE: PRINT"Should never get here!": BEEP: WAITKEY\$ END SELECT PRINT FORMAT\$(a(1),"* ###.0000"); PRINT FORMAT\$(a(2),"* ###.0000"); PRINT FORMAT\$(le," 0.0000; -0.0000");CHR\$(177);FORMAT\$(dle,"0.0000"); PRINT" ";dyn\$;n&+2 IF ddw&=1 THEN 'Plot the corresponding color GRAPHIC SET PIXEL(b&+xp&,yp&) ELSE GRAPHIC BOX(b&+xp&,yp&)-(b&+xp&+ddw&,yp&+ddw&),0,,-1 END IF GRAPHIC REDRAW IF done&>0 AND dyn\$<>"U" THEN INCR done& NEXT yp& GRAPHIC COLOR %BLACK GRAPHIC BOX (wo&+b&, ho&) - (wo&+b&+w&, ho&+h&+1) FOR i&=1 TO xticks&-1 j&=b&+i&*w&/xticks& GRAPHIC LINE(j&,0)-(j&,h&/36) GRAPHIC LINE(j&,h&)-(j&,h&-h&/36) NEXT i& FOR i&=1 TO yticks&-1 j&=i&*h&/yticks& GRAPHIC LINE(b&,j&)-(b&+w&/36,j&) GRAPHIC LINE(b&+w&,j&)-(b&+w&-w&/36,j&) NEXT i& GRAPHIC REDRAW IF INKEY\$=\$ESC THEN EXIT LOOP NEXT xp& GRAPHIC REDRAW GRAPHIC SAVE EXE.NAME\$+".bmp" IF done&=1 THEN EXIT LOOP LOOP x1&=b& 'Label the different regions in a contrasting color lastc&=%YELLOW 'A color not used FOR xp&=b& TO w&+b& GRAPHIC GET PIXEL(xp&,h&/2) TO c& IF c&<>cold& THEN x2&=xp& xpo&=x1&+(x2&-x1&)/2-6 IF x2&-x1&>25 AND cold&<>lastc& THEN GRAPHIC SET POS(xpo&,h&/2-12) SELECT CASE cold& CASE %RGB_WHITE: dyn\$="U": GRAPHIC COLOR %BLACK,-2 CASE %RGB_BLACK: dyn\$="C": GRAPHIC COLOR %WHITE,-2 CASE %RGB_MAGENTA: dyn\$="D": GRAPHIC COLOR %WHITE,-2 CASE %RGB_CYAN: dyn\$="P": GRAPHIC COLOR %BLACK,-2 CASE %RGB_GOLD: dyn\$="Q": GRAPHIC COLOR %BLACK,-2 CASE %RGB_RED: dyn\$="S": GRAPHIC COLOR %WHITE,-2 END SELECT GRAPHIC PRINT dyn\$; lastc&=cold& END IF x1&=xp& cold&=c& END IF NEXT xp& GRAPHIC SAVE EXE.NAME\$+".bmp" WAITKEY\$ GRAPHIC WINDOW END CASE"R" lemin=0.002## 'Minimum LE to qualify as hyperchaotic prec=1## 'Relative precision in calculated LE bail&=1e5 'Bailout condition (maximum # of iterations) nn&=ne&*(ne&+1) 'total number of odes (nonlinear + linear) irate&=10 'integration steps per reorthonormalization rd&=8 'resolution of dots (1 for maximum resolution) DESKTOP GET SIZE TO ddw&, dh& CONSOLE SET LOC w&+b&+6,0 CLS GRAPHIC WINDOW EXE.FULL\$, 0, 0, w&+b&, h&+30 TO hWin??? GRAPHIC ATTACH hWin???, 0, REDRAW GRAPHIC COLOR %BLACK, %WHITE GRAPHIC CLEAR GRAPHIC SET FONT fhndl& xm\$=TRIM\$(STR\$(amax(2))) IF LEFT\$(xm\$,1)="." THEN xm\$="0"+xm\$ REPLACE"-." WITH"-0." IN xm\$ GRAPHIC TEXT SIZE xm\$ TO xmw, xmh GRAPHIC SET POS(b&-5-xmw,0) GRAPHIC PRINT xm\$ xm\$=ay\$ 'a(2) - Vertical axis GRAPHIC TEXT SIZE xm\$ TO xmw, xmh GRAPHIC SET POS(b&-5-xmw,h&/2-13) GRAPHIC PRINT xm\$ xm\$=TRIM\$(STR\$(amin(2))) REPLACE"-." WITH"-0." IN xm\$ IF LEFT\$(xm\$,1)="." THEN xm\$="0"+xm\$ GRAPHIC TEXT SIZE xm\$ TO xmw, xmh GRAPHIC SET POS(b&-5-xmw,h&-xmh) GRAPHIC PRINT xm\$ xm\$=TRIM\$(STR\$(amin(1))) IF LEFT\$(xm\$,1)="." THEN xm\$="0"+xm\$ REPLACE"-." WITH"-0." IN xm\$ GRAPHIC TEXT SIZE xm\$ TO xmw, xmh GRAPHIC SET POS(b&,h&) GRAPHIC PRINT xm\$ xm\$=ax\$ 'a(1) - Horizontal axis GRAPHIC TEXT SIZE xm\$ TO xmw, xmh GRAPHIC SET POS(b&+w&/2-xmw/2,h&) GRAPHIC PRINT xm\$ xm\$=TRIM\$(STR\$(amax(1))) IF LEFT\$(xm\$,1)="." THEN xm\$="0"+xm\$ REPLACE"-." WITH"-0." IN xm\$ GRAPHIC TEXT SIZE xm\$ TO xmw, xmh GRAPHIC SET POS(b&+w&-xmw,h&) GRAPHIC PRINT xm\$ GRAPHIC BOX (b&, 0) - (b&+w&, h&+1) 'GRAPHIC RENDER EXE.NAME\$+".bmp",(0,0)-(w&+b&-1,h&+30-1) GRAPHIC REDRAW CONSOLE SET FOCUS ddw&=1024: done&=0 DO IF ddw&>1 THEN ddw&=ddw&/2 ELSE done&=1 'Box size (gradually improve the resolution) FOR xp&=0 TO w&-ddw& STEP ddw& FOR yp&=0 TO h&-ddw& STEP ddw& IF done&>0 THEN GRAPHIC GET PIXEL(b&+xp&,yp&) TO c& IF c&<>%RGB_WHITE THEN ITERATE 'This point has already been colored GRAPHIC GET PIXEL(b&+xp&,yp&-1) TO cn& 'North point GRAPHIC GET PIXEL(b&+xp&,yp&+1) TO cs& 'South point GRAPHIC GET PIXEL(b&+xp&+1,yp&) TO ce& 'East point GRAPHIC GET PIXEL(b&+xp&-1,yp&) TO cw& 'West point IF cn&=c& AND cs&=c& AND ce&=c& AND cw&=c& THEN ITERATE 'This point is surrounded by white ELSE IF ddw&xbound THEN 'Unbounded orbit dyn\$="U" 'EXIT for END IF 'construct new orthonormal basis by gram-schmidt znorm(1)=0## 'normalize the first vector FOR j&=1 TO ne& znorm(1)=znorm(1)+x(ne&*j&+1)^2 NEXT j& znorm(1)=SQR(znorm(1)) FOR j&=1 TO ne& x(ne&*j&+1)=x(ne&*j&+1)/znorm(1) NEXT j& 'generate a new orthonormal set FOR j&=2 TO ne& FOR k&=1 TO j&-1 'make j-1 gsr coefficients gsc(k&)=0## FOR l&=1 TO ne& gsc(k&)=gsc(k&)+x(ne&*l&+j&)*x(ne&*l&+k&) NEXT l& NEXT k& FOR k&=1 TO ne& 'construct a new vector FOR l&=1 TO j&-1 x(ne&*k&+j&)=x(ne&*k&+j&)-gsc(l&)*x(ne&*k&+l&) NEXT l& NEXT k& znorm(j&)=0## 'calculate the vector's norm FOR k&=1 TO ne& znorm(j&)=znorm(j&)+x(ne&*k&+j&)^2 NEXT k& znorm(j&)=SQR(znorm(j&)) FOR k&=1 TO ne& 'normalize the new vector x(ne&*k&+j&)=x(ne&*k&+j&)/znorm(j&) NEXT k& NEXT j& FOR k&=1 TO ne& 'update the running vector magnitudes IF znorm(k&)>0 THEN ltot(k&)=ltot(k&)+LOG(znorm(k&)) NEXT k& INCR m& IF (m& MOD nm&)=0 THEN 'normalize exponent and print every nm& iterations 'PRINT m&;"LEs ="; lsum=0##: kmax&=0 FOR k&=1 TO ne& le(k&)=ltot(k&)/t lsum=lsum+le(k&) IF lsum>0 THEN lsum0=lsum: kmax&=k& 'PRINT USING\$("##.#### ", le(k&)); NEXT k& IF ltot(1)>0## AND kmax&10*nm& OR 2*m&>bail&) AND le(1)<-lemin THEN dyn\$="S": EXIT LOOP minle=1## FOR k&=1 TO ne& IF ABS(le(k&))=bail& THEN 'Give up and take best guess SELECT CASE kmin& CASE 1: dyn\$="P" CASE 2: dyn\$="C" CASE >2: dyn\$="H" END SELECT EXIT LOOP END IF IF minle>lemin THEN ITERATE 'Keep going until one LE is zero IF le(2)>lemin THEN dyn\$="H": EXIT LOOP IF le(2)lemin THEN dyn\$="C": EXIT LOOP IF ABS(le(1))0 AND dyn\$<>"U" THEN INCR done& NEXT yp& GRAPHIC COLOR %BLACK GRAPHIC BOX (wo&+b&, ho&) - (wo&+b&+w&, ho&+h&+1) FOR i&=1 TO xticks&-1 j&=b&+i&*w&/xticks& GRAPHIC LINE(j&,0)-(j&,h&/36) GRAPHIC LINE(j&,h&)-(j&,h&-h&/36) NEXT i& FOR i&=1 TO yticks&-1 j&=i&*h&/yticks& GRAPHIC LINE(b&,j&)-(b&+w&/36,j&) GRAPHIC LINE(b&+w&,j&)-(b&+w&-w&/36,j&) NEXT i& GRAPHIC REDRAW IF INKEY\$=\$ESC THEN EXIT LOOP NEXT xp& GRAPHIC SAVE EXE.NAME\$+".bmp" IF done&=1 THEN EXIT LOOP LOOP x1&=b& 'Label the different regions in a contrasting color lastc&=%YELLOW 'A color not used FOR xp&=b& TO w&+b& GRAPHIC GET PIXEL(xp&,h&/2) TO c& IF c&<>cold& THEN x2&=xp& xpo&=x1&+(x2&-x1&)/2-6 IF x2&-x1&>25 AND cold&<>lastc& THEN GRAPHIC SET POS(xpo&,h&/2-12) SELECT CASE cold& CASE %RGB_RED: dyn\$="H": GRAPHIC COLOR %WHITE,-2 CASE %RGB_WHITE: dyn\$="U": GRAPHIC COLOR %BLACK,-2 CASE %RGB_BLACK: dyn\$="C": GRAPHIC COLOR %WHITE,-2 CASE %RGB_CYAN: dyn\$="P": GRAPHIC COLOR %BLACK,-2 CASE %RGB_GOLD: dyn\$="Q": GRAPHIC COLOR %BLACK,-2 CASE %RGB_GREEN: dyn\$="S": GRAPHIC COLOR %WHITE,-2 END SELECT GRAPHIC PRINT dyn\$; lastc&=cold& END IF x1&=xp& cold&=c& END IF NEXT xp& GRAPHIC COLOR %BLACK FOR i&=1 TO xticks&-1 j&=b&+i&*w&/xticks& GRAPHIC LINE(j&,0)-(j&,h&/36) GRAPHIC LINE(j&,h&)-(j&,h&-h&/36) NEXT i& FOR i&=1 TO yticks&-1 j&=i&*h&/yticks& GRAPHIC LINE(b&,j&)-(b&+w&/36,j&) GRAPHIC LINE(b&+w&,j&)-(b&+w&-w&/36,j&) NEXT i& GRAPHIC REDRAW GRAPHIC SAVE EXE.NAME\$+".bmp" WAITKEY\$ GRAPHIC WINDOW END CASE"S" DESKTOP GET SIZE TO ddw&, dh& CONSOLE SET LOC w&+b&+6,0 CLS GRAPHIC WINDOW EXE.FULL\$, 0, 0, w&+b&, h&+30 TO hWin??? GRAPHIC ATTACH hWin???, 0, REDRAW GRAPHIC COLOR %BLACK, %WHITE GRAPHIC CLEAR GRAPHIC SET FONT fhndl& xm\$=TRIM\$(STR\$(amax(2))) IF LEFT\$(xm\$,1)="." THEN xm\$="0"+xm\$ REPLACE"-." WITH"-0." IN xm\$ GRAPHIC TEXT SIZE xm\$ TO xmw, xmh GRAPHIC SET POS(b&-5-xmw,0) GRAPHIC PRINT xm\$ xm\$=ay\$ 'a(2) - Vertical axis GRAPHIC TEXT SIZE xm\$ TO xmw, xmh GRAPHIC SET POS(b&-5-xmw,h&/2-13) GRAPHIC PRINT xm\$ xm\$=TRIM\$(STR\$(amin(2))) REPLACE"-." WITH"-0." IN xm\$ IF LEFT\$(xm\$,1)="." THEN xm\$="0"+xm\$ GRAPHIC TEXT SIZE xm\$ TO xmw, xmh GRAPHIC SET POS(b&-5-xmw,h&-xmh) GRAPHIC PRINT xm\$ xm\$=TRIM\$(STR\$(amin(1))) IF LEFT\$(xm\$,1)="." THEN xm\$="0"+xm\$ REPLACE"-." WITH"-0." IN xm\$ GRAPHIC TEXT SIZE xm\$ TO xmw, xmh GRAPHIC SET POS(b&,h&) GRAPHIC PRINT xm\$ xm\$=ax\$ 'a(1) - Horizontal axis GRAPHIC TEXT SIZE xm\$ TO xmw, xmh GRAPHIC SET POS(b&+w&/2-xmw/2,h&) GRAPHIC PRINT xm\$ xm\$=TRIM\$(STR\$(amax(1))) IF LEFT\$(xm\$,1)="." THEN xm\$="0"+xm\$ REPLACE"-." WITH"-0." IN xm\$ GRAPHIC TEXT SIZE xm\$ TO xmw, xmh GRAPHIC SET POS(b&+w&-xmw,h&) GRAPHIC PRINT xm\$ GRAPHIC BOX (b&, 0) - (b&+w&, h&+1) GRAPHIC REDRAW CONSOLE SET FOCUS FOR yp&=1 TO h&-1: FOR xp&=b&+1 TO b&+w&-2 a(2)=amin(2)+(amax(2)-amin(2))*(h&-yp&)/h& a(1)=amin(1)+(amax(1)-amin(1))*(xp&-b&)/w& xeq(1)=0##: xeq(2)=0##: xeq(3)=0## 'Equilibrium point 'xeq(3)=-1##: xeq(2)=sqr(-a(2)): xeq(1)=xeq(2) CALL jacobian(xeq(),dxdt(),ne&,dr,jacob(),h) CALL balanc(jacob(), ne&) 'Balance Jacobian matrix for improved accuracy CALL elmhes(jacob(), ne&) 'Reduce Jacobian to Hessenberg form CALL hqr(jacob(), ne&, wr(), wi()) 'Calculate eigenvalues wrmin=1##: wrmax=-wrmin: wimax=0## FOR i&=1 TO ne& wrmax=MAX(wrmax,wr(i&)) wrmin=MIN(wrmin,wr(i&)) wimax=MAX(wimax,ABS(wi(i&))) NEXT i& IF wrmax<0## THEN s\$="S" ELSE s\$="U" IF wrmax*wrmin<0## THEN s\$="X" IF wimax>0## THEN s\$=s\$+"F" ELSE s\$=s\$+"N" SELECT CASE s\$ 'Select colors for different stabilities CASE "SF": COLOR 7: GRAPHIC COLOR %RGB_WHITE CASE "SN": COLOR 15: GRAPHIC COLOR %RGB_BLACK CASE "UF": COLOR 13: GRAPHIC COLOR %RGB_MAGENTA CASE "UN": COLOR 11: GRAPHIC COLOR %RGB_CYAN CASE "XF": COLOR 7: GRAPHIC COLOR %RGB_GOLD CASE "XN": COLOR 12: GRAPHIC COLOR %RGB_RED END SELECT pr\$=FORMAT\$(a(1),"* ###.0000") pr\$=pr\$+FORMAT\$(a(2),"* ###.0000") FOR i&=1 TO ne& eigen\$=USING\$("* ##.#######",wr(i&)) IF ROUND(wi(i&),7)<>0## THEN eigen\$=eigen\$+" +"+USING\$("* ##.#######",wi(i&))+"i" END IF REPLACE"+ -" WITH "- " IN eigen\$ REPLACE"+ " WITH "+ " IN eigen\$ pr\$=pr\$+eigen\$ NEXT i& pr\$=pr\$+" "+s\$ rd&=0 'Radius of dot to plot rde&=MIN(rd&,xp&-b&-1,yp&-1,b&+w&-xp&-1,h&-yp&) IF rde&>0# THEN GRAPHIC ELLIPSE(xp&-rde&,yp&-rde&)-(xp&+rde&,yp&+rde&),,-1 ELSE GRAPHIC SET PIXEL(xp&,yp&) END IF NEXT xp& 'print pr\$ GRAPHIC REDRAW IF INKEY\$=\$ESC THEN EXIT FOR NEXT yp& x1&=b& 'Label the different regions in a contrasting color lastc&=%YELLOW 'A color not used FOR xp&=b& TO w&+b& GRAPHIC GET PIXEL(xp&,h&/2) TO c& IF c&<>cold& THEN x2&=xp& xpo&=x1&+(x2&-x1&)/2-6 IF x2&-x1&>25 AND cold&<>lastc& THEN GRAPHIC SET POS(xpo&,h&/2-12) SELECT CASE cold& CASE %RGB_WHITE: s\$="SF": GRAPHIC COLOR %BLACK,-2 CASE %RGB_BLACK: s\$="SN": GRAPHIC COLOR %WHITE,-2 CASE %RGB_MAGENTA: s\$="UF": GRAPHIC COLOR %WHITE,-2 CASE %RGB_CYAN: s\$="UN": GRAPHIC COLOR %BLACK,-2 CASE %RGB_GOLD: s\$="XF": GRAPHIC COLOR %BLACK,-2 CASE %RGB_RED: s\$="XN": GRAPHIC COLOR %WHITE,-2 END SELECT GRAPHIC PRINT s\$; lastc&=cold& END IF x1&=xp& cold&=c& END IF NEXT xp& GRAPHIC COLOR %BLACK FOR i&=1 TO xticks&-1 j&=b&+i&*w&/xticks& GRAPHIC LINE(j&,0)-(j&,h&/36) GRAPHIC LINE(j&,h&)-(j&,h&-h&/36) NEXT i& FOR i&=1 TO yticks&-1 j&=i&*h&/yticks& GRAPHIC LINE(b&,j&)-(b&+w&/36,j&) GRAPHIC LINE(b&+w&,j&)-(b&+w&-w&/36,j&) NEXT i& GRAPHIC REDRAW GRAPHIC SAVE EXE.NAME\$+".bmp" WAITKEY\$ GRAPHIC WINDOW END CASE"H" bail&=1e4 'Bailout condition (maximum # of iterations) dmax&=5 'Maximum number of digits of parameter resolution DESKTOP GET SIZE TO ddw&, dh& CONSOLE SET LOC w&+b&+6,0 GRAPHIC WINDOW EXE.FULL\$, 0, 0, w&+b&, h&+30 TO hWin??? GRAPHIC ATTACH hWin???, 0, REDRAW GRAPHIC COLOR %BLACK, %WHITE GRAPHIC CLEAR GRAPHIC SET FONT fhndl& xm\$=TRIM\$(STR\$(amax(2))) IF LEFT\$(xm\$,1)="." THEN xm\$="0"+xm\$ REPLACE"-." WITH"-0." IN xm\$ GRAPHIC TEXT SIZE xm\$ TO xmw, xmh GRAPHIC SET POS(b&-5-xmw,0) GRAPHIC PRINT xm\$ xm\$=ay\$ 'a(2) - Vertical axis GRAPHIC TEXT SIZE xm\$ TO xmw, xmh GRAPHIC SET POS(b&-5-xmw,h&/2-13) GRAPHIC PRINT xm\$ xm\$=TRIM\$(STR\$(amin(2))) REPLACE"-." WITH"-0." IN xm\$ IF LEFT\$(xm\$,1)="." THEN xm\$="0"+xm\$ GRAPHIC TEXT SIZE xm\$ TO xmw, xmh GRAPHIC SET POS(b&-5-xmw,h&-xmh) GRAPHIC PRINT xm\$ xm\$=TRIM\$(STR\$(amin(1))) IF LEFT\$(xm\$,1)="." THEN xm\$="0"+xm\$ REPLACE"-." WITH"-0." IN xm\$ GRAPHIC TEXT SIZE xm\$ TO xmw, xmh GRAPHIC SET POS(b&,h&) GRAPHIC PRINT xm\$ xm\$=ax\$ 'a(1) - Horizontal axis GRAPHIC TEXT SIZE xm\$ TO xmw, xmh GRAPHIC SET POS(b&+w&/2-xmw/2,h&) GRAPHIC PRINT xm\$ xm\$=TRIM\$(STR\$(amax(1))) IF LEFT\$(xm\$,1)="." THEN xm\$="0"+xm\$ REPLACE"-." WITH"-0." IN xm\$ GRAPHIC TEXT SIZE xm\$ TO xmw, xmh GRAPHIC SET POS(b&+w&-xmw,h&) GRAPHIC PRINT xm\$ GRAPHIC BOX (b&, 0) - (b&+w&, h&+1) 'GRAPHIC RENDER exe.name\$+".bmp",(0,0)-(w&+b&-1,h&+30-1) GRAPHIC REDRAW CONSOLE SET FOCUS CLS PRINT"Finding parameters that minimize |dxdt|..." PRINT q\$=" |dxdt| " FOR i&=1 TO np& q\$=q\$+SPACE\$(dmax&)+" a("+TRIM\$(STR\$(i&))+")" NEXT i& PRINT q\$ dxdtbest=xbound WHILE INKEY\$<>\$ESC FOR i&=1 TO ne& x(i&)=gauss NEXT i& yp=h&*RND xp=b&+w&*RND GRAPHIC GET PIXEL(xp,yp) TO c& IF c&<>%WHITE THEN ITERATE a(2)=ROUND(amin(2)+(amax(2)-amin(2))*(h&-yp)/h&,dmax&) a(1)=ROUND(amin(1)+(amax(1)-amin(1))*(xp-b&)/w&,dmax&) 'a(1)=MAX(a(1),0.5): a(2)=MAX(a(2),0.5) '?CONSTRAINT ON PARAMETERS dyn\$="": dxdtplot=xbound FOR n&=1 TO bail& CALL rk5(x(),dxdt(),ne&,h,x()) dxdtsq=0## FOR i&=1 TO ne& IF ABS(x(i&))>xbound THEN dyn\$="U" 'Unbounded solution dxdtsq=dxdtsq+dxdt(i&)^2 NEXT i& dxdt=SQR(dxdtsq) IF dxdt<1##/bail& THEN dyn\$="S" 'Stable equilibrium IF n&>bail&/2 THEN dxdtplot=MIN(dxdtplot,dxdt) IF n&>bail&/2## AND dxdtdxdtold THEN 'Success! dxdtbest=dxdtold PRINT USING\$(" ##.####^^^^",dxdtbest); FOR i&=1 TO np& PRINT USING\$("####."+STRING\$(dmax&,"#"),a(i&)); abest(i&)=a(i&) NEXT i& PRINT dyn\$="H" END IF dxdtold=dxdt IF n&=bail& THEN dyn\$="B" IF dyn\$<>"" THEN EXIT FOR NEXT n& IF dyn\$="U" OR dyn\$="B" THEN ITERATE IF dxdtplot>0.1## THEN ITERATE 'Don't plot this point (too far from equilibrium) rd&=3 'Radius of dot to plot rde&=MIN(rd&,xp-b&-1,yp-1,b&+w&-xp-1,h&-yp) IF rde&>0# THEN GRAPHIC ELLIPSE(xp-rde&,yp-rde&)-(xp+rde&,yp+rde&),,-1 ELSE GRAPHIC SET PIXEL(xp,yp) END IF GRAPHIC REDRAW ti=TIMER IF titold+300 THEN 'Save image every 300 seconds GRAPHIC SAVE EXE.NAME\$+".bmp" told=ti END IF WEND CASE"A","L","K" lemin=0.001## 'Minimum LE to qualify as chaotic prec=0.1## 'Relative precision in calculated LE dkymin=2.01## 'Minimun DKY to qualify as chaotic (starting value) bail&=1e7 'Bailout condition (maximum # of iterations) eps=0.1## 'Size of search space IF q\$="A" THEN what\$="DMIN" IF q\$="L" THEN what\$="LE" IF q\$="K" THEN what\$="DKY" ff&=FREEFILE 'Start over with new search.txt file OPEN EXE.NAME\$+".txt" FOR OUTPUT AS ff& CLOSE ff& DESKTOP GET SIZE TO ddw&, dh& CONSOLE SET LOC w&+b&+6,0 CLS IF what\$="DKY" AND ne&>3 THEN PRINT "Warning: Cannot reliable Dky in";ne&;"dimensions!" GRAPHIC WINDOW EXE.FULL\$, 0, 0, w&+b&, h&+24 TO hWin??? GRAPHIC ATTACH hWin???, 0, REDRAW GRAPHIC COLOR %BLACK, %WHITE CONSOLE SET FOCUS FOR i&=1 TO ne& xobest(i&)=gauss NEXT i& FOR i&=1 TO np& abest(i&)=(amin(i&)+amax(i&))/2## d&(i&)=dmax& NEXT i& trials&=0 dbest=0## DO eps=0.99999*eps 'Slowly shrink size of search space INCR trials& FOR i&=1 TO np& IF np&*RND<4 THEN di&=MIN(dmax&,d&(i&)+1) ELSE di&=d&(i&) a(i&)=ROUND(abest(i&)+eps*MAX(1,ABS(abest(i&)))*(gauss-nudge*SGN(CINT(abest(i&)))),di&) 'a(i&)=ABS(a(i&)) '?ONLY POSITIVE PARAMETERS!!!!!!!!!!!!!! NEXT i& dif=0##: FOR i&=1 TO np&: dif=dif+ABS(abest(i&)-a(i&)): NEXT i& IF dif=0## THEN eps=10##*eps: ITERATE 'Neighborhood shrank too much! prec=MIN(prec,SQR(eps)) 'Improve precsion when neighborhood is small RESET xe(),dx(),lehist() FOR j&=1 TO ne& xo(j&)=ROUND(xobest(j&)+eps*gauss,2) 'Choose a random initial condition 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: xom=1e6: xtot=0##: xabs=0##: ytot=0##: yabs=0##: ztot=0##: zabs=0## dmin=xbound DO 'n&=1: Orbits are separating too much 'n&=2: LE is too small or negative 'n&=3: Phase space is expanding 'n&=4: Orbits have collapsed together 'n&=5: Bailout condition reached 'n&=6: Orbit is unbounded 'n&=7: Stable fixed point 'n&=8: Attractor is too small 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 n&=6: EXIT LOOP 'Orbit is unbounded dx(i&)=xe(i&)-x(i&) rs=rs+dx(i&)*dx(i&) NEXT i& IF rs>dr THEN n&=1: EXIT LOOP 'Orbits are separating too much 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 dmin=MIN(dmin,ABS(x(1))) ltot=ltot+LOG(rs) INCR ntot& xtot=xtot+x(1): xabs=xabs+ABS(x(1)): asym=xtot/xabs 'Calculate asymmetry IF ne&>1 THEN ytot=ytot+x(2): yabs=yabs+ABS(x(2)): asym=asym+ytot/yabs 'if ne&>2 then ztot=ztot+x(3): zabs=zabs+ABS(x(3)): asym=asym+ztot/zabs 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 x1hist(ntot& MOD nm&)=x(1) lehist(ntot& MOD nm&)=le IF ntot&>=nm& AND (ntot& MOD nm&)=nm&-1 THEN'See if LE has converged CALL lyapext(x1hist(),x1,dx1) IF dx1<0.1 THEN n&=8: EXIT LOOP 'Attractor is too small CALL lyapext(lehist(),le,dle) IF lesum-dle>ABS(lemin) THEN n&=3: EXIT LOOP 'Phase space is expanding dky=ne&-1+le/(le-lesum) ddky=ABS(dle/(le-lesum)) IF what\$="LE" THEN IF le+dle/preclemin THEN 'Chaos! lemin=le-dle n&=0 EXIT LOOP END IF END IF IF what\$="DKY" THEN 'Only valid for ne&<4! IF dky+ddky/precdkymin THEN 'Chaos dkymin=dky-ddky n&=0 EXIT LOOP END IF END IF IF what\$="DMIN" THEN IF le+dle/preclemin THEN 'Chaos! n&=9 IF dmin>dbest AND ABS(asym)>0.1## THEN dbest=dmin: n&=0 EXIT LOOP END IF END IF END IF END IF ELSE n&=4: EXIT LOOP 'Orbits have collapsed together END IF IF n&>bail& THEN n&=5: EXIT LOOP 'Bailout condition reached CALL derivs(x(),dxdt(),ne&) dxdt=0## FOR i&=1 TO ne& dxdt=dxdt+ABS(dxdt(i&))/(1+ABS(x(i&))) NEXT i& IF dxdt<1e-4 THEN n&=7: EXIT LOOP 'Stable fixed point LOOP IF INKEY\$=\$ESC THEN EXIT LOOP IF n&>0 THEN 'Solution is not interesting COLOR 14 PRINT TRIM\$(STR\$(n&)); ELSE 'An improved chaotic case was found eps=1.1##*eps 'Increase size of search space PRINT CALL plot(w&,h&,b&,x(),dxdt(),ne&,h,0,0) IF ABS(x(1))>xbound OR ABS(x(2))>xbound THEN EXIT IF 'Graph was off scale BEEP INCR cases& GRAPHIC SAVE EXE.NAME\$+".bmp" z\$=STRING\$(dmax&,"0") CALL best(ebest&,xo(),le,dle,lesum,ntot&,z\$) COLOR 15 PRINT FORMAT\$(cases&,"* #"); FOR i&=1 TO np& PRINT FORMAT\$(a(i&)," 0."+z\$+"; -0."+z\$); NEXT i& PRINT FORMAT\$(le," 0.0000; -0.0000");CHR\$(177); PRINT FORMAT\$(dle,"0.0000"); IF ne&=3 THEN PRINT " 0"; PRINT FORMAT\$(lesum-le," 0.0000; -0.0000"); PRINT FORMAT\$(dky," 0.0000 "); ELSE PRINT FORMAT\$(lesum," 0.0000 ; -0.0000 "); END IF FOR i&=1 TO np&: abest(i&)=a(i&): NEXT i& FOR i&=1 TO ne&: xobest(i&)=xo(i&): PRINT xo(i&);: NEXT i& PRINT ntot&+1;" eps =";CSNG(eps);" asym =";CSNG(asym);" dmin =";CSNG(dmin) END IF LOOP PRINT PRINT"Trials =";trials& GRAPHIC REDRAW GRAPHIC SAVE EXE.NAME\$+".bmp" WAITKEY\$ GRAPHIC WINDOW END CASE ELSE BEEP END SELECT LOOP END FUNCTION SUB map (x(), f(), n&) 'Replace SUB DERIVS with this for map 'Returns the map function f-x(i&) of x(i&) f(1) = -3.5##*x(1)*x(1) + x(2) f(2) = 0.3## - ABS(x(1)) 'Don't mess with the following lines! FOR i&=1 TO n& 'For iterated map f(i&)=f(i&)-x(i&) NEXT i& END SUB SUB DERIVS (x(), dxdt(), n&) 'Returns the time derivatives dxdt(i&) of x(i&) dxdt(1) = a(1)*x(1) + a(2)*x(3) dxdt(2) = x(1)*x(3) - x(2) dxdt(3) = -x(1) + x(2) END SUB SUB DERIVALL(x(), dxdt(), nn&) n&=SQR(1+4##*nn&)/2## - 0.5## CALL derivs(x(), dxdt(), n&) IF n&<4 THEN EXIT SUB 'Cannot be hyperchaotic! ' Linearized equations: DIM jacob(n&,n&) CALL jacobian(x(),dxdt(),n&,1e-8##,jacob(),0.01##) 'jacob(1,1)=-1##: jacob(2,1)=1## 'We know the Jacobian! 'jacob(1,2)=-x(3)*x(3)*x(3): jacob(3,2)=-3##*x(1)*x(3)*x(3): jacob(4,2)=1## 'jacob(1,3)=x(2)*x(2)*x(2): jacob(2,3)=3##*x(1)*x(2)*x(2) 'jacob(2,4)=-a(2) FOR i&=1 TO n& FOR j&=1 TO n& xdot=0## FOR k&=1 TO n& xdot=xdot+jacob(k&,j&)*x(n&*k&+i&) NEXT k& dxdt(n&*j&+i&)=xdot NEXT j& NEXT i& END SUB SUB RK5 (x(), dxdt(), n&, h, xnew()) 'Calls RK4 (4th order Runge-Kutta) repeatedly with adjustable step size eallow=1e-8 '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&,100) '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&)) NEXT i& 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 DERIVALL(X(), DXDT(), N&) FOR I& = 1 TO N& XT(I&) = X(I&) + HH * DXDT(I&) NEXT I& CALL DERIVALL(XT(), DXT(), N&) FOR I& = 1 TO N& XT(I&) = X(I&) + HH * DXT(I&) NEXT I& CALL DERIVALL(XT(), DXM(), N&) FOR I& = 1 TO N& XT(I&) = X(I&) + H * DXM(I&) DXM(I&) = DXT(I&) + DXM(I&) NEXT I& CALL DERIVALL(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 FUNCTION gauss 'Returns a normally (Gaussian) distributed random deviate 'with mean of zero and standard deviation of 1.0 DO V1 = 2## * RND - 1## V2 = 2## * RND - 1## R = V1 * V1 + V2 * V2 LOOP WHILE R >= 1## OR R = 0## GAUSS = V2 * SQR(-2## * LOG(R) / R) END FUNCTION 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)0## THEN q1\$=q1\$+STR\$(i&) 'NEXT i& 'q1\$=q1\$+SPACE\$(MAX(0,2*np&-LEN(q1\$)))+FORMAT\$(ebest&,"* #") FOR i&=1 TO np& q1\$=q1\$+FORMAT\$(a(i&)," 0."+z\$+"; -0."+z\$) NEXT i& found&=0 IF ISFILE(EXE.NAME\$+".txt") THEN ff&=FREEFILE OPEN EXE.NAME\$+".txt" FOR INPUT AS ff& WHILE NOT EOF(ff&) LINE INPUT #ff&, q2\$ IF LEFT\$(q2\$,LEN(q1\$))=q1\$ THEN found&=-1 WEND CLOSE ff& END IF IF found& THEN EXIT SUB 'We already recorded this one OPEN EXE.NAME\$+".txt" FOR APPEND AS ff& PRINT #ff&, q1\$; PRINT #ff&, FORMAT\$(le," 0.0000; -0.0000");CHR\$(177); PRINT #ff&, FORMAT\$(dle,"0.0000"); IF ne&=3 THEN PRINT #ff&, " 0"; PRINT #ff&, FORMAT\$(lesum-le," 0.0000; -0.0000"); PRINT #ff&, FORMAT\$(2+le/(le-lesum)," 0.0000"); ELSE PRINT #ff&, FORMAT\$(lesum," 0.0000; -0.0000"); END IF PRINT #ff&, ntot&+1; FOR i&=1 TO ne&: PRINT #ff&, xo(i&);: NEXT i& IF maxdots& THEN PRINT #ff&, " dots =";dots& ELSE PRINT #ff&,"" CLOSE ff& END SUB SUB plot(w&,h&,b&,x(),dxdt(),ne&,h,wo&,ho&) 'Plots state-space trajectory np&=w&*h&/3 'Number of iterates to plot os=0.05## 'Amount of overscan xm1=1e6: ym1=xm1: xm2=-xm1: ym2=xm2 DIM xo(ne&) htry=h FOR n&=-np& TO np& 'Find tentative range of values CALL rk5(x(),dxdt(),ne&,htry,x()) IF n&<1 THEN ITERATE 'Not on attractor yet xm1=MIN(xm1,x(1)): ym1=MIN(ym1,x(2)) xm2=MAX(xm2,x(1)): ym2=MAX(ym2,x(2)) NEXT n& xmo&=1+INT(xm2-xm1): ymo&=1+INT(ym2-ym1) FOR i&=1 TO ne&: xo(i&)=x(i&): NEXT i& xm1=1e6: ym1=xm1: xm2=-xm1: ym2=xm2: zmin=1e6: zmax=-zmin FOR n&=1 TO np& 'Find exact range of values CALL rk5(x(),dxdt(),ne&,htry,x()) xm1=MIN(xm1,x(1)): ym1=MIN(ym1,x(2)) xm2=MAX(xm2,x(1)): ym2=MAX(ym2,x(2)) IF ne&>2 THEN zmin=MIN(zmin,x(3)): zmax=MAX(zmax,x(3)) v=ABS(2*dxdt(1)*w&/xmo&)+ABS(2*dxdt(2)*h&/ymo&) htry=MIN(1/v,h) 'Move at 1 pixel per time step NEXT n& xm1&=INT(xm1-os*(xm2-xm1)): ym1&=INT(ym1-os*(ym2-ym1)) xm2&=CEIL(xm2+os*(xm2-xm1)): ym2&=CEIL(ym2+os*(ym2-ym1)) xm&=xm2&-xm1&: ym&=ym2&-ym1& 'xm&=MAX(xm&,ym&): ym&=MAX(xm&,ym&) '?FORCE SCALES TO BE THE SAME IF xm&=0 OR ym&=0 THEN x(1)=1e99: EXIT SUB 'A variable has collapsed GRAPHIC SET FONT fhndl& GRAPHIC CLEAR SELECT CASE xm& CASE 0 TO 2: xticks&=10 CASE 3 TO 20: xticks&=xm& CASE 21 TO 40: xticks&=CEIL(xm&/2): xm&=2*xticks& CASE 41 TO 100: xticks&=CEIL(xm&/5): xm&=5*xticks& CASE 101 TO 1000: xticks&=CEIL(xm&/10): xm&=10*xticks& CASE > 1000: xticks&=100: xm&=100*CEIL(xm&/100) END SELECT SELECT CASE ym& CASE 0 TO 2: yticks&=10 CASE 3 TO 20: yticks&=ym& CASE 21 TO 40: yticks&=CEIL(ym&/2): ym&=2*yticks& CASE 41 TO 100: yticks&=CEIL(ym&/5): ym&=5*yticks& CASE 101 TO 1000: yticks&=CEIL(ym&/10): ym&=10*yticks& CASE > 1000: yticks&=100: ym&=100*CEIL(ym&/100) END SELECT xm1&=(xm1&+xm2&-xm&)/2: ym1&=(ym1&+ym2&-ym&)/2 xm2&=xm1&+xm&: ym2&=ym1&+ym& xm\$=TRIM\$(STR\$(ym2&)) GRAPHIC TEXT SIZE xm\$ TO xmw, xmh GRAPHIC SET POS(MAX(wo&+b&-5-xmw,0),ho&) GRAPHIC PRINT xm\$ xm\$="Y" '?V OR Y GRAPHIC TEXT SIZE xm\$ TO xmw, xmh GRAPHIC SET POS(MAX(wo&+b&-5-xmw,0),ho&+h&/2-xmh/2) GRAPHIC PRINT xm\$ xm\$=TRIM\$(STR\$(ym1&)) GRAPHIC TEXT SIZE xm\$ TO xmw, xmh GRAPHIC SET POS(MAX(wo&+b&-5-xmw,0),ho&+h&-20) GRAPHIC PRINT xm\$ xm\$=TRIM\$(STR\$(xm1&)) GRAPHIC TEXT SIZE xm\$ TO xmw, xmh GRAPHIC SET POS(wo&+b&,ho&+h&) GRAPHIC PRINT xm\$ xm\$="X" GRAPHIC TEXT SIZE xm\$ TO xmw, xmh GRAPHIC SET POS(wo&+b&+w&/2-xmw/2,ho&+h&) GRAPHIC PRINT xm\$ xm\$=TRIM\$(STR\$(xm2&)) GRAPHIC TEXT SIZE xm\$ TO xmw, xmh GRAPHIC SET POS(wo&+b&+w&-xmw,ho&+h&) GRAPHIC PRINT xm\$ FOR i&=1 TO xticks&-1 j&=wo&+b&+i&*w&/xticks& GRAPHIC LINE(j&,ho&)-(j&,ho&+h&/36) GRAPHIC LINE(j&,ho&+h&)-(j&,ho&+h&-h&/36) NEXT i& FOR i&=1 TO yticks&-1 j&=ho&+i&*h&/yticks& GRAPHIC LINE(wo&+b&,j&)-(wo&+b&+w&/36,j&) GRAPHIC LINE(wo&+b&+w&-1,j&)-(wo&+b&+w&-w&/36-1,j&) NEXT i& GRAPHIC BOX (wo&+b&, ho&) - (wo&+b&+w&, ho&+h&+1) FOR i&=1 TO ne&: x(i&)=xo(i&): NEXT i& htry=h dots&=0 FOR n&=1 TO np& CALL rk5(x(),dxdt(),ne&,htry,x()) IF n&<1 THEN ITERATE CALL project(x(),xm1&,xm2&,ym1&,ym2&,zmin,zmax,b&,h&,w&,xp,yp,zp,"plane") IF xp2 THEN zp=255*(x(3)-zmin)/(zmax-zmin) ELSE zp=255 dx=xp-INT(xp) 'Antialias dy=yp-INT(yp) g1=180-zp*dx*dy g2=180-zp*dx*(1-dy) g3=180-zp*(1-dx)*dy g4=180-zp*(1-dx)*(1-dy) c1&=RGB(g1,g1,g1) c2&=RGB(g2,g2,g2) c3&=RGB(g3,g3,g3) c4&=RGB(g4,g4,g4) GRAPHIC GET PIXEL(wo&+INT(xp)+1,ho&+INT(yp)+1) TO b1& IF c1&2 THEN 'Make shadow gs=0: FOR i&=1 TO 12: gs=gs+RND: NEXT i& 'Approximation to a Gaussian xp=xp+(w&/600)*zp/gs gs=0: FOR i&=1 TO 12: gs=gs+RND: NEXT i& 'Approximation to a Gaussian yp=yp+(h&/600)*zp/gs GRAPHIC GET PIXEL(wo&+xp,ho&+yp) TO b1& IF b1&=%WHITE AND xp0## AND dx<8## dx=10##*dx WEND ticks&=0 FOR i&=1 TO maxticks& IF ABS(i&*CINT(dx/i&)-dx)<0.001 THEN ticks&=i& NEXT i& END FUNCTION SUB mouseloc(w&,h&,b&,amin(),amax(),ax\$,ay\$) 'prints the position of the graphic mouse when it's clicked TXT.END 'Close any previos text window GRAPHIC WINDOW CLICK TO click&, x!, y! 'Mouse position IF click&=0 THEN EXIT SUB GRAPHIC GET LOC TO xw&, yw& 'Upper-left corner of graphic window a1!=amin(1)+(amax(1)-amin(1))*(x!-b&)/w& a2!=amin(2)+(amax(2)-amin(2))*(h&-y!)/h& a\$=ax\$+" = "+TRIM\$(STR\$(a1!))+", "+ay\$+" = "+TRIM\$(STR\$(a2!)) TXT.WINDOW("", x!+xw&+10, y!+yw&+10,1,LEN(a\$)) TO hWin??? TXT.PRINT a\$; t=TIMER 'Timeout after 10 seconds DO GRAPHIC WINDOW CLICK TO click&, x!, y! LOOP UNTIL click& OR TIMER>t+10 END SUB SUB jacobian(x(),dxdt(),n&,dr,jacob(),h) 'Calculates Jacobian jacob(i&,j&) of derivs at current x() 'Takes numerical partial derivatives by perturbing +/- dr/2 FOR i&=1 TO n& xsave=x(i&) x(i&)=x(i&)-0.5##*dr IF h=1## THEN CALL map(x(),dxdt(),n&) ELSE CALL derivs(x(),dxdt(),n&) FOR j&=1 TO n& jacob(i&,j&)=dxdt(j&) NEXT j& x(i&)=x(i&)+dr IF h=1## THEN CALL map(x(),dxdt(),n&) ELSE CALL derivs(x(),dxdt(),n&) FOR j&=1 TO n& jacob(i&,j&)=(dxdt(j&)-jacob(i&,j&))/dr NEXT j& x(i&)=xsave NEXT i& END SUB FUNCTION det(d()) 'Calculates the determinant of the square matrix d() n& = UBOUND(d,1) DIM a(1 TO n&*n&) IF n& <> UBOUND(d,2) THEN PRINT "Array must be square to have a determinant." FUNCTION = 0## EXIT FUNCTION END IF s = 1## DO b = 0## FOR i& = 1 TO n& FOR j& = 1 TO n& aa = ABS(d(i&,j&)) IF aa < b THEN ITERATE p& = i& q& = j& b = aa NEXT j& NEXT i& w = 1## + 4##*INT((p&+q&)/2##) - 2##*(p&+q&) s = s*d(p&,q&)*w IF s = 0## THEN EXIT LOOP v& = 1 FOR i& = 1 TO n& IF i& = p& THEN ITERATE FOR j& = 1 TO n& IF j& = q& THEN ITERATE d(i&,j&) = d(i&,j&) - d(p&,j&)/d(p&,q&)*d(i&,q&) a(v&) = d(i&,j&) v& = v& + 1 NEXT j& NEXT i& n& = n& - 1 IF n& = 0 THEN EXIT LOOP v& = 1 FOR i& = 1 TO n& FOR j& = 1 TO n& d(i&,j&) = a(v&) v& = v& + 1 NEXT j& NEXT i& LOOP FUNCTION = s END FUNCTION SUB balanc (A(), N&) 'Balances the Jacobian matrix for improved calculation of eigenvalues RADIX = 2## SQRDX = 4## DO LAST& = 1 FOR I& = 1 TO N& C = 0## R = 0## FOR J& = 1 TO N& IF J& <> I& THEN C = C + ABS(A(J&, I&)) R = R + ABS(A(I&, J&)) END IF NEXT J& IF C <> 0## AND R <> 0## THEN G = R / RADIX F = 1## S = C + R WHILE C < G F = F * RADIX C = C * SQRDX WEND G = R * RADIX WHILE C > G F = F / RADIX C = C / SQRDX WEND IF (C + R) / F < 0.95## * S THEN LAST = 0## G = 1## / F FOR J& = 1 TO N& A(I&, J&) = A(I&, J&) * G NEXT J& FOR J& = 1 TO N& A(J&, I&) = A(J&, I&) * F NEXT J& END IF END IF NEXT I& LOOP WHILE LAST& = 0 END SUB SUB elmhes (A(), N&) 'Converts an N& x N& matrix A() to Hessenberg form FOR M& = 2 TO N& - 1 X = 0## I& = M& FOR J& = M& TO N& IF ABS(A(J&, M& - 1)) > ABS(X) THEN X = A(J&, M& - 1) I& = J& END IF NEXT J& IF I& <> M& THEN FOR J& = M& - 1 TO N& Y = A(I&, J&) A(I&, J&) = A(M&, J&) A(M&, J&) = Y NEXT J& FOR J& = 1 TO N& Y = A(J&, I&) A(J&, I&) = A(J&, M&) A(J&, M&) = Y NEXT J& END IF IF X <> 0## THEN FOR I& = M& + 1 TO N& Y = A(I&, M& - 1) IF Y <> 0## THEN Y = Y / X A(I&, M& - 1) = Y FOR J& = M& TO N& A(I&, J&) = A(I&, J&) - Y * A(M&, J&) NEXT J& FOR J& = 1 TO N& A(J&, M&) = A(J&, M&) + Y * A(J&, I&) NEXT J& END IF NEXT I& END IF NEXT M& END SUB SUB hqr (A(), N&, WR(), WI()) 'Calculates eigenvalues of Hessenberg matrix A() ANORM = ABS(A(1, 1)) FOR I& = 2 TO N& FOR J& = I& - 1 TO N& ANORM = ANORM + ABS(A(I&, J&)) NEXT J& NEXT I& NN& = N& T = 0## WHILE NN& >= 1 ITS& = 0 DO DONE& = -1 FOR L& = NN& TO 2 STEP -1 S = ABS(A(L& - 1, L& - 1)) + ABS(A(L&, L&)) IF S = 0## THEN S = ANORM IF ABS(A(L&, L& - 1)) + S = S THEN EXIT FOR NEXT L& IF ABS(A(L&, L& - 1)) + S <> S THEN L& = 1 X = A(NN&, NN&) IF L& = NN& THEN WR(NN&) = X + T WI(NN&) = 0## NN& = NN& - 1 ELSE Y = A(NN& - 1, NN& - 1) W = A(NN&, NN& - 1) * A(NN& - 1, NN&) IF L& = NN& - 1 THEN P = 0.5## * (Y - X) Q = P ^ 2 + W Z = SQR(ABS(Q)) X = X + T IF Q >= 0## THEN Z = P + ABS(Z) * SGN(P) WR(NN&) = X + Z WR(NN& - 1) = WR(NN&) IF Z <> 0## THEN WR(NN&) = X - W / Z WI(NN&) = 0## WI(NN& - 1) = 0## ELSE WR(NN&) = X + P WR(NN& - 1) = WR(NN&) WI(NN&) = Z WI(NN& - 1) = -Z END IF NN& = NN& - 2 ELSE IF ITS& = 30 THEN PRINT "Too many iterations": EXIT SUB IF ITS& = 10 OR ITS& = 20 THEN T = T + X FOR I& = 1 TO NN& A(I&, I&) = A(I&, I&) - X NEXT I& S = ABS(A(NN&, NN& - 1)) + ABS(A(NN& - 1, NN& - 2)) X = 0.75## * S Y = X W = -0.4375## * S ^ 2 END IF ITS& = ITS& + 1 FOR M& = NN& - 2 TO L& STEP -1 Z = A(M&, M&) R = X - Z S = Y - Z P = (R * S - W) / A(M& + 1, M&) + A(M&, M& + 1) Q = A(M& + 1, M& + 1) - Z - R - S R = A(M& + 2, M& + 1) S = ABS(P) + ABS(Q) + ABS(R) P = P / S Q = Q / S R = R / S IF M& = L& THEN EXIT FOR U = ABS(A(M&, M& - 1)) * (ABS(Q) + ABS(R)) V = ABS(P) * (ABS(A(M& - 1, M& - 1)) + ABS(Z) + ABS(A(M& + 1, M& + 1))) IF U + V = V THEN EXIT FOR NEXT M& FOR I& = M& + 2 TO NN& A(I&, I& - 2) = 0## IF I& <> M& + 2 THEN A(I&, I& - 3) = 0## NEXT I& FOR K& = M& TO NN& - 1 IF K& <> M& THEN P = A(K&, K& - 1) Q = A(K& + 1, K& - 1) R = 0## IF K& <> NN& - 1 THEN R = A(K& + 2, K& - 1) X = ABS(P) + ABS(Q) + ABS(R) IF X <> 0## THEN P = P / X Q = Q / X R = R / X END IF END IF S = P ^ 2 + Q ^ 2 + R ^2 S = SQR(S) * SGN(P) IF S <> 0## THEN IF K& = M& THEN IF L& <> M& THEN A(K&, K& - 1) = -A(K&, K& - 1) ELSE A(K&, K& - 1) = -S * X END IF P = P + S X = P / S Y = Q / S Z = R / S Q = Q / P R = R / P FOR J& = K& TO NN& P = A(K&, J&) + Q * A(K& + 1, J&) IF K& <> NN& - 1 THEN P = P + R * A(K& + 2, J&) A(K& + 2, J&) = A(K& + 2, J&) - P * Z END IF A(K& + 1, J&) = A(K& + 1, J&) - P * Y A(K&, J&) = A(K&, J&) - P * X NEXT J& NTMP& = NN& IF K& + 3 < NN& THEN NTMP& = K& + 3 FOR I& = L& TO NTMP& P = X * A(I&, K&) + Y * A(I&, K& + 1) IF K& <> NN& - 1 THEN P = P + Z * A(I&, K& + 2) A(I&, K& + 2) = A(I&, K& + 2) - P * R END IF A(I&, K& + 1) = A(I&, K& + 1) - P * Q A(I&, K&) = A(I&, K&) - P NEXT I& END IF NEXT K& DONE& = 0 END IF END IF LOOP WHILE NOT DONE& WEND END SUB