'******************************************************************** '* Differential equations of order N * '* by Runge-Kutta method of order 4 * '* with boundary conditions * '* ---------------------------------------------------------------- * '* Reference: "Analyse en Turbo Pascal versions 5.5 et 6.0 by Marc * '* DUCAMP et Alain REVERCHON - Eyrolles, Paris 1991" * '* * '* Basic Version By J-P Moreau, Paris * '* (www.jpmoreau.fr) * '* ---------------------------------------------------------------- * '* SAMPLE RUN: * '* * '* Example: integrate y"=(2y'y'+yy)/y from x=4 to x=6 * '* with initial conditions: y(4)=2 and y'(4)=-2 * '* and final condition y(6)=2 * '* * '* Exact solution is: y = 2cos(1)/cos(x-5) * '* * '* DIFFERENTIAL EQUATION WITH 1 VARIABLE OF ORDER N * '* of type y(n) = f(y(n-1),...,y',y,x) * '* * '* order of equation: ? 2 * '* begin value x : ? 4 * '* end value x : ? 6 * '* y0 value at x0 : ? 2 * '* y1 value at x0 : ? -2 * '* Desired value at xf: ? 2 * '* number of points : ? 11 * '* finesse : ? 20 * '* Max. iterations : ? 200 * '* * '* Last iteration: * '* u =-3.114815521240235 du =-7.62939453125D-07 * '* dyb= 8.19541996488482D-08 67 steps * '* * '* X Y * '* ------------------------ * '* 4.000000 2.000000 * '* 4.200000 1.551018 * '* 4.400000 1.309291 * '* 4.600000 1.173217 * '* 4.800000 1.102583 * '* 5.000000 1.080605 * '* 5.200000 1.102583 * '* 5.400000 1.173217 * '* 5.600000 1.309291 * '* 5.800000 1.551018 * '* 6.000000 2.000000 * '* * '******************************************************************** DEFINT I-N DEFDBL A-H, O-Z 'ifi,i,ndata,iordre: INTEGER DIM yi(10), t(50) DIM ta(10), tb(10), tc(10), td(10), y(10), z(10) CLS PRINT PRINT " DIFFERENTIAL EQUATION WITH 1 VARIABLE OF ORDER N" PRINT " of type y(n) = f(y(n-1),...,y',y,x)" PRINT " with boundary conditions" PRINT PRINT " order of equation: "; : INPUT iordre PRINT PRINT " begin value x : "; : INPUT xi PRINT " end value x : "; : INPUT xf FOR i = 0 TO iordre - 1 PRINT " y"; i; " value at x0 : "; : INPUT yi(i) NEXT i PRINT " Desired value at xf: "; : INPUT yb PRINT " number of points : "; : INPUT m PRINT " finesse : "; : INPUT ifi PRINT " Max. iterations : "; : INPUT itermax niter = 1 u = yi(1) 'First derivative initial value du = ABS(u) / 10# ' estimate initial error on desired final value ' call subroutine equadiffn with no printing iflag = 0: GOSUB 2000 olddyb = ABS(t(m) - yb) dyb = olddyb ' Iterate until desired final value yb is obtained WHILE dyb > .0000001# AND niter < itermax u = u + du: yi(1) = u iflag = 0: GOSUB 2000 dyb = ABS(t(m) - yb) IF ABS(dyb) > ABS(olddyb) THEN 'overshoot? IF niter = 1 THEN du = -du ELSE du = -du / 2# END IF END IF olddyb = dyb niter = niter + 1 WEND PRINT PRINT " Last iteration:" PRINT " u ="; u; " du ="; du PRINT " dyb="; dyb; " "; niter; " steps" ' call subroutine equadiffn and print final results iflag = 1: GOSUB 2000 END 'Example: y"=(2y'y'+yy)/y 1000 'FUNCTION fp(x,y()) IF ABS(y(0)) < .0000000001# THEN y(0) = .0000000001# fp = (2# * y(1) * y(1) + y(0) * y(0)) / y(0) RETURN '*************************************************************************** '* SOLVING DIFFERENTIAL EQUATIONS WITH 1 VARIABLE OF ORDER N * '* of type y(n) = f(y(n-1),y(n-2),...,y',y,x) * '* ----------------------------------------------------------------------- * '* INPUTS: * '* fp Equation to integrate * '* xi, xf Begin, end values of variable x * '* Yi Begin values at xi (of f and derivatives) * '* m Number of points to calculate * '* n Order of differential equation * '* fi finesse (number of intermediary points) * '* * '* OUTPUTS: * '* t real vector storing m results for function y * '* ----------------------------------------------------------------------- * '* EXAMPLE: y" = (2 y'y' + yy) / y with y(4) = 2, y'(4) = -2*tan(1) * '* Exact solution: y = (2 cos(1))/cos(x-5) * '*************************************************************************** 2000 'Subroutine Equadiffn 'h,x,a,b,c,d : double 'ta,tb,tc,td,y,z : Table; 'i,j,k,ni,n1,n2 : integer n = iordre IF ifi < 1 THEN RETURN h = (xf - xi) / ifi / (m - 1) n1 = n - 1: n2 = n - 2 t(1) = yi(0) FOR k = 0 TO n1 y(k) = yi(k): z(k) = yi(k) NEXT k FOR i = 1 TO m ni = (i - 1) * ifi - 1 FOR j = 1 TO ifi x = xi + h * (ni + j) FOR k = 0 TO n1 y(k) = z(k) NEXT k GOSUB 1000 a = h * fp FOR k = 0 TO n2 ta(k) = h * y(k + 1): y(k) = z(k) + ta(k) / 2# NEXT k y(n1) = z(n1) + a / 2# x = x + h / 2 GOSUB 1000 b = h * fp FOR k = 0 TO n2 tb(k) = h * y(k + 1): y(k) = z(k) + tb(k) / 2# NEXT k y(n1) = z(n1) + b / 2# GOSUB 1000 c = h * fp FOR k = 0 TO n2 tc(k) = h * y(k + 1): y(k) = z(k) + tc(k) NEXT k y(n1) = z(n1) + c x = x + h / 2 GOSUB 1000 d = h * fp FOR k = 0 TO n2 z(k) = z(k) + (ta(k) + 2# * tb(k) + 2# * tc(k) + h * y(k + 1)) / 6# NEXT k z(n1) = z(n1) + (a + b + b + c + c + d) / 6# NEXT j t(i + 1) = z(0) NEXT i 'call subroutine Display if required IF iflag <> 0 THEN GOSUB 3000 RETURN 3000 'Subroutine Display h = (xf - xi) / (m - 1) x = xi - h PRINT PRINT " X Y " PRINT " ----------------------" FOR i = 1 TO m x = x + h PRINT USING " ##.###### ##.######"; x; t(i) NEXT i PRINT " ----------------------" RETURN 'End of file tshoot.bas