'Program DDELEODE.BAS (c) 2006 by J. C. Sprott 'Calculates Lyapunov exponent for the Mackey-Glass equation (Runge-Kutta method) 'Use PowerBasic Console Compiler #DEBUG ERROR OFF DEFEXT a-z FUNCTION PBMAIN() CLS n&=1e2 h=23##/n& 'Iteration step size for tau=23 DIM x(n&),xe(n&),xnew(n&),xenew(n&),dx(n&),dxdt(n&) dr=1e-10 'Separation of orbits FOR i&=0 TO n&: x(i&)=0.9: xe(i&)=x(i&): NEXT i& xe(0)=x(0)+dr 'Nearby initial condition dr2=dr*dr tp&=1e6/n& 'Number of iterations between prints DO CALL RK4 (x(), dxdt(), n&, h, xnew()) CALL RK4 (xe(), dxdt(), n&, h, xenew()) INCR its&& dx2=0 FOR i&=0 TO n& x(i&)=xnew(i&) xe(i&)=xenew(i&) dx=x(i&)-xe(i&) dx2=dx2+dx*dx dx(i&)=dx NEXT i& rs=SQR(dx2/dr2) FOR i&=0 TO n&: xe(i&)=x(i&)+dx(i&)/rs: NEXT i& IF its&&>tp& THEN ltot=ltot+LOG(rs) IF (its&& MOD tp&)=0 THEN le=ltot/(h*(its&&-tp&)) PRINT FORMAT\$(le," #.#######") END IF IF INKEY\$=CHR\$(27) THEN EXIT LOOP END IF LOOP END FUNCTION SUB DERIVS (x(), dxdt(), n&) 'Returns the time derivatives dxdt(i&) of x(i&) dxdt(0)=0.2##*x(n&)/(1##+x(n&)^10)-0.1##*x(0) FOR i&=1 TO n&-1 dxdt(i&)=0.5##*n&*(x(i&-1)-x(i&+1))/23## NEXT i& dxdt(n&)=n&*(x(i&-1)-x(i&))/23## 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& = 0 TO N& XT(I&) = X(I&) + HH * DXDT(I&) NEXT I& CALL DERIVS(XT(), DXT(), N&) FOR I& = 0 TO N& XT(I&) = X(I&) + HH * DXT(I&) NEXT I& CALL DERIVS(XT(), DXM(), N&) FOR I& = 0 TO N& XT(I&) = X(I&) + H * DXM(I&) DXM(I&) = DXT(I&) + DXM(I&) NEXT I& CALL DERIVS(XT(), DXT(), N&) FOR I& = 0 TO N& XNEW(I&) = X(I&) + H6 * (DXDT(I&) + DXT(I&) + 2## * DXM(I&)) NEXT I& ERASE DXM, DXT, XT END SUB