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