Difficulties with RK4 VBA Function. Function only returns #REF! value

Question

I was wondering if I could get some help finding where I'm going wrong in my function. I've been trying to write a function to use the RK4 method (Runge-Kutta 4) to evaluate an ODE. I've tried changing a few things, but no matter what I do, the function will only return "#REF!". I've attached my RK4 code below (with the corresponding code for the differential equation used in the function). I'm trying to apply Euler method for the first few points, and then transition into using the RK4 method.

 Function RK44(x0 As Double, y0 As Double, n As Integer, xtarg As Double) As Double

 Dim k1 As Double, k2 As Double, k3 As Double, k4 As Double
 Dim ym As Double, ym2 As Double, ye As Double
 Dim Slope As Double, yi As Double, xi As Double, yold As Double

 h = (xtarg - x0) / n
 yi = y0

For i = 1 To n
'apply Euler Method to get y at the end of the interval
yold = yi
k1 = dCdt(yold)
ym = y + k1 * h / 2

k2 = dCdt(ym)
ym2 = y + k2 * h / 2

k3 = dCdt(ym2)
ye = y + k3 * h

k4 = dCdt(ye)

Slope = (1 / 6) * (k1 + 2 * k2 + 2 * k3 + k4)

yi = yold + Slope * h

xi = x0 + i * h
Next i

RK44 = yi
End Function

Here is my function,

Function dCdt(y As Variant)
dCdt = ((0.02) - (0.1) * (y))
End Function

To debug your function, call it from a VBA sub - you can then see where the problem is. Also very useful if you could supply sample arguments for RK44 — Tim Williams, Apr 19 '20 at 05:34

John Alexiou · Answer 1 · 2020-04-19T06:09:44.773

Try this:

Function dCdt(y As Variant) As Variant
    dCdt = ((0.02) - (0.1) * (y))
End Function

Function SIMRK44(ByVal x0 As Double, ByVal y0 As Double, ByVal n As Integer, ByVal xtarg As Double) As Double

    Dim k1 As Double, k2 As Double, k3 As Double, k4 As Double
    Dim ym As Double, ym2 As Double, ye As Double
    Dim Slope As Double, yi As Double, xi As Double, h As Double

    h = (xtarg - x0) / n
    xi = x0
    yi = y0

    For i = 1 To n
        'apply Euler Method to get y at the end of the interval
        k1 = dCdt(yi)
        ym = yi + k1 * h / 2#

        k2 = dCdt(ym)
        ym2 = yi + k2 * h / 2#

        k3 = dCdt(ym2)
        ye = yi + k3 * h

        k4 = dCdt(ye)

        Slope = (1 / 6#) * (k1 + 2# * k2 + 2# * k3 + k4)

        yi = yi + Slope * h
        xi = xi + h
    Next i

    SIMRK44 = yi
End Function

The main problem is the name RK44 refers to a cell so Excel gets confused. Rename the function to SIMRK44 and the problem goes away.

I tweaked the code to be a true simulation where the results of each step depend on the previous step. There was some confusion going on in your code with yold which isn't needed.

Also for better readability each value that is supposed to be a real coefficient and not an integer, I used literals. So 2# instead of 2. Internally the compiler makes the change anyway, it just becomes clearer that 1/6# is not integer division.

Also, VBA defaults function arguments to be passed ByRef which might cause side effects. You must put ByVal before each argument to show Excel that the function will not modify the arguments. Only then it can be used as a UDF.

If I was asked to write a RK4 integrator, then I would stick to close to the standard pseudocode given in textbooks. Simpler and easier to understand what the intent is.

Function f(ByVal x As Double, ByVal y As Double) As Double
    f = 0.02 - 0.1 * y
End Function

Function SIMRK4(ByVal x0 As Double, ByVal y0 As Double, ByVal n As Long, ByVal xtarg As Double) As Double
    Dim xi As Double, yi As Double, h As Double
    Dim k1 As Double, k2 As Double, k3 As Double, k4 As Double
    h = (xtarg - x0) / n

    xi = x0
    yi = y0

    Do
        'Ensure we don't overstep
        If xi + h > xtarg Then
            h = xtarg - xi
        End If
        k1 = f(xi, yi)
        k2 = f(xi + h / 2#, yi + h / 2# * k1)
        k3 = f(xi + h / 2#, yi + h / 2# * k2)
        k4 = f(xi + h, yi + h * k3)

        xi = xi + h
        yi = yi + (h / 6#) * (k1 + 2# * k2 + 2# * k3 + k4)
    Loop Until xi = xtarg
    SIMRK4 = yi
End Function

PS. Also, the code above could handle variable time steps as the number of iterations isn't fixed.

Difficulties with RK4 VBA Function. Function only returns #REF! value

1 Answers1