# 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.

### 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.

dy_{1}/dt = -2*y_{1}+ y_{2}- cos(t) dy_{2}/dt = 2*y_{1}- 3*y_{2}+ 3*cos(t) - sin(t) (y_{1}= 1 and y_{2}= 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 y_{1} = y and y_{2} = dy/dt.

The following result was obtained by running this program.