'Program FIG13.BAS calculates probability distribution function and variance versus time 'by J. C. Sprott 'Compile with the PowerBASIC Console Compiler DEFEXT a-z 'do all calculations in extended (80-bit) precision FUNCTION PBMAIN() n&=3 'number of nonlinear odes h=0.05 'stepsize of integration tmax&=8192 'time to spend at each initial condition s=100## 'scale factor on graph pi=4#*ATN(1##) CLS CONSOLE NAME "Fig 13" CONSOLE SET LOC 808,0 GRAPHIC WINDOW "Fig 13", 0, 0, 800, 640 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, REDRAW CONSOLE SET FOCUS DIM x(n&), dxdt(n&), p&(900), vartot(100) OPEN"fig15.dat" FOR INPUT AS #2 FOR i&=0 TO LOG2(tmax&): INPUT#2,ntot&,t&,vartot(i&): NEXT i& CLOSE RANDOMIZE TIMER OPEN"fig13.dat" FOR APPEND AS #1 WHILE INKEY\$<>CHR\$(27) ' try different random initial conditions x(1)=0.1-0.2*RND: x(2)=0.1-0.2*RND: x(3)=0.1-0.2*RND INCR ntot& OPEN"fig15.dat" FOR OUTPUT AS #2 FOR t&=1 TO tmax&+1 FOR i&=1 TO 1/h CALL rk4(x(),dxdt(),n&,h,x()) xp&=400+400*(x(2)-.707*x(1))/s yp&=320-320*(x(3)-.707*x(1))/s GRAPHIC SET PIXEL(xp&,yp&) NEXT i& IF t&=4000 THEN 'This is for figure 13 FOR i&=1 TO n& PRINT#1,TRIM\$(STR\$(CLNG(x(i&)))) NEXT i& END IF logt&=INT(LOG2(t&)) IF LOG2(t&)=logt& THEN 'This is for figure 15 FOR i&=1 TO n& vartot(logt&)=vartot(logt&)+x(i&)*x(i&) NEXT i& PRINT#2,ntot&,t&,vartot(logt&) END IF NEXT t& CLOSE #2 INCR j& LOCATE 1,1 PRINT j&; GRAPHIC REDRAW GRAPHIC CLEAR WEND CLOSE #1 CLS n&=0: v=0: kur=0 OPEN"fig13.dat" FOR INPUT AS #1 WHILE NOT EOF(1) INPUT#1,x i&=x+400 IF i&<0 THEN i&=0 IF i&>800 THEN i&=800 INCR p&(i&) IF i&>0 AND i&<799 THEN pmax&=MAX(pmax&,p&(i&)) v=v+x*x INCR n& WEND CLOSE sd=SQR(v/n&) OPEN"fig13.dat" FOR INPUT AS #1 WHILE NOT EOF(1) INPUT#1,x kur=kur+(x/sd)^4 WEND CLOSE kur=kur/n&-3 PRINT"N =";n& PRINT"sd =";sd PRINT"kur =";kur FOR i&=1 TO 799 j&=639-639*p&(i&)/pmax& IF i&=1 THEN GRAPHIC SET PIXEL(i&,j&) ELSE GRAPHIC LINE-(i&,j&) NEXT i& GRAPHIC REDRAW GRAPHIC BOX(0,0)-(800,640) GRAPHIC LINE(400,0)-(400,640) GRAPHIC COLOR RGB(255,0,0) FOR i&=1 TO 799 x=i&-400 pdf=n&*EXP(-x*x/(2*sd*sd))/(sd*SQR(2*pi)) j&=639-639*pdf/pmax& IF i&=1 THEN GRAPHIC SET PIXEL(i&,j&) ELSE GRAPHIC LINE-(i&,j&) NEXT i& GRAPHIC REDRAW GRAPHIC SAVE"fig13.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&)=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