'Program WOLFFLOW.BAS calculates the spectrum of LEs for the Lorenz attractor 'Ported from Wolf's Fortran code 'by J. C. Sprott 'Compile with the PowerBASIC Console Compiler DEFEXT a-z 'do all calculations in extended (80-bit) precision FUNCTION PBMAIN() CLS n&=3 'number of nonlinear odes nn&=n&*(n&+1) 'total number of odes (nonlinear + linear) DIM x(nn&),dxdt(nn&),v(nn&),A(nn&),B(nn&),C(nn&) DIM D(nn&),ltot(n&),znorm(n&),gsc(n&) irate&=20 'integration steps per reorthonormalization h=0.01 'stepsize of integration io&=5000 'number of iterations between printouts v(1)=1.0 'initial conditions for nonlinear ODEs v(2)=1.0 'must be within the basin of attraction v(3)=1.0 t=0## FOR i&=n&+1 TO nn& v(i&)=0## 'Don't mess with these; they are problem independent! NEXT i& FOR i&=1 TO n& v((n&+1)*i&)=1## ltot(i&)=0## NEXT i& DO 'This is a built-in 4th order Runge-Kutta integrator: FOR j&=1 TO irate& FOR i&=1 TO nn& x(i&)=v(i&) NEXT i& CALL DERIVS(x(), dxdt(), n&) FOR i&=1 TO nn& A(i&)=dxdt(i&) NEXT i& FOR i&=1 TO nn& x(i&)=v(i&)+(h*A(i&))/2## NEXT i& CALL DERIVS(x(), dxdt(), n&) FOR i&=1 TO nn& B(i&)=dxdt(i&) NEXT i& FOR i&=1 TO nn& x(i&)=v(i&)+(h*B(i&))/2## NEXT i& CALL DERIVS(x(), dxdt(), n&) FOR i&=1 TO nn& C(i&)=dxdt(i&) NEXT i& FOR i&=1 TO nn& x(i&)=v(i&)+h*C(i&) NEXT i& CALL DERIVS(x(), dxdt(), n&) FOR i&=1 TO nn& D(i&)=dxdt(i&) NEXT i& FOR i&=1 TO nn& v(i&)=v(i&)+h*(A(i&)+D(i&)+2##*(B(i&)+C(i&)))/6## NEXT i& t=t+h NEXT j& 'construct new orthonormal basis by gram-schmidt znorm(1)=0## 'normalize the first vector FOR j&=1 TO n& znorm(1)=znorm(1)+v(n&*j&+1)^2 NEXT j& znorm(1)=SQR(znorm(1)) FOR j&=1 TO n& v(n&*j&+1)=v(n&*j&+1)/znorm(1) NEXT j& 'generate a new orthonormal set FOR j&=2 TO n& FOR k&=1 TO j&-1 'make j-1 gsr coefficients gsc(k&)=0## FOR l&=1 TO n& gsc(k&)=gsc(k&)+v(n&*l&+j&)*v(n&*l&+k&) NEXT l& NEXT k& FOR k&=1 TO n& 'construct a new vector FOR l&=1 TO j&-1 v(n&*k&+j&)=v(n&*k&+j&)-gsc(l&)*v(n&*k&+l&) NEXT l& NEXT k& znorm(j&)=0## 'calculate the vector's norm FOR k&=1 TO n& znorm(j&)=znorm(j&)+v(n&*k&+j&)^2 NEXT k& znorm(j&)=SQR(znorm(j&)) FOR k&=1 TO n& 'normalize the new vector v(n&*k&+j&)=v(n&*k&+j&)/znorm(j&) NEXT k& NEXT j& FOR k&=1 TO n& 'update the running vector magnitudes IF znorm(k&)>0 THEN ltot(k&)=ltot(k&)+LOG(znorm(k&)) NEXT k& INCR m& IF (m& MOD io&)=0 THEN 'normalize exponent and print every io& iterations PRINT "Time =";CQUD(t);TAB(17);"LEs ="; lsum=0##: kmax&=0 FOR k&=1 TO n& le=ltot(k&)/t lsum=lsum+le IF lsum>0 THEN lsum0=lsum: kmax&=k& PRINT USING\$("##.######## ", le); NEXT k& IF ltot(1)>0 THEN dky=kmax&-t*lsum0/ltot(kmax&+1) ELSE dky=0 PRINT USING\$(" Dky =##.########", dky) 'Kaplan-Yorke dimension IF INKEY\$=CHR\$(27) THEN EXIT LOOP 'exit when key is pressed END IF LOOP END FUNCTION SUB DERIVS(x(), dxdt(), n&) b=8##/3## sg=10## r=28## ' Nonlinear Lorenz equations: dxdt(1)=sg*(x(2)-x(1)) dxdt(2)=-x(1)*x(3)+r*x(1)-x(2) dxdt(3)=x(1)*x(2)-b*x(3) ' Linearized Lorenz equations: FOR i&=0 TO 2 dxdt(4+i&)=sg*(x(7+i&)-x(4+i&)) dxdt(7+i&)=(r-x(3))*x(4+i&)-x(7+i&)-x(1)*x(10+i&) dxdt(10+i&)=x(2)*x(4+i&)+x(1)*x(7+i&)-b*x(10+i&) NEXT i& END SUB