12. Ordinary Differential Equations

12.1 Overview

Ordinary differential equations (ODEs) frequently appear in scientific and engineering computations. Since they are rarely solvable analytically, numerical methods are often used.

12.1.1 Initial Value Problems

A first-order ordinary differential equation is expressed as follows:
\[
y’ = \frac{dy}{dt} = f(t, y)
\] Here, \(y\) is the unknown function \(y(t)\) dependent on \(t\). The initial value problem of an ordinary differential equation involves determining the function value \(y(t_f)\) at the final point \(t_f\), given the function value at the initial point \(t = t_0\), denoted as \(y = y_0 = y(t_0)\) (known as the initial condition).

12.1.2 First-Order System of Equations

The initial value problem for a first-order ordinary differential equation can be extended to a system of equations by replacing the scalar \(y\) with the vector \(\boldsymbol{y}\):
\[
\boldsymbol{y’} = \boldsymbol{f}(\boldsymbol{y}, t) \\
\boldsymbol{y_0} = \boldsymbol{y}(t_0) \\
\] For example, in the case of two variables, the system can be written in scalar form as follows:
\[
y_1′ = f_1(y_1, y_2, t) \\
y_2′ = f_2(y_1, y_2, t) \\
\] Two initial conditions are required for the scalar equations:
\[
{y_1}_0 = y_1(t_0) \\
{y_2}_0 = y_2(t_0) \\
\]

12.1.3 Higher-Order Ordinary Differential Equations

Higher-order ordinary differential equations, meaning those involving second or higher derivatives, can be rewritten as a system of first-order ordinary differential equations. For example, consider the following differential equation:
\[
y” = f(y, y’, t)
\] By defining \(y_1 = y\) and \(y_2 = y’\), the equation can be transformed into a system of first-order ordinary differential equations:
\[
y_1′ = y_2 \\
y_2′ = f(y_1, y_2, t) \\
\] In this case, two initial conditions are required:
\[
y_1(t_0) = y(t_0) = y_0 \\
y_2(t_0) = y'(t_0) = y’_0 \\
\] Using the same approach, systems of third-order or higher ordinary differential equations can also be solved.

12.2 Solution Methods for Initial Value Problems of Ordinary Differential Equations

12.2.1 Euler’s Method

The variable \(t\) is discretized as \(t_0, t_1, \dots , t_n, \dots\) The intervals between these points are called step sizes and are represented by \(h_n\):
\[
h_n = t_{n+1} – t_n
\] The step size does not necessarily have to be constant, but here we assume it is independent of \(n\) and simply denote it as \(h\). The approximate value of the unknown function \(y(t)\) at \(t_n\) (\(y(t_n)\)) is represented as \(y_n\).

Here, the derivative is approximated using finite differences as follows:
\[
y’_n = f(t_n, y_n) \simeq (y_{n+1} – y_n) / h
\] From this, the following approximation formula is obtained:
\[
y_{n+1} = y_n + hf(t_n, y_n)
\] Using this equation, starting from the initial value \(y_0\), successive approximate solutions \(y_1, y_2, \dots\) at each point are determined. This method is known as Euler’s method.

The Euler’s method formula corresponds to truncating the second and higher order terms in the Taylor series expansion of \(y\), making it a first-order numerical method. The error is of the order of \(h^2\) (\(O(h^2)\)).

12.2.2 Runge-Kutta Method

The Runge-Kutta method is a technique for obtaining accuracy equivalent to a higher-order Taylor expansion without explicitly computing derivatives of \(f()\). By evaluating \(f(t, y)\) at several points \(t\) and \(y\), and forming a weighted average to match the terms up to order \(p\) in the Taylor expansion, a \(p\)-th order formula is constructed.

The general form of the Runge-Kutta method is as follows:
\[
\begin{align}
& y_{n+1} = y_n + h\sum_{i=1}^s b_ik_i \\
& k_i = f(t_n + c_ih, y_n + h\sum_{j=1}^s a_{ij}k_j) \\
\end{align}
\] The parameters \(a_{ij}, b_i, c_i\) are coefficients. Here, we assume \(a_{ij} = 0 (j >= i), c_i = \sum_{j=1}^s a_{ij}\). The value \(s\) is referred to as the stage number.

12.2.2.1 Heun’s Method

For \(s = 2\), the following equations hold to match the terms up to second order in the Taylor expansion:
\[
\begin{align}
& b_1 + b_2 = 1 \\
& b_2c_2 = 1/2 \\
& b_2a_{21} = 1/2 \\
\end{align}
\] Since there are four unknowns and only three equations, there is one degree of freedom. A commonly used approach sets \(c_2 = 1\), leading to the formula known as Heun’s method. It is not typically referred to as a second-order Runge-Kutta method.
\[
\begin{align}
& y_{n+1} = y_n + h(k_1 + k_2)/2 \\
& k_1 = f(t_n, y_n) \\
& k_2 = f(t_n + h, y_n + hk_1) \\
\end{align}
\] This is a second-order formula, with an error term of \(O(h^3)\).

12.2.2.2 Fourth-Order Runge-Kutta Method

For \(s = 4\), a fourth-order formula can be constructed. In this case, there are 12 parameters and 8 equations, leading to 4 degrees of freedom, allowing for various formulations. The most commonly used formula is as follows:
\[
\begin{align}
& y_{n+1} = y_n + h(k_1 + 2k_2 + 2k_3 + k_4)/6 \\
& k_1 = f(t_n, y_n) \\
& k_2 = f(t_n + h/2, y_n + hk_1/2) \\
& k_3 = f(t_n + h/2, y_n + hk_2/2) \\
& k_4 = f(t_n + h, y_n + hk_3) \\
\end{align}
\] This is often referred to simply as the Runge-Kutta method. It is a fourth-order formula, with an error term of \(O(h^5)\).

12.2.3 Embedded Runge-Kutta Methods

The numerical methods presented so far have an error of \(O(h^p)\), meaning that reducing \(h\) improves accuracy. However, decreasing \(h\) also increases the number of \(f()\) evaluations needed to obtain the desired solution, creating a trade-off.

If the accuracy of the obtained result can be estimated, it becomes possible to select an appropriate (necessary and sufficient) step size \(h\) according to the required precision, avoiding unnecessary computation time.

To improve accuracy estimation along with solution computation, formulas have been developed to efficiently estimate errors. The embedded Runge-Kutta method allows for the calculation of two different order approximations within a single formula. The higher-order approximation serves as the computed result, while the other is used for error estimation.

Using this approach, an “adaptive program” can be created to automatically adjust the step size \(h\) based on the estimated error.

As an example, the formula for the embedded Runge-Kutta method known as the 5(4)th order Runge-Kutta-Fehlberg method is shown below. The adaptive program RKF45 using this formula, has been widely used for many years.
\[
\begin{align}
& y_{n+1} = y_n + h((16/135)k_1 + (6656/12825)k_3 + (28561/56430)k_4 – (9/50)k_5 + (2/55)k_6) \space (Fifth-order formula for computation) \\
& y^*_{n+1} = y_n + h((25/216)k_1 + (1408/2565)k_3 + (2197/4104)k_4 – (1/5)k_5) \space (Fourth-order formula for error evaluation) \\
& k_1 = f(t_n, y_n) \\
& k_2 = f(t_n + (1/4)h, y_n + (1/4)hk_1) \\
& k_3 = f(t_n + (3/8)h, y_n + (1/32)h(3k_1 + 9k_2)) \\
& k_4 = f(t_n + (12/13)h, y_n + (1/2197)h(1932k_1 – 7200k_2 + 7296k_3)) \\
& k_5 = f(t_n + h, y_n + h((439/216)k_1 – 8k_2 + (3680/513)k_3 – (845/4104)k_4)) \\
& k_6 = f(t_n + (1/2)h, y_n + h((-8/27)k_1 + 2k_2 – (3544/2565)k_3 + (1859/4104)k_4 – (11/40)k_5)) \\
\end{align}
\] The fifth-order formula is used for computation. The automatic adjustment of \(h\) is performed using the difference from the fourth-order formula for error evaluation, given by \(\delta_{n+1} = ||y_{n+1} – y^*_{n+1}||\).

12.2.4 Comparison of Solution Methods

We will perform calculations using the solution methods described above.

Example 1

Solve the following initial value problem:
\[
\frac{dy}{dt} = 1/(2 – t)^2 \space (y(0) = 0.5)
\] This example has an algebraic solution, given as:
\[
y = 1/(2 – t)
\] The figure below presents a logarithmic scale where the horizontal axis represents error (with greater precision further to the right), and the vertical axis represents computational cost (measured by the number of function calls for \(f(t, y))\). A method is considered more efficient if its graph is located toward the lower right.

Euler’s method is a first-order formula, Heun’s method is second-order, and the Runge-Kutta method is fourth-order. Each of these methods exhibits a linear relationship on a logarithmic plot, with the slope reflecting the respective order.

Derkfa is an adaptive program based on the 5(4)th order Runge-Kutta-Fehlberg method, used in XLPack. Due to the overhead required for automatic step size adjustment, its efficiency is slightly lower than that of a pure fifth-order formula.

12.2.5 Second-Order Ordinary Differential Equations

A special form of second-order ordinary differential equations often appears:
\[
y” = f(y, t) \space (does not depend on y’)
\] Second-order ordinary differential equations can be solved by transforming them into a system of first-order equations, as previously explained. However, for special cases like this, formulas have been developed to solve them directly in their second-order form, such as the Nyström method.

12.3 Solving Ordinary Differential Equations Using XLPack

Derkfa is an adaptive program based on the 5(4)th order Runge-Kutta-Fehlberg method, similar to RKF45, but with an added dense output feature (described later).

Dopn43 is an adaptive program using the 4(3)th order Runge-Kutta-Nyström method. It is designed for solving second-order systems of ordinary differential equations that do not depend on y’.

Both Derkfa and Dopn43 can also be accessed through the XLPack solver.

12.3.1 Single Ordinary Differential Equation

An example of a single (non-system) ordinary differential equation is explained using Example 1.

12.3.1.1 Solving with a VBA Program (1)

Below is a sample program using Derkfa.

Const Ndata = 19, RowR = 5, ColR = 1
Dim Neval As Long

Sub F(N As Long, T As Double, Y() As Double, Yp() As Double)
    Yp(0) = 1 / (2 - T) ^ 2
    Neval = Neval + 1
End Sub

Sub Start()
    Const N = 1
    Dim Y(N - 1) As Double, T As Double, Tout As Double, Tend As Double
    Dim RTol(0) As Double, ATol(0) As Double
    Dim Mode As Long, Info As Long, I As Long
    RTol(0) = 0.00000001 '** 1.0e-8
    ATol(0) = RTol(0)
    Mode = 2
    T = Cells(RowR, ColR)
    Y(0) = Cells(RowR, ColR + 1)
    Tend = Cells(RowR + Ndata, ColR)
    Info = 0
    For I = 1 To Ndata
        Tout = Cells(RowR + I, ColR)
        Neval = 0
        Call Derkfa(N, AddressOf F, T, Y(), Tout, Tend, RTol(), ATol(), Mode, Info)
        If Info <> 0 And Info <> Mode Then
            Cells(RowR + I, ColR + 2) = "Error: " + Str(Info)
            Exit For
        End If
        Cells(RowR + I, ColR + 1) = Y(0)
        Cells(RowR + I, ColR + 2) = Neval
    Next
End Sub

Derkfa automatically adjusts the step size to obtain the solution based on the target accuracy. The differential equation is prepared as a subroutine F, and by specifying the initial value and required accuracy, the solution \(y\) can be computed.

In this example case, since this is not a system of equations, Y() and Yp() could be simple variables instead of arrays. However, for consistency with VBA’s argument, they are declared as arrays and only Y(0) and Yp(0) are used.

Mode = 2 is the standard operation mode, where calculations begin from the initial value T and proceed to the final value Tend. The step size is automatically adjusted to optimize accuracy based on the required precision. If Tout approaches the optimal step size, the step size is adjusted so that the next step aligns precisely with Tout. After this adjustment, the calculation proceeds, and the intermediate results are returned. The computation continues by setting Tout to new value and calling the function again until reaching Tend.

Neval is a variable used solely for counting the number of function evaluations.

Executing this program produces the following results.

The correct result has been obtained with an accuracy of eight decimal places. Additionally, it is observed that the number of function evaluations increases (step size is reduced) both immediately after the start and when \(t\) becomes large and changes rapidly, in order to maintain accuracy.

12.3.1.2 Solving with a VBA Program (2)

Below is an example program using the reverse communication version (RCI) of Derkfa_r.

Sub Start_r()
    Const N = 1
    Dim Y(N - 1) As Double, T As Double, Tout As Double, Tend As Double
    Dim RTol(0) As Double, ATol(0) As Double
    Dim Mode As Long, Info As Long, I As Long
    Dim TT As Double, YY(0) As Double, YYp(0) As Double, IRev As Long
    RTol(0) = 0.00000001 '** 1.0e-8
    ATol(0) = RTol(0)
    Mode = 2
    T = Cells(RowR, ColR)
    Y(0) = Cells(RowR, ColR + 1)
    Tend = Cells(RowR + Ndata, ColR)
    Info = 0
    For I = 1 To Ndata
        Tout = Cells(RowR + I, ColR)
        Neval = 0
        IRev = 0
        Do
            Call Derkfa_r(N, T, Y(), Tout, Tend, RTol(), ATol(), Mode, Info, TT, YY(), YYp(), IRev)
            If IRev = 1 Then
                YYp(0) = 1 / (2 - TT) ^ 2
                Neval = Neval + 1
            End If
        Loop While IRev <> 0
        If Info <> 0 And Info <> Mode Then
            Cells(RowR + I, ColR + 2) = "Error: " + Str(Info)
            Exit For
        End If
        Cells(RowR + I, ColR + 1) = Y(0)
        Cells(RowR + I, ColR + 2) = Neval
    Next
End Sub

Instead of providing the objective function as an external subroutine, when IRev = 1, the function value is calculated using TT and YY and stored in YYp(), after which Derkfa_r is called again.

For more details on RCI, please refer to this link.

Executing this program yields the same results as above.

12.3.1.3 Solving with a VBA Program (3)

Here is an example using the dense output feature of Derkfa. The dense output feature allows the program to return intermediate results without modifying the step size, ensuring that the optimal step size is maintained and reducing the total number of function evaluations. Intermediate solutions are obtained through interpolation.

Note – For more details on dense output, please refer to this link (click on the “Dense Output” tab near the middle of the index).

To use the dense output feature in the program, simply set Mode = 3 instead of Mode = 2. Running the modified version of the program from section 12.3.1.1 with Mode = 3 yields the following results.

There are intervals where the function evaluation count Neval = 0, clearly indicating that interpolation is being used to obtain the solution in those cases.

The dense output feature is particularly useful in situations where high accuracy is not required, but finer intervals for intermediate output are desired – such as for generating smooth graphs.

12.3.1.4 Solving Using a Solver

It is also possible to solve the equation using the XLPack solver add-in’s “ODE” function. The formula (=1/(2-B4)^2) is entered in cell B6.

For more information about the solver, refer to this link.

12.3.2 System of Ordinary Differential Equations

An example of a system of ordinary differential equations is explained below.

Example 2

Solve the following initial value problem:
\[
\frac{d^2y1}{dt^2} = -y1/r \\
\frac{d^2y2}{dt^2} = -y2/r \\
\] where,
\[
r = (y1^2 + y2^2)^{3/2} \\
\] At \(t = 0\),
\[
y_1 = 1, y_2 = 0, dy_1/dt = 0, dy_2/dt = \sqrt{3}
\]

12.3.2.1 Solving with a VBA Program (1)

This example represents a system of two differential equations. However, to solve it using Derkfa, it is treated as a system of four equations as explained in section 12.1.3, where \(y_3 = y_1′, y_4 = y_2’\).

Below is an example program.

Const Ndata = 19, RowR = 5, ColR = 1
Dim Neval As Long

Sub F2(N As Long, T As Double, Y() As Double, Yp() As Double)
    Dim R As Double
    R = (Y(0) ^ 2 + Y(1) ^ 2) ^ (3 / 2)
    Yp(0) = Y(2)
    Yp(1) = Y(3)
    Yp(2) = -Y(0) / R
    Yp(3) = -Y(1) / R
    Neval = Neval + 1
End Sub

Sub Start_2()
    Const N = 4
    Dim Y(N - 1) As Double, T As Double, Tout As Double, Tend As Double
    Dim RTol(0) As Double, ATol(0) As Double
    Dim Mode As Long, Info As Long, I As Long
    RTol(0) = 0.00000001 '**1e-8
    ATol(0) = RTol(0)
    Mode = 2
    T = Cells(RowR, ColR)
    Y(0) = Cells(RowR, ColR + 1)
    Y(1) = Cells(RowR, ColR + 2)
    Y(2) = Cells(RowR, ColR + 3)
    Y(3) = Cells(RowR, ColR + 4)
    Tend = Cells(RowR + Ndata, ColR)
    Info = 0
    For I = 1 To Ndata
        Tout = Cells(RowR + I, ColR)
        Neval = 0
        Call Derkfa(N, AddressOf F2, T, Y(), Tout, Tend, RTol(), ATol(), Mode, Info)
        If Info <> 0 And Info <> Mode Then
            Cells(RowR + I, ColR + 2) = "Error: " + Str(Info)
            Exit For
        End If
        Cells(RowR + I, ColR + 1) = Y(0)
        Cells(RowR + I, ColR + 2) = Y(1)
        Cells(RowR + I, ColR + 3) = Y(2)
        Cells(RowR + I, ColR + 4) = Y(3)
        Cells(RowR + I, ColR + 5) = Neval
    Next
End Sub

The only difference from the single ordinary differential equation case (Example 1) is that multiple values (four in this case) are assigned to the arrays Y() and Yp().

Executing this program produces the following results.

12.3.2.2 Solving with a VBA Program (2)

Since Example 2 does not include first-order terms, Dopn43 can be used, allowing it to be handled as a system of two equations. Below is an example program.

Const Ndata = 19, RowR = 5, ColR = 1
Dim Neval As Long

Sub F3(N As Long, T As Double, Y() As Double, Yp2() As Double)
    Dim R As Double
    R = (Y(0) ^ 2 + Y(1) ^ 2) ^ (3 / 2)
    Yp2(0) = -Y(0) / R
    Yp2(1) = -Y(1) / R
    Neval = Neval + 1
End Sub

Sub Start_3()
    Const N = 2
    Dim T As Double, Y(N - 1) As Double, Yp(N - 1) As Double, Tout As Double, Tend As Double
    Dim RTol(0) As Double, ATol(0) As Double
    Dim Mode As Long, Info As Long, I As Long
    RTol(0) = 0.00000001 '1.0e-8
    ATol(0) = RTol(0)
    Mode = 2
    T = Cells(RowR, ColR)
    Y(0) = Cells(RowR, ColR + 1)
    Y(1) = Cells(RowR, ColR + 2)
    Yp(0) = Cells(RowR, ColR + 3)
    Yp(1) = Cells(RowR, ColR + 4)
    Tend = Cells(RowR + Ndata, ColR)
    Tend = 12
    Info = 0
    For I = 1 To Ndata
        Tout = Cells(RowR + I, ColR)
        Neval = 0
        Call Dopn43(N, AddressOf F3, T, Y(), Yp(), Tout, Tend, RTol(), ATol(), Mode, Info)
        If Info <> 0 And Info <> Mode Then
            Cells(RowR + I, ColR + 2) = "Error: " + Str(Info)
            Exit For
        End If
        Cells(RowR + I, ColR + 1) = Y(0)
        Cells(RowR + I, ColR + 2) = Y(1)
        Cells(RowR + I, ColR + 3) = Yp(0)
        Cells(RowR + I, ColR + 4) = Yp(1)
        Cells(RowR + I, ColR + 5) = Neval
    Next
End Sub

Dopn43 is used in the same way as Derkfa, except for the difference in order.

Executing this program produces the following results:

12.3.2.3 Solving Using a Solver

It is also possible to solve the equation using the XLPack solver add-in’s “2nd order ODE” function.

For more details on the solver, please refer to this link.

Additionally, the downloadable example worksheet includes a case where the solver add-in is used for “ODE”.