Team:St Andrews/project/modelling/models/RK4

Differential Solver (RK4)

Contents 1 Introduction 2 Mathematical basis 3 Implementation 4 References

Introduction

Our method of solving differential equations is based on Fourth Order Runge-Kutta Method. This technique is the most widely used way of numerically solving differential equations and various methods of implementation were looked at. Most of our coding has been based on the work of Aberdeen 2009 iGEM team who used the same method in their modelling. We would like to thank the team for their work and making it available to others like us for future use.

Mathematical basis

While the strict mathematical derivation of the Runge-Kutta method is available( see [http://www.ss.ncu.edu.tw/~lyu/lecture_files_en/lyu_NSSP_Notes/Lyu_NSSP_AppendixC.pdf here for details] ), we omit it here and instead give a brief explanation of the principles behind the technique.

There are several key formulae in the Fourth Order Runge-Kutta algorithm, which are:

The iteration of the x-values is done very simply adding a fixed step-size (h) at each iteration, thus resulting in a constant increase in the x-value according to the value of h chosen for a particular purpose.

The y-value iteration contains a much more elegant set of relationships. It should be clear to you that the y-iteration formula is in fact a weighted average of the four value k values, k₁, k₂, k₃ and k₄. Even at a first glance it should also be apparent that a 'weight' of 2/6 is given to k₂ and k₃ while a smaller weighting of 1/6 is attributed to k₁ and k₂. What then do these k-value correspond to geometrically?

Look first at k₁. It is equal to the function f(x_n,y_n) multiplied by the step-size h, which is just the Euler prediction for Δy. It is the vertical distance from the previous point to the next Euler-predicted point.

Then consider k₂. It should be noted that the x-value at which the function is being evaluated is halfway across the prediction interval and the y-value consists of the current y-value plus half of the Euler predicted Δy. So then the function is evaluated at a point lying halfway between the current point and the Euler-predicted next point. When we evaluate the function at this point what is produced is an estimate of the slope of the solution curve at this halfway point on curve. Multiplying this result by h gives a prediction of the y-jump made by the solution across the entire interval. However, this prediction is made not on the slope of the solution at the left end of the interval but on the estimated slope halfway to the next Euler-predicted point.

A comparison can be made between k₃ and k₂, with the only difference in their formulation being the replacement of k₁ by k₂. Here the function evaluated at this point gives an estimate of the solution slope at the midpoint of the prediction interval but one which is based in the y-jump predicted by k₂ rather than the Euler-predicted one. When multiplied by h we get a further estimate of the y-jump made by the solution across the whole interval.

Finally, k₄ evaluates the function at the right-hand side of the prediction interval. The y-value at which the function is evaluated, y_n+k₃ is an estimate of the y-value at the right-hand side based on the estimate of the y-jump made by k₃. In its entirey, k₄ gives a final estimate of the y-jump made by the solution across the entire width of the prediction interval.

So to provide a summary: each k_i yields an estimate of the y-jump made by the solution across the entire width of the prediction interval h, with each using the previous k_i to make its estimate. The Runge-Kutta formula can be viewed as the y-value of the current point plus a weighted average of four different predictions for the slope.

Implementation

The code which performs the operations in the Runge-Kutta algorithm is shown below.

EXPLANATION OF OPERATION?

The operation of the algorithm is thus; k

void RungeKutta(double y[], double dy[]){
           double k1[var], k2[var], k3[var], k4[var], y2[var], y3[var], y4[var], *pode;
           double h = 0.025;
           pode = ode(y,dy);
           
           //Calculates k1 (for all variables) -> The slope at the start of the interval (h)
           for (int i = 0; i < var; ++i) {
                       k1[i] = *(pode + i) * h;
           }

           //  Calculates k2 (for all variables) -> The slope at the midpoint of the interval (h),
           for (int i = 0; i < var; ++i) {
                       y2[i] =y[i] + k1[i] * 0.5 ;
           }
           pode = ode(y2,dy);
           for (int i = 0; i < var; ++i) {
                       k2[i] = *(pode + i) * h;
           }

           // Calculates k3 (for all variables) -> The slope at the midpoint of... using the y value (y2) determined from k2.
           for (int i = 0; i < var; i++) {
                       y3[i]=y[i] + k2[i] * 0.5;
           }
           pode = ode(y3,dy);
           for (int i = 0; i < var; i++) {
                       k3[i] = *(pode + i) * h;
           }

           // Calculates k4 (for all variables) -> The slope at the end of the interval (h).
           for (int i = 0; i < var; i++) {
                       y4[i]= y[i] + k3[i];
           }
           pode=ode(y4,dy);
           for (int i = 0; i < var; i++) {
                       k4[i] = *(pode + i) * h;
           }

           // Calculates the new y values (for all variables).
           for (int i = 0; i < var; i++){
                       y[i] +=((k1[i] + 2*k2[i] + 2*k3[i] + k4[i])/6);
           }
}

References

[1] Barker.C.A, "Numerical Methods for Solving Differential Equations,The Runge-Kutta Method, Theoretical Introduction", San Joaquin Delta College, 2009 http://calculuslab.deltacollege.edu/ODE/7-C-3/7-C-3-h.html, [Accessed 13 September 2010]

[2] "Numerical Recipes in Fortran 77: The Art of Scientific Computing", Cambridge University Press, 1992, pp 704-707