'*********************************************************** '* SOLVING DIFFERENTIAL EQUATIONS OF ORDER 1 * '* (Prediction-correction method) * '* ------------------------------------------------------- * '* Reference: * '* "Analyse en Turbo Pascal versions 5.5 et 6.0 By Marc * '* DUCAMP et Alain REVERCHON - Eyrolles, Paris 1991" * '* [BIBLI 03]. * '* * '* Basic version by J-P Moreau, Paris * '* (www.jpmoreau.fr) * '* ------------------------------------------------------- * '* SAMPLE RUN: * '* (Integrate y'=(4xy+4x*sqrt(y))/(1+xý) from x=0 to x=10 * '* with initial condition y(0)=1) * '* * '* Input x begin: 0 * '* Input x end : 10 * '* Input y begin: 1 * '* Input number of points: 10 * '* Input finesse: 10 * '* * '* X Y * '* ----------------------------- * '* 1.000000 8.999984 * '* 2.000000 80.999881 * '* 3.000000 360.999496 * '* 4.000000 1088.998512 * '* 5.000000 2600.996484 * '* 6.000000 5328.992838 * '* 7.000000 9800.986875 * '* 8.000000 16640.977767 * '* 9.000000 26568.964559 * '* 10.000000 40400.946171 * '* ----------------------------- * '* * * '*********************************************************** defint i-n defdbl a-h,o-z MAXDATA=100 dim tx(MAXDATA),ty(MAXDATA) 'integer ifi, n cls print input " Input x begin: ", xi input " Input x end : ", xf input " Input y begin: ", yi input " Input number of points: ", n input " Input finesse: ", ifi gosub 3000 'call equadiff_pc(tx,ty,xi,xf,yi,n,ifi) ' print results F\$="#####.###### #####.######" print print " X Y " print " -----------------------------" for i=1 to n print using F\$; tx(i); ty(i) next i print " -----------------------------" print end 'of main program 1000 'Function y'=(4xy+4x*sqrt(y))/(1+xý) if yy>0 then f = (4*xx*yy+4*xx*sqr(yy))/(1+xx*xx) else f = 0 end if return 2000 'classical Runge-Kutta method of order 4 xx=x : yy=y : gosub 1000 : a=h*f xx=x+h/2 : yy=y+a/2 : gosub 1000 : b=h*f yy=y+b/2 : gosub 1000 : c=h*f x=x+h : xx=x : yy=y+c : gosub 1000 : d=h*f y = y + (a+b+b+c+c+d)/6 return '********************************************************** '* Prediction-correction method * '* ------------------------------------------------------ * '* INPUTS: * '* xi begin x value * '* xf end x value * '* yi begin y value (at x=xi) * '* n number of points to calculate * '* ifi finesse (number of intermediate * '* points (for example 10) * '* OUTPUTS: * '* tx table of x values (m values) * '* ty table of y values (m values) * '* * '* DESCRIPTION: * '* The algorithm has the following steps: * '* 1. calculate y2, y3, y4 using a classical Runge- * '* Kutta method of order 4. * '* 2. from point (x4, y4), first estimate y(n+1) by * '* formula: * '* y(n+1)=y(n) + h/24(55y'(n)-59y'(n-1)+37y'(n-2) * '* -9y'(n-3) * '* then continue with formula: * '* y(n+1)=y(n) + h/24(9y'(n+1)+19y'(n)-5y'(n-1) * '* +y'(n-2), * '* noting that y'(n+1)=f(x(n+1),y(n+1)) with the * '* estimated value of y(n+1), until convergence is * '* obtained. * '********************************************************** 3000 'Subroutine equadiff_pc(tx, ty, xi, xf, yi, n, ifi) dim p(3) z=yi m=n if (m>MAXDATA) or (ifi<1) then return h = (xf-xi)/ifi/m xx=xi : yy=yi : gosub 1000 : p(3)=f tx(0)=xi : ty(0)=yi k=0 for i=1 to m ni = (i-1)*ifi-1 for j=1 to ifi x=xi+h*(ni+j) k=k+1 if (k<4) then 'call runge_kutta(h,x,z) y=z : gosub 2000 : z=y xx=x : yy=z : gosub 1000 : p(3-k)=f else x=x+h w=z+(h/24)*(55*p(0)-59*p(1)+37*p(2)-9*p(3)) 'continue 10 y=w xx=x : yy=y : gosub 1000 w=z+(h/24)*(9*f+19*p(0)-5*p(1)+p(2)) if abs(y-w) > 1e-10 then goto 10 z=w : p(3)=p(2) : p(2)=p(1) p(1)=p(0) : xx=x : yy= z : gosub 1000 : p(0)=f end if next j 'j loop tx(i)=x : ty(i)=z next i 'i loop return ' end of file teqdifpc.bas