'******************************************************************** '* Differential equations of order N * '* by Runge-Kutta method of order 4 * '* ---------------------------------------------------------------- * '* 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: * '* * '* Example: integrate y"=(2y'y'+yy)/y from x=4 to x=6 * '* with initial conditions: y(4)=2 and y'(4)=-2tan(1) * '* * '* 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 : ? -3.114815 * '* number of points : ? 11 * '* finesse : ? 20 * '* * '* 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 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 " number of points : "; : input m print " finesse : "; : input ifi 'call subroutine equadiffn gosub 2000 END 'Example: y"=(2y'y'+yy)/y 1000 'FUNCTION fp(x,y()) if abs(y(0))<1e-12 then y(0)=1e-12 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 intermediate 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 Affiche gosub 3000 Return 3000 'Subroutine Affiche h=(xf-xi)/(m-1) x=xi-h Cls print print " X Y " print "----------------------" for i=1 to m x=x+h print using " ##.###### ##.######"; x; t(i) next i return 'End of file teqdifn.bas