'Program VERIFY.BAS (c) 2008 by J. C. Sprott 'Verifies, plots, and calculates LEs and optimal initial conditions for a system of ODEs 'Use PowerBasic Console Compiler #DEBUG ERROR OFF DEFEXT a-z FUNCTION PBMAIN() ne&=4 'Number of equations dr=1e-8## 'Separation of orbits h=.01## 'Iteration step size nm&=1e4 'Number of iterates per test w&=800 'Width of plot h&=800 'Height of plot b&=40 'Border size dr2=dr*dr DIM x(ne&),dxdt(ne&) x(1)=0.9##: x(2)=2.4##: x(3)=-1.8##: x(4)=4## 'Initial conditions DIM pi AS GLOBAL EXT pi=3.141592653589794## zs=0## 'Take the Poincare section at this value of x(3) CONSOLE NAME EXE.NAME\$ DO CONSOLE SET FOCUS CLS PRINT"Press G for graph" PRINT" L for Lyapunov exponent" PRINT" P for Poincare section" PRINT" S for chaotic sea" PRINT" I for initial condition" PRINT" T for time series of Poincare plot" PRINT" B for basin of attractor" PRINT" W for spatiotemporal plot" PRINT" V for values of the variables" PRINT" A for average values and " PRINT" M for map data" PRINT" X for time series x(t)" PRINT" to end program: "; WHILE q\$="" q\$=UCASE\$(INKEY\$) IF q\$=\$ESC THEN EXIT FUNCTION IF INSTR("GLPSITBWVAMX",q\$)=0 THEN EXIT LOOP WEND SELECT CASE q\$ CASE "G" 'Graph attractor DESKTOP GET SIZE TO ddw&, dh& CONSOLE SET LOC w&+b&+6,0 CLS GRAPHIC WINDOW EXE.NAME\$, 0, 0, w&+b&, h&+24*SGN(b&) TO hWin??? GRAPHIC ATTACH hWin???, 0, REDRAW GRAPHIC COLOR %BLACK, %WHITE GRAPHIC SET FOCUS: GRAPHIC REDRAW: CONSOLE SET FOCUS PRINT"Press to save image to "+EXE.NAME\$+".bmp: "; WHILE INKEY\$<>\$ESC GRAPHIC CLEAR CALL plot(w&,h&,b&,x(),dxdt(),ne&,h,0,0) WEND BEEP GRAPHIC SAVE EXE.NAME\$+".bmp" SHELL"bmp2eps "+EXE.NAME\$ CALL initcond(x(),dxdt(),ne&,h) CASE "L" 'Calculate Lyapunov exponents CLS REDIM xe(ne&),dx(ne&),lehist(nm&-1) FOR j&=1 TO ne& xe(j&)=x(j&)+dr/SQR(ne&) 'Perturb nearby initial condition NEXT j& n&=0: ltot=0##: tr=0##: ntot&=0 DO INCR n& CALL rk5(x(),dxdt(),ne&,h,x()) CALL rk5(xe(),dxdt(),ne&,h,xe()) rs=0## FOR i&=1 TO ne& 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&) IF ABS(x(i&))>1e11 THEN PRINT"Orbit escaped: x(";TRIM\$(STR\$(i&));") =";x(i&): WAITKEY\$: EXIT LOOP NEXT i& 'IF rs>dr THEN PRINT"Orbits separating too much: rs =";rs: WAITKEY\$: 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 INKEY\$=\$ESC THEN EXIT SELECT PRINT FORMAT\$(n&+2,"* ########"); PRINT FORMAT\$(le," 0.0000; -0.0000");CHR\$(177);FORMAT\$(dle,"0.0000"); IF ne&=3 THEN PRINT " 0"; PRINT FORMAT\$(lesum-le," 0.0000; -0.0000"); PRINT FORMAT\$(2+le/(le-lesum)," 0.0000"); ELSE PRINT FORMAT\$(lesum," 0.0000; -0.0000"); END IF m=1.4##: q=4##: p=0##: mu=4##: k=1##: e=(q-1##)*mu/(q*q): emu=e/mu e=(q-1##)*mu/q^2 ham=0.5##*m*x(2)*x(2)+m*(p-x(4))^2/(8##*mu*mu*x(1)*x(1)) ham=ham+m*e*e*x(1)*x(1)/8##+0.25##*m*e*(p+x(4))/mu ham=ham+(1##-m/q)*(q-2##)*x(4)/(2##*q)+k/((q-1##)*x(1)) ham=ham+(1##-m/q)*(x(2)*SIN(2##*mu*x(3))-0.5##*(e*x(1)+(p-x(4))/(mu*x(1)))*COS(2##*mu*x(3)))*SQR(x(4)) PRINT" ";ham/2.52415086202632## '?CHANGE THE NORMALIZATION END IF END IF ELSE n&=0: PRINT"Orbits collapsed together: rs=";rs: WAITKEY\$: EXIT LOOP END IF IF n&>2e9 THEN PRINT"Done!": WAITKEY\$: EXIT LOOP 'Stop before n& overflows LOOP CASE "P" 'Make Poincare section for attractor DESKTOP GET SIZE TO ddw&, dh& CONSOLE SET LOC w&+b&+6,0 CLS GRAPHIC WINDOW EXE.NAME\$, 0, 0, w&+b&, h&+b&/2 TO hWin??? GRAPHIC ATTACH hWin???, 0, REDRAW GRAPHIC COLOR %BLACK, %WHITE GRAPHIC SET FOCUS: GRAPHIC REDRAW: CONSOLE SET FOCUS PRINT"Press to save image to "+EXE.NAME\$+".bmp: "; GRAPHIC CLEAR CALL poincare(w&,h&,b&,x(),dxdt(),ne&,h,zs,0,0) GRAPHIC SAVE EXE.NAME\$+".bmp" SHELL"bmp2eps "+EXE.NAME\$ WAITKEY\$ CASE "S" 'Make Poincare section for chaotic sea DESKTOP GET SIZE TO ddw&, dh& CONSOLE SET LOC w&+b&+6,0 CLS GRAPHIC WINDOW EXE.NAME\$, 0, 0, w&+b&, h&+b&/2 TO hWin??? GRAPHIC ATTACH hWin???, 0, REDRAW GRAPHIC COLOR %BLACK, %WHITE GRAPHIC SET FOCUS: GRAPHIC REDRAW: CONSOLE SET FOCUS PRINT"Press to save image to "+EXE.NAME\$+".bmp: "; GRAPHIC CLEAR CALL sea(w&,h&,b&,x(),dxdt(),ne&,h,zs,0,0) GRAPHIC SAVE EXE.NAME\$+".bmp" SHELL"bmp2eps "+EXE.NAME\$ WAITKEY\$ CASE "I" CALL initcond(x(),dxdt(),ne&,h) CASE "T" 'Make Bifurcation diagram versus time DESKTOP GET SIZE TO ddw&, dh& CONSOLE SET LOC w&+b&+6,0 CLS GRAPHIC WINDOW EXE.NAME\$, 0, 0, w&+b&, b&/8+h&+b&/2 TO hWin??? GRAPHIC ATTACH hWin???, 0, REDRAW GRAPHIC COLOR %BLACK, %WHITE GRAPHIC SET FOCUS: GRAPHIC REDRAW: CONSOLE SET FOCUS PRINT"Press to save image to "+EXE.NAME\$+".bmp: "; CALL bifurcation(w&,h&,b&,x(),dxdt(),ne&,h,zs) GRAPHIC SAVE EXE.NAME\$+".bmp" SHELL"bmp2eps "+EXE.NAME\$ CASE "B" 'Plot basin of attraction DESKTOP GET SIZE TO ddw&, dh& CONSOLE SET LOC w&+b&+6,0 CLS GRAPHIC WINDOW EXE.NAME\$, 0, 0, w&+b&, h&+b&/2 TO hWin??? GRAPHIC ATTACH hWin???, 0, REDRAW GRAPHIC COLOR %BLACK, %WHITE GRAPHIC SET FOCUS: GRAPHIC REDRAW: CONSOLE SET FOCUS PRINT"Press to save image to "+EXE.NAME\$+".bmp: "; GRAPHIC CLEAR CALL basin(w&,h&,b&,x(),dxdt(),ne&,h,zs,0,0) GRAPHIC SAVE EXE.NAME\$+".bmp" SHELL"bmp2eps "+EXE.NAME\$ WAITKEY\$ CASE "W" 'Make spatiotemporal plot DESKTOP GET SIZE TO ddw&, dh& CONSOLE SET LOC w&+b&+6,0 CLS GRAPHIC WINDOW EXE.NAME\$, 0, 0, w&+b&, b&/8+h&+b&/2 TO hWin??? GRAPHIC ATTACH hWin???, 0, REDRAW GRAPHIC COLOR %BLACK, %WHITE GRAPHIC SET FOCUS: GRAPHIC REDRAW: CONSOLE SET FOCUS PRINT"Press to save image to "+EXE.NAME\$+".bmp: "; GRAPHIC CLEAR CALL spatiotemporal(w&,h&,b&,x(),dxdt(),ne&,h,wo&,ho&) GRAPHIC SAVE EXE.NAME\$+".bmp" SHELL"bmp2eps "+EXE.NAME\$ WAITKEY\$ CASE "V" 'Display value of variables CLS PRINT SPACE\$(5); FOR i&=1 TO ne& av\$="x("+TRIM\$(STR\$(i&))+")" PRINT av\$+SPACE\$(13-LEN(av\$)); NEXT i& PRINT WHILE q\$<>\$ESC FOR i&=1 TO ne& PRINT FORMAT\$(x(i&),"* ####.0#####"); NEXT i& PRINT CALL RK5 (x(), dxdt(), ne&, h, x()) q\$=WAITKEY\$ WEND CASE "A" 'Display average value of variables nm&=0: nn&=1e5 REDIM xav(ne&) CLS PRINT"Iterations",; FOR i&=1 TO ne& av\$="" PRINT av\$+SPACE\$(11-LEN(av\$)); NEXT i& PRINT"Size" WHILE INKEY\$<>\$ESC xmin=1e99: xmax=-xmin FOR i&=1 TO nn& CALL RK5 (x(), dxdt(), ne&, h, x()) FOR j&=1 TO ne& xav(j&)=xav(j&)+x(j&) xmin=MIN(xmin,x(j&)) xmax=MAX(xmax,x(j&)) NEXT j& NEXT i& INCR nm& its&&=nm&*nn& PRINT USING\$("##########",its&&); FOR i&=1 TO ne& PRINT USING\$(" ##.#^^^^",xav(i&)/its&&); NEXT i& PRINT USING\$(" ##.#^^^^",xmax-xmin) WEND CASE "M" 'Save time-series data for maximum x(1) nimax&=10000 'Number of data points to save DIM xfit(nimax&-1) CLS filenum%=FREEFILE OPEN EXE.NAME\$+".dat" FOR OUTPUT AS filenum% ni&=0 WHILE INKEY\$<>\$ESC CALL RK5 (x(), dxdt(), ne&, h, x()) IF xold>xolder AND xold>x(1) THEN 'A local maximum was crossed CALL extremum (xolder, xold, x(1), xe, ye) PRINT ni&+1;FORMAT\$(ye,"* ####.0#####") PRINT #filenum%, ye xfit(ni&)=ye INCR ni& q\$=INKEY\$ IF ni&=nimax& THEN q\$=\$ESC END IF xolder=xold xold=x(1) WEND CLOSE filenum% CALL quadfit (xfit(), afit, bfit, cfit, rsq) PRINT PRINT " Xnew = ";FORMAT\$(afit,"##.####"); PRINT "X^2 ";FORMAT\$(bfit,"+ ##.####;- ##.####"); PRINT "X ";FORMAT\$(cfit,"+ ##.####;- ##.####") PRINT " Rsquared ="rsq ERASE xfit WAITKEY\$ CASE "X" 'Make plot of time series x versus time DESKTOP GET SIZE TO ddw&, dh& h&=3*w&/4 CONSOLE SET LOC w&+b&+6,0 CLS GRAPHIC WINDOW EXE.NAME\$, 0, 0, w&+b&, b&/8+h&+b&/2 TO hWin??? GRAPHIC ATTACH hWin???, 0, REDRAW GRAPHIC COLOR %BLACK, %WHITE GRAPHIC SET FOCUS: GRAPHIC REDRAW: CONSOLE SET FOCUS PRINT"Press to save image to "+EXE.NAME\$+".bmp: "; CALL timeseries(w&,h&,b&,x(),dxdt(),ne&,h,zs) GRAPHIC SAVE EXE.NAME\$+".bmp" SHELL"bmp2eps "+EXE.NAME\$ WAITKEY\$ END SELECT q\$="" GRAPHIC WINDOW END LOOP END FUNCTION SUB DERIVS (x(), dxdt(), n&) 'Returns the time derivatives dxdt(i&) of x(i&) m=1.4##: q=4##: p=0##: mu=4##: k=1##: e=(q-1##)*mu/(q*q): emu=e/mu sqp=SQR(x(4)): sin2p=SIN(2##*mu*x(3)): cos2p=COS(2##*mu*x(3)): moq=1##-m/q dxdt(1)=m*x(2)+moq*sin2p*sqp dxdt(2)=m*(p-x(4))^2/(4##*mu*mu*x(1)*x(1)*x(1)) dxdt(2)=dxdt(2)-m*e*e*x(1)/4##+k/((q-1##)*x(1)*x(1)) dxdt(2)=dxdt(2)+0.5##*moq*(e-(p-x(4))/(mu*x(1)^2))*cos2p*sqp dxdt(3)=moq*(q-2##)/(2##*q)+0.25##*m*(emu-(p-x(4))/(mu*mu*x(1)*x(1))) c=0.5##*(e*x(1)+(p-3##*x(4))/(mu*x(1)))*cos2p dxdt(3)=dxdt(3)+0.5##*moq*(x(2)*sin2p-c)/sqp c=0.5##*(e*x(1)+(p-x(4))/(mu*x(1)))*sin2p dxdt(4)=-2##*mu*moq*(x(2)*cos2p+c)*sqp 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&,10000) 'Bailout after 10000 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(3)<0 THEN xnew(3)=xnew(3)+pi '?phi IS PERIODIC IF xnew(3)>pi THEN xnew(3)=xnew(3)-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 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)2 THEN zmin=MIN(zmin,x(3)): zmax=MAX(zmax,x(3)) wmin=MIN(wmin,x(1+(3 MOD ne&))): wmax=MAX(wmax,x(1+(3 MOD ne&))) vmin=MIN(vmin,x(1+(4 MOD ne&))): vmax=MAX(vmax,x(1+(4 MOD ne&))) 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 cas&<1 THEN GRAPHIC CLEAR GRAPHIC FONT "MS Sans Serif",12,1 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,wo&),ho&) GRAPHIC PRINT xm\$ xm\$="Y" '?V OR Y GRAPHIC TEXT SIZE xm\$ TO xmw, xmh GRAPHIC SET POS(MAX(wo&+b&-5-xmw,wo&),ho&+h&/2-12) GRAPHIC PRINT xm\$ xm\$=TRIM\$(STR\$(ym1&)) GRAPHIC TEXT SIZE xm\$ TO xmw, xmh GRAPHIC SET POS(MAX(wo&+b&-5-xmw,wo&),ho&+h&-20) GRAPHIC PRINT xm\$ xm\$=TRIM\$(STR\$(xm1&)) GRAPHIC TEXT SIZE xm\$ TO xmw, xmh GRAPHIC SET POS(wo&+b&-5,ho&+h&) GRAPHIC PRINT xm\$ GRAPHIC SET POS(wo&+w&/2+b&-8,ho&+h&) GRAPHIC PRINT"X" 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 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 xp>b& AND xp0 AND yp2 THEN zp=255*(x(3)-zmin)/(zmax-zmin) ELSE zp=255 wp=28*cf&*(x(1+(3 MOD ne&))-wmin)/(wmax-wmin) vp=28*cf&*(x(1+(4 MOD ne&))-vmin)/(vmax-vmin) 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+wp,g1+vp,g1) c2&=RGB(g2+wp,g2+vp,g2) c3&=RGB(g3+wp,g3+vp,g3) c4&=RGB(g4+wp,g4+vp,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 xp2 THEN zmin=MIN(zmin,x(3)): zmax=MAX(zmax,x(3)) wmin=MIN(wmin,dxdt(1)): wmax=MAX(wmax,dxdt(1)) vmin=MIN(vmin,dxdt(2)): vmax=MAX(vmax,dxdt(2)) 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 cas&<1 THEN GRAPHIC CLEAR GRAPHIC FONT "MS Sans Serif",12,1 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,wo&),ho&) GRAPHIC PRINT xm\$ xm\$="Y" '?V OR Y GRAPHIC TEXT SIZE xm\$ TO xmw, xmh GRAPHIC SET POS(MAX(wo&+b&-5-xmw,wo&),ho&+h&/2-12) GRAPHIC PRINT xm\$ xm\$=TRIM\$(STR\$(ym1&)) GRAPHIC TEXT SIZE xm\$ TO xmw, xmh GRAPHIC SET POS(MAX(wo&+b&-5-xmw,wo&),ho&+h&-20) GRAPHIC PRINT xm\$ xm\$=TRIM\$(STR\$(xm1&)) GRAPHIC TEXT SIZE xm\$ TO xmw, xmh GRAPHIC SET POS(wo&+b&-5,ho&+h&) GRAPHIC PRINT xm\$ GRAPHIC SET POS(wo&+w&/2+b&-8,ho&+h&) GRAPHIC PRINT"X" 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) htry=h/10 'A smaller step helps resolve the fine structure (but slower) WHILE INKEY\$<>\$ESC FOR n&=1 TO np& xold=x(1) yold=x(2) zold=x(3) CALL rk5(x(),dxdt(),ne&,htry,x()) 'IF x(3)>zs AND zoldb& AND xp0 AND yp2 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=1/v '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)) 'xm1&=-2: xm2&=2: ym1&=-7: ym2&=7 '?OVERRIDE AUTOSCALING xm&=xm2&-xm1&: ym&=ym2&-ym1& 'xm&=MAX(xm&,ym&): ym&=MAX(xm&,ym&) '?FORCE SCALES TO BE THE SAME GRAPHIC FONT "MS Sans Serif",12,1 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\$ GRAPHIC SET POS(wo&+b&-20,ho&+h&/2-12) GRAPHIC PRINT"Y" '?V OR Y? 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&-5,ho&+h&) GRAPHIC PRINT xm\$ GRAPHIC SET POS(wo&+w&/2+b&-8,ho&+h&) GRAPHIC PRINT"X" 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) htry=h/2 'A smaller step helps resolve the fine structure (but slower) FOR vo=0.01## TO 4.4## STEP 0.2## x(1)=0##: x(2)=vo: x(3)=0## nmax&=5*w&*h& FOR n&=1 TO nmax& xold=x(1) yold=x(2) zold=x(3) CALL rk5(x(),dxdt(),ne&,htry,x()) IF n&100*xm& OR ABS(x(2))>100*ym& THEN EXIT FOR 'unbounded orbit 'IF x(3)>zs AND zoldb& AND xp0 AND yp\$ESC CALL rk5(x(),dxdt(),ne&,htry,x()) d2=0##: e&=0 'Calculate elegance FOR i&=1 TO ne& xo(i&)=ROUND(x(i&),d&) 'IF i&=ne& THEN xo(i&)=0## '?FORCE TIME TO BE ZERO d2=d2+(x(i&)-xo(i&))^2 IF xo(i&)=0 THEN ITERATE 'Don't count zeros IF xo(i&)=1 THEN e&=e&+1: ITERATE e&=e&+LEN(TRIM\$(STR\$(xo(i&))))+1 NEXT i& dxdt=1## 'Avoid equilibrium points CALL rk5(xo(),dxdt(),ne&,h,xnew()) FOR i&=1 TO ne&: dxdt=dxdt*ABS(dxdt(i&)): NEXT i& '?ADD -1 FOR NON-AUTONOMOUS SYSTEM IF dxdt=0## THEN ITERATE IF (e&=ebest& AND d22 THEN zmin=MIN(zmin,x(3)): zmax=MAX(zmax,x(3)) NEXT n& GRAPHIC CLEAR GRAPHIC FONT "MS Sans Serif",14,1 SELECT CASE xm& 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 xm\$=TRIM\$(STR\$(xm&)) GRAPHIC TEXT SIZE xm\$ TO xmw, xmh GRAPHIC SET POS(MAX(b&-5-xmw,wo&),0) GRAPHIC PRINT xm\$ xm\$="Y" '?V OR Y GRAPHIC TEXT SIZE xm\$ TO xmw, xmh GRAPHIC SET POS(MAX(wo&+b&-5-xmw,wo&),ho&+h&/2-12) GRAPHIC PRINT xm\$ GRAPHIC SET POS(MAX(b&-12-xmw,wo&),h&-xmh) GRAPHIC PRINT "-"+xm\$ GRAPHIC SET POS(b&-10,h&) GRAPHIC PRINT "0" GRAPHIC SET POS(w&/2+b&-4,h&) GRAPHIC PRINT"t" GRAPHIC SET POS(b&+w&-xmw-30,h&) GRAPHIC PRINT "1e5" FOR i&=1 TO 9 j&=b&+i&*w&/10 GRAPHIC LINE(j&,0)-(j&,h&/36) GRAPHIC LINE(j&,h&)-(j&,h&-h&/36) NEXT i& FOR i&=1 TO ticks&-1 j&=i&*h&/ticks& GRAPHIC LINE(b&,j&)-(b&+w&/36,j&) GRAPHIC LINE(b&+w&-1,j&)-(b&+w&-w&/36-1,j&) NEXT i& GRAPHIC BOX (b&, 0) - (b&+w&, h&+1) htry=h/1## 'A smaller step helps resolve the fine structure (but slower) n&=0 x(1)=0##: x(2)=1##: x(3)=0## WHILE INKEY\$<>\$ESC INCR n& t=n&*htry xold=x(1) yold=x(2) zold=x(3) CALL rk5(x(),dxdt(),ne&,htry,x()) xp=b&+w&*t/tmax 'IF x(3)>zs AND zoldb& AND xp0 AND yptmax THEN EXIT LOOP END IF IF (n& MOD 1e4)=0 THEN GRAPHIC REDRAW: PRINT t; WEND WAITKEY\$ END SUB SUB basin(w&,h&,b&,x(),dxdt(),ne&,h,zs,wo&,ho&) 'Plots basin of attraction at z = zs nmax&=5e3 'Number of iterations at each initial condition os=0.05## 'Amount of overscan xm1=1e6: ym1=xm1: xm2=-xm1: ym2=xm2: zmin=1e6: zmax=-zmin: vav=0## FOR n&=-100*nmax& TO 100*nmax& 'Find range of values and average CALL rk5(x(),dxdt(),ne&,h,x()) IF n&<1 THEN ITERATE xm1=MIN(xm1,x(1)): ym1=MIN(ym1,x(2)) xm2=MAX(xm2,x(1)): ym2=MAX(ym2,x(2)) xm&=xm2&-xm1& IF ne&>2 THEN zmin=MIN(zmin,x(3)): zmax=MAX(zmax,x(3)) vav=vav+x(2)*x(2) NEXT n& vav=vav/(100*nmax&) 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 cas&<1 THEN GRAPHIC CLEAR GRAPHIC FONT "MS Sans Serif",12,1 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,wo&),ho&) GRAPHIC PRINT xm\$ xm\$="Y" '?V OR Y GRAPHIC TEXT SIZE xm\$ TO xmw, xmh GRAPHIC SET POS(MAX(wo&+b&-5-xmw,wo&),ho&+h&/2-12) GRAPHIC PRINT xm\$ xm\$=TRIM\$(STR\$(ym1&)) GRAPHIC TEXT SIZE xm\$ TO xmw, xmh GRAPHIC SET POS(MAX(wo&+b&-5-xmw,wo&),ho&+h&-20) GRAPHIC PRINT xm\$ xm\$=TRIM\$(STR\$(xm1&)) GRAPHIC TEXT SIZE xm\$ TO xmw, xmh GRAPHIC SET POS(wo&+b&-5,ho&+h&) GRAPHIC PRINT xm\$ GRAPHIC SET POS(wo&+w&/2+b&-8,ho&+h&) GRAPHIC PRINT"X" 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) GRAPHIC REDRAW PRINT dv=.08##*vav: dvtot=0##: nvt&=0 ddw&=1 'Box size FOR xb&=ddw& TO w&-ddw& STEP ddw& FOR yb&=ddw& TO h&-ddw& STEP ddw& x(1)=xm1&+xm&*xb&/w& x(2)=ym1&+ym&*yb&/h& x(3)=zs vavt=0## FOR n&=-nmax& TO nmax& xold=x(1) yold=x(2) zold=x(3) CALL rk5(x(),dxdt(),ne&,h,x()) IF ABS(x(2))>1e9 THEN EXIT FOR 'Orbit escaped IF n&<1 THEN ITERATE 'Not on the attractor yet vavt=vavt+x(2)*x(2) IF vavt>2*vav*nmax& THEN EXIT FOR 'Cannot be on the attractor IF nmax&-n&>500 THEN ITERATE 'Don't need to plot too much 'IF x(3)>zs AND zoldb& AND xp0 AND yp10##*dv OR vavt=0## THEN 'Outside the basin IF ddw&=1 THEN GRAPHIC GET PIXEL(wo&+b&+xb&,yo&+yb&) TO b1& IF b1&<>%BLACK THEN GRAPHIC SET PIXEL(wo&+b&+xb&,ho&+yb&),%GREEN ELSE GRAPHIC BOX(wo&+b&+xb&,ho&+yb&)-(wo&+b&+xb&+ddw&,ho&+yb&+ddw&),0,%GREEN,%GREEN END IF COLOR 2 'Set text color to green if outside basin ELSE COLOR 7 'Otherwise set it to white dvtot=dvtot+dvt INCR nvt& dv=dvtot/nvt& 'Average deviation END IF PRINT xb&,yb&,CLNG(dvt/dv),dv COLOR 7 IF INKEY\$=\$ESC THEN EXIT SUB NEXT yb& GRAPHIC REDRAW NEXT xb& PRINT"Done!" END SUB SUB spatiotemporal(w&,h&,b&,x(),dxdt(),ne&,h,wo&,ho&) 'Makes spatiotemporal plot tmax=1e2 '?MAXIMUM TIME sk&=1 '?PLOT EVERY sk& POINTS cf&=0 '?COLORIFY (IF cf& NOT ZERO) ticks&=10 '?NUMBER OF TICKS ON GRAPH GRAPHIC FONT "MS Sans Serif",14,1 xm\$=TRIM\$(STR\$(ne&)) GRAPHIC TEXT SIZE xm\$ TO xmw, xmh GRAPHIC SET POS(wo&+MAX(b&-5-xmw,wo&),ho&) GRAPHIC PRINT xm\$ xm\$="i" '?i OR X GRAPHIC TEXT SIZE xm\$ TO xmw, xmh GRAPHIC SET POS(MAX(wo&+b&-5-xmw,wo&),ho&+h&/2-12) GRAPHIC PRINT xm\$ GRAPHIC SET POS(wo&+MAX(b&-18,wo&),ho&+h&-xmh) GRAPHIC PRINT "1" t=0## GRAPHIC BOX (wo&+b&, ho&) - (wo&+b&+w&, ho&+h&+1) imax&=MAX(tmax/(h*w&),1) htry=tmax/(imax&*w&) WHILE INKEY\$<>\$ESC IF t<1000 THEN xm\$=TRIM\$(STR\$(t)) ELSE xm\$=TRIM\$(USING\$("##^^^",t)) REPLACE "+" WITH "" IN xm\$ REPLACE "E" WITH "e" IN xm\$ END IF GRAPHIC TEXT SIZE xm\$ TO xmw, xmh GRAPHIC SET POS(wo&+b&,ho&+h&+1) GRAPHIC PRINT xm\$ GRAPHIC SET POS(wo&+w&/2+b&-4,ho&+h&+1) GRAPHIC PRINT"t" t=t+tmax IF t<1000 THEN xm\$=TRIM\$(STR\$(t)) ELSE xm\$=TRIM\$(USING\$("##^^^",t)) REPLACE "+" WITH "" IN xm\$ REPLACE "E" WITH "e" IN xm\$ END IF GRAPHIC TEXT SIZE xm\$ TO xmw, xmh GRAPHIC SET POS(wo&+b&+w&-xmw,ho&+h&+1) GRAPHIC PRINT xm\$ xmin=1e99: xmax=-xmin FOR i&=1 TO ne& STEP sk& 'Find range of values xmin=MIN(xmin,x(i&)) xmax=MAX(xmax,x(i&)) NEXT i& FOR xp&=b&+1 TO b&+w&-2 FOR i&=1 TO imax& CALL rk5(x(),dxdt(),ne&,htry,x()) NEXT i& 'xmin=1e99: xmax=-xmin 'FOR i&=1 TO ne& STEP sk& 'Find range of values ' xmin=MIN(xmin,x(i&)) ' xmax=MAX(xmax,x(i&)) 'NEXT i& FOR yp&=1 TO h&-1 j&=INT(ne&-ne&*yp&/h&) j&=sk&*INT(j&/sk&) '?ACCOUNT FOR ALTERNATE DATA dj=ne&-ne&*yp&/h&-j& dj=dj/sk& '?ACCOUNT FOR ALTERNATE DATA jmin&=1+(j& MOD ne&) jmax&=sk&+(jmin& MOD ne&) '?ACCOUNT FOR ALTERNATE DATA IF jmax&>ne& THEN jmax&=jmax&-ne& u=x(jmax&)*dj+x(jmin&)*(1##-dj) 'Interpolate jmind&=jmin&-cf&: IF jmind&<1 THEN jmind&=jmind&+ne& jmaxd&=jmax&-cf&: IF jmaxd&<1 THEN jmaxd&=jmaxd&+ne& ud=x(jmaxd&)*dj+x(jmind&)*(1##-dj) jminu&=jmin&+cf&: IF jminu&>ne& THEN jminu&=jminu&-ne& jmaxu&=jmax&+cf&: IF jmaxu&>ne& THEN jmaxu&=jmaxu&-ne& uu=x(jmaxu&)*dj+x(jminu&)*(1##-dj) red&=255*(ud-xmin)/(xmax-xmin) grn&=255*(u-xmin)/(xmax-xmin) blu&=255*(uu-xmin)/(xmax-xmin) GRAPHIC SET PIXEL(wo&+xp&,ho&+yp&),RGB(red&,grn&,blu&) NEXT yp& GRAPHIC REDRAW NEXT xp\$ WEND FOR i&=1 TO ticks&-1 j&=b&+i&*w&/ticks& GRAPHIC LINE(wo&+j&,ho&)-(wo&+j&,ho&+h&/36) GRAPHIC LINE(wo&+j&,ho&+h&)-(wo&+j&,ho&+h&-h&/36) j&=i&*h&/ticks& GRAPHIC LINE(wo&+b&,ho&+j&)-(wo&+b&+w&/36,ho&+j&) GRAPHIC LINE(wo&+b&+w&-1,ho&+j&)-(wo&+b&+w&-w&/36-1,ho&+j&) NEXT i& GRAPHIC REDRAW END SUB SUB timeseries(w&,h&,b&,x(),dxdt(),ne&,h,zs) 'Plots time series x(i&) vs. time tmax=1e3 'Maximum time nx&=1 'Number of x values to plot offset=1## 'Vertical offset of each plot (0-1) os=0.05## 'Amount of overscan DIM xo(ne&),ypold(nx&) FOR i&=1 TO ne&: xo(i&)=x(i&): NEXT i& ym1=1e6: ym2=-ym1 FOR n&=1 TO tmax/h 'Find exact range of values CALL rk5(x(),dxdt(),ne&,h,x()) FOR i&=1 TO nx& ym1=MIN(ym1,x(i&)-offset*(nx&-1)) ym2=MAX(ym2,x(i&)) NEXT i& NEXT n& ym1&=INT(ym1-os*(ym2-ym1)) ym2&=CEIL(ym2+os*(ym2-ym1)) ym&=ym2&-ym1& 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 ym1&=(ym1&+ym2&-ym&)/2 ym2&=ym1&+ym& GRAPHIC CLEAR GRAPHIC FONT "MS Sans Serif",14,1 ticks&=10 xm\$=TRIM\$(STR\$(ym2&)) GRAPHIC TEXT SIZE xm\$ TO xmw, xmh GRAPHIC SET POS(wo&+MAX(b&-5-xmw,wo&),ho&) GRAPHIC PRINT xm\$ xm\$="V" GRAPHIC TEXT SIZE xm\$ TO xmw, xmh GRAPHIC SET POS(MAX(wo&+b&-5-xmw,wo&),ho&+h&/2-12) GRAPHIC PRINT xm\$ xm\$=TRIM\$(STR\$(ym1&)) GRAPHIC TEXT SIZE xm\$ TO xmw, xmh GRAPHIC SET POS(wo&+MAX(b&-5-xmw,wo&),ho&+h&-xmh) GRAPHIC PRINT xm\$ t=0## htry=h FOR i&=1 TO ne&: x(i&)=xo(i&): NEXT i& WHILE INKEY\$<>\$ESC GRAPHIC BOX (wo&+b&, ho&) - (wo&+b&+w&, ho&+h&+1),,,%WHITE GRAPHIC BOX (wo&+b&, ho&+h&+1) - (wo&+b&+w&, ho&+h&+b&),,%WHITE,%WHITE FOR i&=1 TO 9 j&=b&+i&*w&/10 GRAPHIC LINE(j&,0)-(j&,w&/36) GRAPHIC LINE(j&,h&)-(j&,h&-w&/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&-1,j&)-(b&+w&-w&/36-1,j&) NEXT i& IF t<1000 OR t>=10*tmax THEN xm\$=TRIM\$(STR\$(t)) ELSE xm\$=TRIM\$(USING\$("##^^^",t)) REPLACE "+" WITH "" IN xm\$ REPLACE "E" WITH "e" IN xm\$ END IF GRAPHIC TEXT SIZE xm\$ TO xmw, xmh GRAPHIC SET POS(wo&+b&,ho&+h&+1) GRAPHIC PRINT xm\$ GRAPHIC SET POS(wo&+w&/2+b&-4,ho&+h&+1) GRAPHIC PRINT"t" IF t<1000 OR t>10*tmax THEN xm\$=TRIM\$(STR\$(t+tmax)) ELSE xm\$=TRIM\$(USING\$("##^^^",t+tmax)) REPLACE "+" WITH "" IN xm\$ REPLACE "E" WITH "e" IN xm\$ END IF GRAPHIC TEXT SIZE xm\$ TO xmw, xmh GRAPHIC SET POS(wo&+b&+w&-xmw,ho&+h&+1) GRAPHIC PRINT xm\$ xmin=1e99: xmax=-xmin FOR i&=1 TO ne& 'Find range of values xmin=MIN(xmin,x(i&)) xmax=MAX(xmax,x(i&)) NEXT i& n&=0 DO INCR n& ti=n&*htry CALL rk5(x(),dxdt(),ne&,htry,x()) xp=b&+w&*ti/tmax FOR i&=1 TO nx& yp=h&*(ym2&-x(i&)+offset*(i&-1))/(ym2&-ym1&) dx=xp-INT(xp) 'Antialias dy=yp-INT(yp) g1=255-255*dx*dy g2=255-255*dx*(1-dy) g3=255-255*(1-dx)*dy g4=255-255*(1-dx)*(1-dy) c1&=RGB(g1,g1,g1) c2&=RGB(g2,g2,g2) c3&=RGB(g3,g3,g3) c4&=RGB(g4,g4,g4) IF xp>b& AND xp0 AND ypb&+1 AND ABS(yp-ypold(i&))>1 THEN GRAPHIC LINE(xpold,ypold(i&))-(xp,yp) xpold=xp: ypold(i&)=yp NEXT i& IF (n& MOD CLNG(tmax/(w&*htry)))=0 THEN GRAPHIC REDRAW IF ti>tmax THEN EXIT LOOP LOOP t=t+tmax IF t>5*tmax THEN EXIT LOOP 'Quit after 5 times through the loop WEND GRAPHIC REDRAW END SUB 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 SUB extremum (yl, yo, yh, xe, ye) 'Given three values yl @ x=-1, yo @ x=0 and yh @ x=1 'Calculate the extremum ye and the value xe at which ye occurs 'Fits the data to a parabola y = yo + bx + cx^2 b = (yh - yl) / 2## c = (yh + yl) / 2## - yo IF c = 0## THEN 'Second derivative is zero (no extremum) xe = 0## ye = yo ELSE xe = -b / (2## * c) ye = yo - b * b / (4## * c) END IF END SUB SUB quadfit (XP(), A, B, C, RSQ) 'Fits a parabola xnew = ax^2 + bx + c to the array xp() D& = 2 'Specify a parabola DIM A(2 * d& + 1), R(D& + 1, D& + 2), T(D& + 2) N& = UBOUND(XP) - LBOUND(XP) 'Have to miss one point A(1) = N& FOR I& = LBOUND(XP) TO UBOUND(XP) - 1 X = XP(I&) Y = XP(I& + 1) FOR J& = 2 TO 2 * D& + 1 A(J&) = A(J&) + X ^ (J&-1) NEXT J& FOR K& = 1 TO D& +1 R(K&, D& + 2) = T(K&) + Y * X ^ (K& - 1) T(K&) = T(K&) + Y * X ^ (K& - 1) NEXT K& T(D& + 2) = T(D& + 2) + Y ^ 2 NEXT I& FOR J& = 1 TO D& + 1 FOR K& = 1 TO D& + 1 R(J&, K&) = A(J& + K& - 1) NEXT K& NEXT J& FOR J& = 1 TO D& + 1 FOR K& = J& TO D& + 1 FOR I& = 1 TO D& + 2 S = R(J&, I&) R(J&, I&) = R(K&, I&) R(K&, I&) = S NEXT I& Z = 1## / R(J&, J&) FOR I& = 1 TO D& + 2 R(J&, I&) = Z * R(J&, I&) NEXT I& FOR K& = 1 TO D& + 1 IF K& = J& THEN ITERATE Z = -R(K&, J&) FOR I& = 1 TO D& + 2 R(K&, I&) = R(K&, I&) + Z * R(J&, I&) NEXT I& NEXT K& NEXT K& NEXT J& P = 0## FOR J& = 2 TO D& + 1 P = P + R(J&, D& + 2) * (T(J&) - A(J&) * T(1) / N&) NEXT J& Q = T(D& + 2) - T(1) ^ 2 / N& Z = Q - P I& = N& - D& - 1 RSQ = P / Q A = R(3, D& + 2) B = R(2, D& + 2) C = R(1, D& + 2) END SUB