'Program examine.bas examines the solution of the system in params.dat #DEBUG ERROR OFF DEFEXT a-z GLOBAL a() AS EXT DECLARE SUB DERIVS (x(), dxdt(), n&) DECLARE SUB RK4 (x(), dxdt(), n&, h, xnew()) FUNCTION PBMAIN CONSOLE NAME"Examine" CONSOLE SCREEN 50,80 CLS: CURSOR OFF OPEN"params.dat" FOR INPUT AS #1 INPUT #1, seedo& INPUT #1, n& CLOSE #1 m&=n&*(n&+1) 'Number of parameters h=0.05 'Step size dr=1e-12 'Perturbation for LE calculation dr2=dr*dr DIM x(n&), xnew(n&), dxdt(n&), xe(n&), xenew(n&), dx(n&), a(m&) datafile$=COMMAND$: IF datafile$="" THEN datafile$="params.dat" OPEN datafile$ FOR INPUT AS #1 INPUT #1, seedo& INPUT #1, n& FOR i&=1 TO n&: INPUT #1,x(i&): NEXT i& FOR i&=1 TO m&: INPUT #1,a(i&): NEXT i& CLOSE #1 FOR i&=1 TO n&: xe(i&)=x(i&)+SQR(dr2/n&): NEXT i& ltot=0: lsumtot=0: t&=0: xmin=1 OPEN"biomass.dat" FOR OUTPUT AS #1 OPEN"diverse.dat" FOR OUTPUT AS #2 WHILE INKEY$="" RK4 x(), dxdt(), n&, h, xnew() RK4 xe(), dxdt(), n&, h, xenew() rs2=0 FOR i&=1 TO n& xmin=MIN(x(i&),xmin) 'calculate minimum value of x so far dx(i&)=xenew(i&)-xnew(i&) rs2=rs2+dx(i&)*dx(i&) NEXT i& rs=SQR(rs2/dr2) 'calculate largest Lyapunov exponent and sum of exponents FOR i&=1 TO n&: xe(i&)=xnew(i&)+dx(i&)/rs: NEXT i& INCR t& IF t&>1000 THEN ltot=ltot+LOG(rs) le=ltot/(h*(t&-1000)) FOR i&=1 TO n& sum=0 FOR j&=1 TO n&: sum=sum+a(n&*i&+j&)*x(j&): NEXT j& lsumtot=lsumtot+a(i&)*(1##-a((n&+1)*i&)*x(i&)-sum) NEXT i& lsum=lsumtot/(t&-1000) END IF bm=0 FOR i&=1 TO n& jold&=SCREENX*x(i&) x(i&)=xnew(i&) bm=bm+x(i&) IF i&>SCREENY-2 THEN ITERATE j&=SCREENX*x(i&) IF j&=jold& THEN ITERATE LOCATE i&,1: PRINT SPACE$(SCREENX); LOCATE i&,1: PRINT STRING$(j&, 219); NEXT i& IF (t& MOD 999)=0 THEN LOCATE 50,1: PRINT "LE =";ROUND(le,8); LOCATE 50,21: PRINT CHR$(228);" LE =";ROUND(lsum,8); LOCATE 50,42: PRINT "Xmin =";ROUND(xmin,8); END IF PRINT#1,CSNG(bm/n&) 'save total biomass cm=1 FOR i&=1 TO n& IF n&>1 AND bm>0 THEN cm=cm-ABS(n&*x(i&)/bm-1)/(2*n&-2) END IF NEXT i& PRINT#2,CSNG(cm) 'save biodiversity WEND CLOSE CLS END FUNCTION SUB DERIVS (x(), dxdt(), n&) 'Returns the time derivatives dxdt(i&) of x(i&) FOR i&=1 TO n& asum=0 FOR j&=1 TO n& asum=asum+a(i&*n&+j&)*x(j&) NEXT j& dxdt(i&)=a(i&)*x(i&)*(1##-asum) 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