'******************************************************************** '* Differential equations of order 1 * '* by Runge-Kutta method of order 4 * '* * '* Basic version by J-P Moreau, Paris * '* (www.jpmoreau.fr) * '* ---------------------------------------------------------------- * '* Reference: "Analyse en Turbo Pascal versions 5.5 et 6.0 By Marc * '* DUCAMP et Alain REVERCHON - Eyrolles, Paris 1991" * '* [BIBLI 03]. * '* ---------------------------------------------------------------- * '* SAMPLE RUN: * '* * '* Example: integrate y'=4x*(y+sqrt(y))/1+x^2 from x=0 to x=1 * '* * '* DIFFERENTIAL EQUATION WITH 1 VARIABLE OF ORDER 1 * '* of type y' = f(x,y) * '* * '* begin value x : ? 0 * '* end value x : ? 1 * '* y value at x0 : ? 1 * '* number of points: ? 11 * '* finesse : ? 10 * '* * '* X Y * '* ------------------------- * '* 0.000000 1.000000 * '* 0.100000 1.040400 * '* 0.200000 1.166400 * '* 0.300000 1.392400 * '* 0.400000 1.742400 * '* 0.500000 2.250000 * '* 0.600000 2.958400 * '* 0.700000 3.920400 * '* 0.800000 5.198400 * '* 0.900000 6.864400 * '* 1.000000 9.000000 * '* * '******************************************************************** defint i-n defdbl a-h,o-z 'ifi:INTEGER Cls print print " DIFFERENTIAL EQUATION WITH 1 VARIABLE OF ORDER 1" print " of type y' = f(x,y)" print print " begin value x : "; : input xi print " end value x : "; : input xf print " y value at x0 : "; : input yi print " number of points: "; : input m print " finesse : "; : input ifi dim t(m+1) 'equadiff1(fp,t,xi,xf,yi,m,ifi) gosub 2000 END 1000 'Example: y'=4x(y+rac(y))/(1+x²) fp=4#*xx*(yy+sqr(yy))/(1#+xx*xx) return '*************************************************************************** '* SOLVING DIFFERENTIAL EQUATIONS WITH 1 VARIABLE OF ORDER 1 * '* of type y' = f(x,y) * '* ----------------------------------------------------------------------- * '* INPUTS: * '* fp Given equation to integrate (see test program) * '* xi, xf Begin, end values of variable x * '* yi Begin Value of y at x=xi * '* m Number of points to calculate * '* fi finesse (number of intermediate points) * '* * '* OUTPUTS: * '* t : real vector storing m results for function y * '*************************************************************************** 2000 'Subroutine Equadiff1 if ifi < 1 then return h = (xf - xi) / ifi / (m-1) y = yi t(1) = yi for i = 1 to m ni = (i - 1) * ifi - 1 for j = 1 to ifi x = xi + h * (ni + j) xx=x : yy=y gosub 1000 'calculate current fp(x,y) a = h * fp xx=x+h/2# : yy=y+a/2# gosub 1000 'calculate fp(x+h/2,y+a/2) b = h * fp yy=y+b/2# gosub 1000 'calculate fp(x+h/2,y+b/2) c = h * fp x = x + h xx=x : yy=y+c gosub 1000 'calculate fp(x,y+c) d = h * fp y = y + (a + b + b + c + c + d) / 6 next j t(i+1) = y next i 'Affiche(t,m,xi,xf) 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 teqdif1.bas