Here is a presentation of the algorithm and code for the solutions of systems of linear ordinary differential equations by means of the 4th-order Runge-Kutta method, along with some comments.
Contents
- Two reference texts
- The method formulae
- The algorithm for the 4th-order Runge-Kutta method applied to a system of linear ode's
Two reference texts
Much of the background can be found, partially, in well known textbooks. However, neither the code nor the algorithm is usually presented.
Actually, the book of Numerical Methods for Engineers by Chapra, SC and Canale, RP presents an algorithm for this case but have been slightly modified to our purpose. As you may know, an algorith is a general procedure to implement a numerical method so that it can be implemented in any language (that is the problem).
The book of Differential Equations by Zill, DG is also a good source for examples concerning the numerical solution of systems of ode. In other wordsm, you may choose some exercises from Zill's book and test your program.
These two books are recommended for the reader if a wider view is required. There are several other options to improve the 4th-order Runge-Kutta method but I think this a good point of start.
The method formulae
For short, the 4th-order Runge-Kutta method is better presented if a general set of two ode's is used. If more ode are at hand, the formulae can be easily extended.
Then, consider the system of linear ode's,
$\dfrac{dy_1}{dx}=f(x,y_1,y_2)$ Eq. (01)$\dfrac{dy_2}{dx}=g(x,y_1,y_2)$ Eq. (02)
where $f$ and $g$ are two linear functions on $y_1$ and $y_2$. The form of $f$ and $g$ on x is not relevant. Equations (01-02) are subjected to the following conditions,
$y_1(x=x_0)=y_{1i}$ Eq. (03)$y_2(x=x_0)=y_{2i}$ Eq. (04)
Using the 4th-order Runge-Kutta method, the system of Eqs. (01-02) is found to have the solution,
$y_{1,i+1}=y_{1,i} + \dfrac{1}{6}\left[ k_{1,1}+2\left( k_{1,2}+k_{1,3} \right) + k_{1,4} \right]$$y_{2,i+1}=y_{2,i} + \dfrac{1}{6}\left[ k_{2,1}+2\left( k_{2,2}+k_{2,3} \right) + k_{2,4} \right]$ Eq. (05)
where the parameters $k_{1,1}$ through $k_{2,4}$ are calculated from,
$k_{1,1}=f(x_i,y_{1,i},y_{2,i})$
$k_{2,1}=g(x_i,y_{1,i},y_{2,i})$
$k_{1,2}=f\left(x_i+\dfrac{1}{2}h,y_{1,i}+\dfrac{1}{2}k_{1,1},y_{2,i}+\dfrac{1}{2}k_{2,1}\right)$
$k_{2,2}=g\left(x_i+\dfrac{1}{2}h,y_{1,i}+\dfrac{1}{2}k_{1,1},y_{2,i}+\dfrac{1}{2}k_{2,1}\right)$
$k_{1,3}=f\left(x_i+\dfrac{1}{2}h,y_{1,i}+\dfrac{1}{2}k_{1,2},y_{2,i}+\dfrac{1}{2}k_{2,2}\right)$
$k_{2,3}=g\left(x_i+\dfrac{1}{2}h,y_{1,i}+\dfrac{1}{2}k_{1,2},y_{2,i}+\dfrac{1}{2}k_{2,2}\right)$
$k_{1,4}=f(x_i+h,y_{1,i}+k_{1,3}h,y_{2,i}+k_{2,3}h)$
$k_{2,4}=g(x_i+h,y_{1,i}+k_{1,3}h,y_{2,i}+k_{2,3}h)$ Eq. (06)
As you can see, Eqs. (05-06) are an extension of the formulae for the implementation of the method on a single ode.
The algorithm for the 4th-order Runge-Kutta method applied to a system of linear ode's
As mentioned before, this algorithm is based on that presented in the book of Chapra and Canale. This is as follows
Main program
Program
Assign known data as parameters
n = number of ode's
m = number of desired iterations
yi = initial values for the n dependent variables
xi = initial value of the independent variable
xf = desired final value of computation for the independent varaible
dx = (xf - xi)/m <- Step size calculation
Load variables
x = xi <- Initial value of x
Do for i = 1,n
y(i) = yi(i)
End do
Loop for iterations
Do for i = 1,m
xend = x + dx
If (xend > xf) Then
xend = xf
End if
h = dx
Call INTEGRATOR(x,y,n,h,xend)
Print results for every iteration
If (x >= xf) Exit
End do
End program
Subroutine for advancing a step
Subroutine INTEGRATOR(x,y,n,h,xend)
Do
If (xend - x < h) Then
h = xend - x
End if
Call RK4(x,y,n,h)
If (x >= xend) Exit
End do
Return
End subroutine
Subroutine for 4th-order Runge-Kutta applied to a system of ode's
Subroutine RK4(x,y,n,h)
Call DERIVS(x,y,k1)
Do for i = 1,n
ym(i) = y(i) + k1(i)*h/2
End do
Call DERIVS(x + h/2,ym,k2)
Do for i = 1,n
ym(i) = y(i) + k2(i)*h/2
End do
Call DERIVS(x + h/2,ym,k3)
Do for i = 1,n
ye(i) = y(i) + k3(i)*h
End do
Call DERIVS(x + h,ye,k4)
Do for i = 1,n
Slope(i) = (k1(i) + 2*(k2(i) + k3(i)) + k4(i))/6
y(i) = y(i) + Slope(i)*h
End do
x = x + h
Return
End subroutine
Subroutine to include the system of ode's
Subroutine DERIVS(x,y,k)
dydx(1) = ...
dydx(2) = ...
k1(1) = ...
k1(2) = ...
etc.
end subroutine
The code for the 4th-order Runge-Kutta method for systems of oden in F90
At this stage, going from the algorithm should be easy but sometimes it can be challenging. Also, every programmer may arrive to his/her own implementation so that different accuracies, efficiencies and code length can be found.
The implementation presented here is based, as much as possible, on the algorithm already presented. For this, three different f90 files, using modules, were created:
- one for the main program (main.f90)
- another one for the subroutines (MySubroutines.f90)
- and another one for the data types (nrtypes.f90).
How does the program works?
Well, it is not difficult to observe from the algorithm that the main program calls all subroutines via Integrator. Therefore, it is reasonable to plan having all subroutines in a single file: MySubroutines.f90.
Also, depending on what kind of compiler you are using the series of commands to compile and execute the program may differ. In my case, I am using a commercial compiler so that the whole thing is very easy.
The code
You may follow the links below to download the source code files mentioned previously,
The set of ode's solved by this program is the same presented by Chapra and Canale in their book,
$\dfrac{dy_1}{dx}=-0.5y_1$$\dfrac{dy_2}{dx}=4-0.3y_2-0.1y_1$
(which you will find in the DERIVS subroutine already) subject to the conditions,
$y_1(x=0)=4$which you will find in the main.f90 file. The solution of the equations is performed with a step size $h=0.5$ and is presented from $x=0$ to $x=4$.
These program will output the result which is the same presented too in the book of Chapra and Canale.
 |
| Fig. f90 results and those presented by Chapra and Canale in their book for comparison. |
Some improvements to the algorithm
- First, the new system I try to solve have a wider range of integration from $x=-5$ to $x=5$. What I found was that as the program was taking steps it went into the negative values rather than to $x=5$. In order to fix this problem I had to perform some manual calculations to follow the loginc of the program which i did not want.
- What I found was that xout was useless and that it was the cause of the problem. What I did was to erase xout and replace it in the main program loop dor iterations by dx (which is the step size). It works!
No comments:
Post a Comment