'******************************************************************** '* 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

```
```

```
```

```
```

EXPLANATION FILE OF PROGRAM TEQUDIF1 ==================================== The Runge-Kutta Method ---------------------- We present here the Runge-Kutta method of order 4 to integrate an ODE of order 1: Y' = F(X, Y) The development of Y around x coincidates with its Taylor development n of order 4: y = y + h y' + (h^2/2) y" + (h^3/6) y"' + (h^4/24) y"" n+1 n n n n n Other orders may also be used for the Runge-Kutta method, for instance the Euler's method is a Runge-Kutta of order one. Order 4 is often used because it is a good compromise between speed and accuracy. From a current point (x , y ), the next point is determined by: n n K = h f(x , y ) 0 n n k = h f(x + h/2, y + K /2) 1 n n 0 k = h f(x + h/2, y + K /2) 2 n n 1 k = h f(x + h, y + K ) 3 n n 2 and y = y + (1/6)(K + 2 K + K ) n+1 n 0 1 3 It can be proved that the error at each step is about h^5. The only drawback is the time of calculation: at each step, four function evaluations are necessary. From [BIBLI 03]. ------------------------------------------------------------- End of file rkutta.txt.

## No comments:

## Post a Comment