'Program FIG5.BAS calculates standard deviation and kurtosis 'by J. C. Sprott 'Compile with the PowerBASIC Console Compiler DEFEXT a-z 'do all calculations in extended (80-bit) precision GLOBAL b AS EXT FUNCTION PBMAIN() n&=3 'number of nonlinear odes h=0.01 'stepsize of integration tmax=5e5 'time to spend at each b value CONSOLE NAME "Fig 5" CONSOLE SET LOC 648,0 GRAPHIC WINDOW "Fig 5", 0, 0, 640, 960 TO hWin1& GRAPHIC ATTACH hWin1&, 0 GRAPHIC WIDTH 1 GRAPHIC COLOR RGB(0,0,0), RGB(255,255,255) GRAPHIC CLEAR GRAPHIC ATTACH hWin1&, 0 GRAPHIC BOX(40,0)-(640,440) GRAPHIC BOX(40,480)-(640,920) GRAPHIC STYLE 0 FOR i&=1 TO 10 xp&=40+600*i&/11 GRAPHIC LINE(xp&,0)-(xp&,20) GRAPHIC LINE(xp&,420)-(xp&,440) GRAPHIC LINE(xp&,480)-(xp&,500) GRAPHIC LINE(xp&,900)-(xp&,920) NEXT i& FOR i&=1 TO 9 yp&=440*i&/10 GRAPHIC LINE(40,yp&)-(60,yp&) GRAPHIC LINE(620,yp&)-(640,yp&) NEXT i& FOR i&=1 TO 9 yp&=480+440*i&/10 GRAPHIC LINE(40,yp&)-(60,yp&) GRAPHIC LINE(620,yp&)-(640,yp&) NEXT i& DIM x(n&), dxdt(n&) FOR i&=1 TO n&: x(i&)=RND: NEXT i& bmin=0.11## bmax=0## FOR ib&=1 TO 600 b=bmin+(bmax-bmin)*ib&/600 v=0: k=0 FOR t=0 TO tmax/h CALL rk4(x(),dxdt(),n&,h,x()) d2=0 FOR i&=1 TO n&: d2=d2+x(i&)*x(i&): NEXT i& v=v+d2 var=v/t k=k+d2*d2/var^2 NEXT t sd=SQR(var) kur=k/t-3 PRINT b,sd,kur yp&=440-440*sd/20 IF ib&=1 THEN yp1&=yp& GRAPHIC LINE(39+ib&,yp1&)-(40+ib&,yp&) yp1&=yp& yp&=480-440*kur/2 IF ib&=1 THEN yp2&=yp& GRAPHIC LINE(39+ib&,yp2&)-(40+ib&,yp&) yp2&=yp& NEXT ib& GRAPHIC SAVE"fig5.bmp" BEEP WAITKEY\$ END FUNCTION SUB DERIVS (x(), dxdt(), n&) 'Returns the time derivatives dxdt(i&) of x(i&) FOR i&=1 TO n& i1&=1+i& MOD n& dxdt(i&)=-b*x(i&)+SIN(x(i1&)) NEXT i& 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