12. Ordinary Differential Equations

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

Example

Now, we consider the following equation.

    dy/dt = 1/(2 - t)2   (y(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

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

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 NMax = 10, RowR = 4, ColR = 1
    Dim N As Long, Info As Long
    Dim Y(NMax) 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.000000000001 '** 1.0e-12
    ATol(0) = RTol(0)
    T = Cells(RowR, ColR)
    Y(0) = Cells(RowR, ColR + 1)
    Info = 0
    For I = 1 To 19
        Tout = Cells(RowR + I, ColR)
        Neval = 0
        Call Derkf(N, AddressOf F, T, Y(), Tout, RTol(), ATol(), Info)
        If Info <> 1 And Info <> 3 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 4(5)-th order Runge-Kutta-Fehlberg adaptive routine which will automatically adjust the step size to obtain a solution satisfying given target accuracy. User 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.

12-02

11 to 12 digits of the result are correct. We see that number of function evaluations increases to satisfy target accuracy 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.

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 NMax = 10, RowR = 4, ColR = 1
    Dim N As Long, Info As Long
    Dim Y(NMax) As Double, T As Double, Tout As Double
    Dim RTol(0) As Double, ATol(0) As Double, I As Integer
    Dim TT As Double, YY(NMax - 1) As Double, YYp(NMax - 1) As Double, IRev As Long
    N = 1
    RTol(0) = 0.000000000001 '** 1.0e-12
    ATol(0) = RTol(0)
    T = Cells(RowR, ColR)
    Y(0) = Cells(RowR, ColR + 1)
    Info = 0
    For I = 1 To 19
        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 Call F(N, TT, YY(), YYp())
        Loop While IRev <> 0
        If Info <> 1 And Info <> 3 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 function, when returned with IRev = 1, users 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.

Top