Solving a differential equation

• Download demo project - 10.7 Kb
• Download source - 7.0 Kb

Introduction

This article describes how to numerically solve a simple ordinary differential equation with an initial condition. First, a solution of the first order equation is found with the help of the fourth-order Runge-Kutta method. Second, Nyström modification of the Runge-Kutta method is applied to find a solution of the second order differential equation. This article concentrates not on the numerical procedures themselves, but on a way in which derivatives can be passed into integration formulas. Three such methods are illustrated here: a method of virtual function, an interface method and a delegate method.

Background

Let us think about an electrical circuit composed from a resistor R and capacitor C. There is a charge q on the capacitor and the current i flows through the circuit.

Using Ohm's and Kirchhoff's laws, we get the formula for current i:

i = -1/(R C) q

Because i = dq/dt, we obtain the first order differential equation:

dq/dt = 1/(R C) q

Knowing the charge q₀ at the instant t₀, we would like to find the time variation of the charge, which is an unknown function q = q(t). Despite the fact that an analytical formula for this function exists, we want rather to get numerical approximations of q = q(t). Among many methods invented to solve such initial problems, fourth-order Runge-Kutta method is very popular. Having a function of the first derivative and an initial condition:

dy/dt = f(t, y), y(t₀) = y₀

the Runge-Kutta method finds the next value y_n+1 from the present value y_n with the help of following equation:

y_n+1 = y_n + h/6 (k₁ + 2 k₂ + 2 k₃ + k₄)

where h is a selected time interval and coefficients k₁ to k₄ are:

k₁ = f(t_n, y_n)
k₂ = f(t_n + h/2, y_n + h/2 k₁)
k₃ = f(t_n + h/2, y_n + h/2 k₂)
k₄ = f(t_n + h, y_n + h k₃)

It is clearly not difficult to write a program which uses this method. Now, let us take another problem. Imagine a wooden ball which was thrown upwards. We want to calculate the time dependency of its height above the initial level.

In the picture above, m is the ball mass, y is its height, v is its vertical speed, G is its weight and D is the drag of the surrounding air. Using Newton's second law of motion, we get the second order differential equation:

d²y/dt² = (-G - D)/m

where the weight G is:

G = m g

and the air drag D can be written as:

D = 1/2 r_a v² c_d A

Here, g is acceleration due the gravity, r_a is the air density, c_d so called the drag coefficient and A is the ball's cross-sectional area. Initial conditions are t₀, y₀ and v₀.

More generally, we have the following function defining a value of second derivative and initial conditions:

d²y/dt² = f(t, y, dy/dt), y(t₀) = y_0, dy/dt(t₀) = (dy/dt)₀

To solve such differential equations numerically, Nyström modification of the fourth-order Runge-Kutta method can be used:

y_n+1 = y_n + h (dy/dt)_n + h²/6 (k₁ + k₂ + k₃)
(dy/dt)_n+1 = (dy/dt)_n + h/6 (k₁ + 2 k₂ + 2 k₃ + k₄)

where coefficients k₁ to k₄ are:

k₁ = f(t_n, y_n, (dy/dt)_n)
k₂ = f(t_n + 1/2 h, y_n + 1/2 h (dy/dt)_n + h²/8 k₁, (dy/dt)_n + h/2 k₁)
k₃ = f(t_n + 1/2 h, y_n + 1/2 h (dy/dt)_n + h²/8 k₂, (dy/dt)_n + 1/2 k₂)
k₄ = f(t_n + h, y_n + h (dy/dt)_n + h²/2 k₃, (dy/dt)_n + h k₃)

Again, this method is easy to program. The only problem which can arise is how to implement the algorithm while keeping it general. In other words, how to enable various equations to be solved while having only one integration class.

Using the code

I start with the fourth-order Runge-Kutta method for a first order differential equation. As we noticed in the previous section, the integration procedure is easy to write. However, because we want to keep the code general, we need a method of passing the first derivative of the unknown function into an integration procedure. In the following paragraphs I will describe three ways in which this can be done.

First derivative as a virtual function

The classic object oriented approach to this problem is to define a base class declaring the shape of a method which calculates the first derivative, and then to derive a class corresponding to a given equation from this base. The new class should rewrite the first derivative method. Here is how such base class can look.

abstract public class FirstDerivativeBase
{
    abstract public double GetValue(double t, double y);
}

The class FirstDerivativeBase is an abstract one. It contains the virtual function GetValue, which simply returns the value of the first derivative of the unknown function for time t and function value y. An object relating to a given differential equation should be derived from FirstDerivativeBase. Here is an example matching the RC circuit from the beginning of previous section.

class CircuitV : FirstDerivativeBase 
{
    double c;
    double r;

    public CircuitV(double c, double r)
    {
       this.c = 1;
       this.r = 1;

       if (c > 0)
          this.c = c;
       if (r > 0)
          this.r = r;
    }

    public override double GetValue(double t, double q)
    {
       double u = q / c;
       double i = -u / r;
       return i;
    }
}

The class CircuitV is derived from FirstDerivativeBase and rewrites its virtual function GetValue. In our case, this function takes the time t and capacitor charge q and, with the help of capacity c and resistance r, calculates electric current i. This is, in fact, the first derivative of the charge. The Runge-Kutta method is implemented in the IntegratorRK4V object.

public class IntegratorRK4V
{
    FirstDerivativeBase derivative;
    double t0, y0, h;
    double d0;

    public IntegratorRK4V(FirstDerivativeBase derivative, 
        double t0, double y0, double h)
    {
       this.derivative = derivative;
       this.t0 = t0;
       this.y0 = y0;
       this.h = h;
       this.d0 = derivative.GetValue(t0, y0);
    }

    public int Step(out double t1, out double y1, out double d1)
    {
       double k1 = d0;
       double k2 = derivative.GetValue(t0 + h / 2, y0 + h / 2 * k1);
       double k3 = derivative.GetValue(t0 + h / 2, y0 + h / 2 * k2);
       double k4 = derivative.GetValue(t0 + h, y0 + h * k3);

       y1 = y0 + h / 6 * (k1 + 2 * k2 + 2 * k3 + k4);
       t1 = t0 + h;
       d1 = derivative.GetValue(t1, y1);

       t0 = t1;
       y0 = y1;
       d0 = d1;

       return 0;
    }
}

The integrator has a pointer to the FirstDerivativeBase object as its member variable. Other member variables are the initial time and function values -- t0 and y0, respectively --and the time interval h. The integrator calculates the next time value t1, function value f1 and first derivative d1 in its Step method. No storage of the calculated values is supported by the integrator. The following code shows how to compute the time variation of the electric current in our RC circuit.

static void Test1()
{
    CircuitV circV = new CircuitV(0.1,500);

    double t0 = 0;
    double tmax = 90;
   double q0 = 25;
    double h = 5;
    double i0 = circV.GetValue(t0, q0);
    IntegratorRK4V integratorRK4V = new IntegratorRK4V(circV,t0,q0,h);

    Console.WriteLine("TEST1, virtual function");
    Console.WriteLine("t,s     Q,C     I,A");
    Console.WriteLine("-------------------");
    Console.WriteLine("{0,3:F0}{1,8:F1}{2,8:F2}",t0,q0,-i0);

    double t, q, i;
    do
    {
       integratorRK4V.Step(out t, out q, out i);
       Console.WriteLine("{0,3:F0}{1,8:F1}{2,8:F2}", t, q, -i);
    } while (t < tmax);
    Console.WriteLine();
}

We start with a declaration of the object of CircuitV type, immediately defining the circuit capacity and resistance in its constructor. t0 and q0 define the initial time and electric charge. h is the time interval and tmax is the maximum time. The integrator is of the IntegratorRK4 type and gets the object calculating the first derivative as its first constructor parameter. Other constructor parameters are initial time, charge and time interval.

The do while cycle provides the calculation of unknown function values. Notice that the calling program should process results. In this case, I simply display the time, charge and current. Because I wanted to show positive values for electric current, I changed the sign of the i variable. Here are results.

time, t (s)	charge, Q (C)	current, I (A)
0	25.0	0.50
5	22.6	0.45
10	20.5	0.41
15	18.5	0.37
20	16.8	0.34
25	15.2	0.30
30	13.7	0.27
35	12.4	0.25
40	11.2	0.22
45	10.2	0.20
50	9.2	0.18
55	8.3	0.17
60	7.5	0.15
65	6.8	0.14
70	6.2	0.12
75	5.6	0.11
80	5.0	0.10
85	4.6	0.09
90	4.1	0.08

First derivative as an interface

There is no major difference between defining the first derivative as a virtual function of a base class or as an interface. Even the syntax is similar. Instead of an abstract base class, we merely have an interface.

interface FirstDerivative
{
    double GetValue(double x, double y);
}

An object relating to a differential equation to be solved inherits the FirstDerivative interface and is obliged to implement it. Notice that the interface name comes before the method name.

class CircuitI : FirstDerivative
{
    double c;
    double r;

    public CircuitI(double c, double r)
    {
       this.c = c;
       this.r = r;
    }

    double FirstDerivative.GetValue(double t, double q)
    {
       double u = q / c;
       double i = -u / r;
       return i;
    }
}

The integrator uses the pointer to an instance of the FirstDerivative interface as its member variable. It calls the GetValue function of the interface to obtain the first derivative value.

class IntegratorRK4I
{
    FirstDerivative firstDerivative;
    double t0, y0, h;
    double d0;

    public IntegratorRK4I(
        FirstDerivative firstDerivative, double t0, double y0, 
        double h, double d0)
    {
       this.t0 = t0;
       this.y0 = y0;
       this.h = h;
       this.d0 = d0;
       this.firstDerivative = firstDerivative;
    }

    public int Step(out double t1, out double y1, out double d1)
    {
       double k1 = d0;
       double k2 = firstDerivative.GetValue(t0 + h / 2, y0 + h / 2 * k1);
       double k3 = firstDerivative.GetValue(t0 + h / 2, y0 + h / 2 * k2);
       double k4 = firstDerivative.GetValue(t0 + h, y0 + h * k3);

       y1 = y0 + h / 6 * (k1 + 2 * k2 + 2 * k3 + k4);
       t1 = t0 + h;
       d1 = firstDerivative.GetValue(t1, y1);

       t0 = t1;
       y0 = y1;
       d0 = d1;

       return 0;
    }
}

And here is how to put all this gear together.

static void Test2()
{
    CircuitI circI = new CircuitI(0.1, 500);
    FirstDerivative firstDerivative = (FirstDerivative)circI;

    double t0 = 0;
    double tmax = 90;
    double q0 = 25;
    double h = 5;
    double i0 = firstDerivative.GetValue(t0, q0);
    IntegratorRK4I integratorRK4I = new IntegratorRK4I(
        firstDerivative, t0, q0, h, i0);

    Console.WriteLine("TEST2, interface");
    Console.WriteLine("t,s     Q,C     I,A");
    Console.WriteLine("-------------------");
    Console.WriteLine("{0,3:F0}{1,8:F1}{2,8:F2}", t0, q0, -i0);

    double t, q, i;
    do
    {
       integratorRK4I.Step(out t, out q, out i);
       Console.WriteLine("{0,3:F0}{1,8:F1}{2,8:F2}", t, q, -i);
    } while (t < tmax);
    Console.WriteLine();
}

Notice how the interface is passed to the integrator constructor as its first parameter:

FirstDerivative firstDerivative = (FirstDerivative)circI;
IntegratorRK4I integratorRK4I = new IntegratorRK4I(firstDerivative, t0, 
    q0, h, i0);

Results should of course be the same as in the previous paragraph when we used the virtual function of the base object.

First derivative as a delegate

This approach may remind of how function pointers were passed into procedures in the C language. However, there is a significant difference here. In C, all functions are external. C# equivalents of external functions are static methods. We will find that a non-static method can be passed as a delegate object, as well. In our case, it is an integrator who declares a delegate type.

class IntegratorRK4D
{
    public delegate double GetDerivative(double t, double y);
    GetDerivative getDerivative;
    double t0, y0, h;
    double d0;

    public IntegratorRK4D(GetDerivative getDerivative, double t0, 
        double y0, double h, double d0)
    {
       this.getDerivative = getDerivative;
       this.t0 = t0;
       this.y0 = y0;
       this.h = h;
       this.d0 = d0;
    }

    public int Step(out double t1, out double y1, out double d1)
    {
       double k1 = d0;
       double k2 = getDerivative(t0 + h / 2, y0 + h / 2 * k1);
       double k3 = getDerivative(t0 + h / 2, y0 + h / 2 * k2);
       double k4 = getDerivative(t0 + h, y0 + h * k3);

       y1 = y0 + h / 6 * (k1 + 2 * k2 + 2 * k3 + k4);
       t1 = t0 + h;
       d1 = getDerivative(t1, y1);

       t0 = t1;
       y0 = y1;
       d0 = d1;

       return 0;
    }
}

The delegate type is called GetDerivative. It describes a method employing two double arguments and returns the double value of the first derivative. The delegate object is then declared as a member variable and is used by the Step method when the derivative is being integrated. The class CircuitD, which provides the first derivative, knows nothing about the delegate. On the other hand, it must expose a method which matches it.

class CircuitD
{
    double c;
    double r;

    public CircuitD(double c, double r)
    {
       this.c = c;
       this.r = r;
    }

    public double GetDerivative(double t, double q)
    {
       double u = q / c;
       double i = -u / r;
       return i;
    }
}

In the following code, see how the non-static circ.GetDerivative method is passed as a delegate object.

static void Test3()
{
    CircuitD circD = new CircuitD(0.1, 500);

    double t0 = 0;
    double tmax = 90;
    double q0 = 25;
    double h = 5;
    double i0 = circD.GetDerivative(t0, q0);
    IntegratorRK4D integratorRK4D = new IntegratorRK4D(circD.GetDerivative, 
        t0, q0, h, i0);

    Console.WriteLine("TEST3, delegate");
    Console.WriteLine("t,s     Q,C     I,A");
    Console.WriteLine("-------------------");
    Console.WriteLine("{0,3:F0}{1,8:F1}{2,8:F2}", t0, q0, -i0);

    double t, q, i;
    do
    {
       integratorRK4D.Step(out t, out q, out i);
       Console.WriteLine("{0,3:F0}{1,8:F1}{2,8:F2}", t, q, -i);
    } while (t < tmax);
    Console.WriteLine();
}

As in the two previous examples, the numerical results are the same.

Second order differential equation

Finally, I will show a solution for a second order differential equation with the initial condition. The method used will be Nyström modification of the fourth-order Runge-Kutta method. Formally, both the Runge-Kutta method and its Nyström modification are similar. The difference is in the derivative function, where the Nyström method works with the second derivative and therefore must consider the first derivative value as an input. Certainly, coefficients k_i are altered. Let us start with the interface defining the way in which the second derivative value is obtained.

interface SecondDerivative
{
    double GetValue(double x, double y, double dy);
}

The interface is called SecondDerivative and its method GetValue returns the value of the second derivative of the unknown function. There are three input arguments here. x is the independent variable -- usually time, but not necessarily -- y is the dependent variable and dy stands for the value of the first derivative dy/dx. A class relating to the differential equation should implement the interface. Let me consider the example of the wooden ball mentioned in the Background section.

class FlyingBall : SecondDerivative
{
    double dB;    // Ball diameter, m
    double roB;   // Ball density, kg/m^3
    double cD;    // Drag coefficient
    double g;     // Gravitational acceleration, m/s^2
    double roA;   // Air density
    double mB;    // Ball mass, kg
    double aB;    // Ball cross-section area, m^2

    public FlyingBall()
    {
       dB = 0.1;
       roB = 600;
       cD = 0.1;
       g = 9.81;
       roA = 1.29;

       double v = Math.PI * Math.Pow(dB, 3) / 6;
       mB = roB * v;
       aB = 0.25 * Math.PI * dB * dB;
    }

    double SecondDerivative.GetValue(double t, double y, double v)
    {
       double f = -mB * g -Math.Sign(v) * 0.5 * cD * roA * v * v * aB;
       double d2y = f / mB;
       return d2y;
    }
}

Values of parameters like the ball diameter, the material density and so on are directly assigned in the constructor. The drag coefficient cD gets the value 1, which suits well to a smooth sphere. Air density roA (1.29 kg/m³) corresponds with the average value at sea level. The method SecondDerivative.GetValue first calculates the total force acting on the ball. Then it divides the force by the ball mass, thus obtaining the ball acceleration. This acceleration is simply the second derivative of the ball's vertical coordinates. The integrator has a similar structure to those used above.

class IntegratorRKN
{
    SecondDerivative secondDerivative;
    double t0, y0, dy0, h;
    double d2y0;

    public IntegratorRKN(SecondDerivative secondDerivative, double t0, 
        double y0, double dy0, double h, double d2y0)
    {
       this.secondDerivative = secondDerivative;
       this.t0 = t0;
       this.y0 = y0;
       this.dy0 = dy0;
       this.h = h;
       this.d2y0 = d2y0;
    }

    public int Step(out double t1, out double y1, out double dy1, 
        out double d2y1)
    {
       double h2 = h * h; 

       double k1 = d2y0;
       double k2 = secondDerivative.GetValue(
           t0 + h/2, y0 + h/2 * dy0 + h2/8 * k1, dy0 + h/2 * k1);
       double k3 = secondDerivative.GetValue(
           t0 + h/2, y0 + h/2 * dy0 + h2/8 * k2, dy0 + h/2 * k2);
       double k4 = secondDerivative.GetValue(
           t0 + h,   y0 +   h * dy0 + h2/2 * k3, dy0 +   h * k3);

       t1 = t0 + h;
       y1 = y0 + h * dy0 + h2/6 * (k1 + k2 + k3);
       dy1 = dy0 + h/6 * (k1 + 2 * k2 + 2 * k3 + k4);
       d2y1 = secondDerivative.GetValue(t1, y1, dy1);

       t0 = t1;
       y0 = y1;
       dy0 = dy1;
       d2y0 = d2y1;

       return 0;
}

The SecondDerivative interface object is declared as the member variable. Its GetValue method is called by the Step method as a part of the integration process. How is the integrator to be used?

static void Test4()
{
    FlyingBall ball = new FlyingBall();
    SecondDerivative acceleration = (SecondDerivative) ball;

    double t0 = 0;
    double tmax = 10;
    double y0 = 0;
    double h = 0.5;
    double v0 = 50;
    double a0 = acceleration.GetValue(t0, y0, v0);
    IntegratorRKN integratorRKN = new IntegratorRKN(acceleration, t0, y0, v0,
        h, a0);

    Console.WriteLine("TEST4, equation of motion");
    Console.WriteLine(" t,s     y,m   v,m/s  a,m/s2");
    Console.WriteLine("----------------------------");
    Console.WriteLine("{0,4:F1}{1,8:F2}{2,8:F2}{3,8:F2}", t0, y0, v0, a0);

    double t, y, v, a;
    do
    {
       integratorRKN.Step(out t, out y, out v, out a);
       Console.WriteLine("{0,4:F1}{1,8:F2}{2,8:F2}{3,8:F2}", t, y, v, a);
    } while (t < tmax);
    Console.WriteLine();
}

First, we declare an object relating to the differential equation. In this case, it would be FlyingBall. Then we take the interface from the FlyingBall object and store it in the acceleration member. The interface object should be passed as the first parameter of the IntegratorRKN constructor. What remains is to call the Step method repeatedly to get a sequence of vertical positions, velocities and accelerations. Here are results.

time, t (s)	position, y (m)	velocity, v (m/s)	acceleration, a (m/s2)
0.0	0.00	50.00	-13.84
0.5	23.31	43.34	-12.84
1.0	43.41	37.13	-12.03
1.5	60.50	31.28	-11.39
2.0	74.74	25.72	-10.88
2.5	86.26	20.38	-10.48
3.0	95.15	15.22	-10.18
3.5	101.50	10.19	-9.98
4.0	105.35	5.23	-9.85
4.5	106.74	0.32	-9.81
5.0	105.67	-4.58	-9.78
5.5	102.16	-9.45	-9.67
6.0	96.24	-14.24	-9.48
6.5	87.95	-18.92	-9.23
7.0	77.35	-23.46	-8.92
7.5	64.52	-27.83	-8.56
8.0	49.55	-32.01	-8.16
8.5	32.45	-35.98	-7.72
9.0	13.60	-39.73	-7.26
9.5	-7.15	-43.25	-6.79
10.0	-29.61	-46.52	-6.32

The graph below was made in Excel and shows a time variation of the ball height above the initial level.

Points of interest

It is difficult to say which approach to passing a derivative-calculating method to an integrator object is better. The usage of a virtual function is almost equivalent to that of an interface. One can argue that an object relating to a differential equation can possibly inherit more interfaces. While C# does not allow derivation of a class from two or more base classes, the interface approach seems to be more potent. For this reason, I used the interface approach in the last example of the second order differential equation. Usage of the delegate approach is also convenient, mainly due to the fact that a non-static method can serve as a delegate object.

A little bit more questionable is the sequence of actions I used in all of the integrators. The Step function supposes that it knows the present value of the derivative. This value therefore becomes a component of the initial conditions, while at the end of the time step a new value of the derivative is calculated. This can be confusing, because it seems more natural to calculate the derivative value at the beginning of the time interval instead of at the end. The reason I did this is that it allows me to utilize the value of the derivative as an output parameter of the Step method.

Conclusion

This article has described how first and second order differential equations with initial conditions can be solved numerically by the fourth-order Runge-Kutta method and Nyström modification thereof. The article was focused on the manner in which a method calculating the value of the first or second derivative can be passed to an integrator object. It was found that none of the three ways examined provides appreciable advantage over the others.

History

• 24 May, 2007 - Original version posted