12. Ordinary Differential Equations

Note – The example worksheets are available for download. See 1. How to obtain and use example worksheets.

Example of solving an initial value problem of ordinary differential equations by XLPack is shown below.

Example 1

Now, we consider the following equation.

  
  dy/dt = 1/(2 - t)2   (ｙ(0) = 0.5)

 
The analytical solution is known for this equation as follows.

  
  y = 1/(2 - t)

 
The graph of this function is as follows.

12-01

Solving problem with VBA program (1)

The following is an example code using VBA library subroutine Derkf.

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()
    Dim N As Long, Info As Long
    Dim Y(0) As Double, T As Double, Tout As Double
    Dim RTol(0) As Double, ATol(0) As Double, I As Integer
    N = 1
    RTol(0) = 0.00000001 '** 1.0e-8
    ATol(0) = RTol(0)
    T = Cells(RowR, ColR)
    Y(0) = Cells(RowR, ColR + 1)
    Info = 0
    For I = 1 To Ndata
        Tout = Cells(RowR + I, ColR)
        Neval = 0
        Call Derkf(N, AddressOf F, T, Y(), Tout, RTol(), ATol(), Info)
        If Info <> 1 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

Derkf is the 4(5)-th order Runge-Kutta-Fehlberg adaptive routine which will automatically adjust the step size to obtain a solution satisfying given target accuracy. You should define the differential equation as the external subroutine and call Derkf with specifying the initial value and the target accuracy. Neval is the variable just to count the number of function evaluations.

The following result was obtained by running this program.

The correct results with 8 digits of requested precision is obtained. We see that number of function evaluations increases to satisfy target accuracy immediately after the start and when t increases and a change of y becomes steeper.

Solving problem with VBA program (2)

The following is an example code using the reverse communication version VBA library subroutine Derkf_r.

Sub Start_r()
    Const Ndata = 19, RowR = 5, ColR = 1
    Dim N As Long, Info As Long
    Dim Y(0) As Double, T As Double, Tout As Double
    Dim RTol(0) As Double, ATol(0) As Double, Neval As Long, I As Integer
    Dim TT As Double, YY(0) As Double, YYp(0) As Double, IRev As Long
    N = 1
    RTol(0) = 0.00000001 '** 1.0e-8
    ATol(0) = RTol(0)
    T = Cells(RowR, ColR)
    Y(0) = Cells(RowR, ColR + 1)
    Info = 0
    For I = 1 To Ndata
        Tout = Cells(RowR + I, ColR)
        Neval = 0
        IRev = 0
        Do
            Call Derkf_r(N, T, Y(), Tout, RTol(), ATol(), 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 <> 1 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 defining the object function as the external VBA function, when returned with IRev = 1, you need to compute the object function value at TT and YY, enter it in YYp() and call Derkf_r again. For further details of reverse communication version, see here.

The same result with above is obtained by this program.

Solving problem with VBA program (3)

The following is an example code using the dense output of Derkf. The dense output is the function to obtain a solution by interpolation if possible so as to reduce the number of function evaluations when the output interval is smaller than the step size of adaptive quadrature.

Const NMax = 2, 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_d()
    Dim N As Long, Info As Long
    Dim Y(0) As Double, T As Double, Tout As Double
    Dim RTol(0) As Double, ATol(0) As Double, I As Integer
    Dim Y1(0) As Double, Tend As Double, Mode As Long, Cont As LongPtr
    N = 1
    RTol(0) = 0.00000001 '** 1.0e-8
    ATol(0) = RTol(0)
    T = Cells(RowR, ColR)
    Y(0) = Cells(RowR, ColR + 1)
    Tend = Cells(RowR + Ndata, ColR)
    Mode = 2
    I = 1
    Tout = Cells(RowR + I, ColR)
    Info = 0
    Neval = 0
    Do
        Call Derkf(N, AddressOf F, T, Y(), Tend, RTol(), ATol(), Info, Mode, Cont)
        If Info <> 1 And Info <> 2 Then
            Cells(RowR + I, ColR + 2) = "Error: " + Str(Info)
            Exit Do
        End If
        While T >= Tout
            Call DerkfInt(N, Tout, Y1(), Cont)
            Cells(RowR + I, ColR + 1) = Y1(0)
            Cells(RowR + I, ColR + 2) = Neval
            I = I + 1
            If I > Ndata Then Exit Do
            Tout = Cells(RowR + I, ColR)
            Neval = 0
        Wend
    Loop
End Sub

The following result was obtained by running this program.

There are intervals where the number of function evaluations (Neval) = 0. We can clearly see that the solution is obtained by interpolation in these intervals.

Solving problem with Solver

“ODEs” of XLPack solver add-in can be used to solve this problem. The formulas (=1/(2-B4)^2) is entered in cell B6.

Please refer to here for how to use solver.

Example 2

Next let’s consider the system of ODEs.

We solve the following system of ODEs.

  

  dy1/dt = -2*y1 + y2 - cos(t)
  dy2/dt = 2*y1 - 3*y2 + 3*cos(t) - sin(t)
  (y1 = 1 and y2 = 2 at t = 0)

 

Solving problem with VBA program (4)

The following is an example code using VBA library subroutine Derkf.

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)
    Yp(0) = -2 * Y(0) + Y(1) - Cos(T)
    Yp(1) = 2 * Y(0) - 3 * Y(1) + 3 * Cos(T) - Sin(T)
    Neval = Neval + 1
End Sub

Sub Start_2()
    Dim N As Long, Info As Long
    Dim Y(1) As Double, T As Double, Tout As Double
    Dim RTol(0) As Double, ATol(0) As Double, I As Integer
    N = 2
    RTol(0) = 0.00000001 '**1e-8
    ATol(0) = RTol(0)
    T = Cells(RowR, ColR)
    Y(0) = Cells(RowR, ColR + 1)
    Y(1) = Cells(RowR, ColR + 2)
    Info = 0
    For I = 1 To Ndata
        Tout = Cells(RowR + I, ColR)
        Neval = 0
        Call Derkf(N, AddressOf F2, T, Y(), Tout, RTol(), ATol(), Info)
        If Info <> 1 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) = Neval
    Next
End Sub

The only difference from the case of a single ODE (Example 1) is that multiple formula values are set in the arrays Y() and Yp() (two in this case). Although Y() and Yp() could be simple variables rather than arrays in the case of Example 1, we declared them as arrays for the consistency with VBA’s argument declarations and used only Y(0) and Yp(0).

In order to solve a second-order ODE, it is sometimes converted to the equivalent system of ODEs such that y1 = y and y2 = dy/dt.

The following result was obtained by running this program.

Published: Oct 11, 2015 (Last update: Nov 4, 2022)