Programming RK4, need to print values but within a loop, choosing cell a5 first, and then next cycle it will print to a6, etc

Question

I was wondering whether someone would be able to help me out with a problem I am having with a printing function in excel VBA. I am trying to print out the outputs of an ODE using RK4, and I would like to know how to loop a printvalue function, so that every loop an increase in the cell value occurs.

Also wondering why the line PrintValues (x.ToString()) is saying that the x is an invalid qualifier? If there are any other mistakes you feel are in the code, please let me know, although I am fairly certain everything works fine Any help would be greatly appreciated :)

Option Explicit

Public Function RK_KCTest() 'ByVal rateConst As Double, ByVal equilConst As Double, ByVal initConc As Double, ByVal epsilon As Double, ByVal initFlow As Double, ByVal dv As Double, ByVal xi As Double, ByVal vi As Double, ByVal vf As Double, ByVal vout As Double) As Double
'NOTE: This is the ultimate function, that will essentially take all the function's calculations below it, and accumate them into a matrix which it will then print to the user, displaying the results

    Dim rateConst As Double
    rateConst = GetValue("F9")
    Dim equilConst As Double
    equilConst = GetValue("F7")
    Dim initConc As Double
    initConc = GetValue("F5")
    Dim epsilon As Double
    epsilon = GetValue("F8")
    Dim initFlow As Double
    initFlow = GetValue("F6")
    Dim dv As Double
    dv = GetValue("F11")
    Dim xi As Double
    xi = GetValue("F12")
    Dim vi As Double
    vi = GetValue("F10")
    Dim vf As Double
    vf = GetValue("F13")
    Dim vout As Double
    vout = GetValue("F15")

    Dim i As Integer
    i = 0
    Dim m As Integer
    m = 0
    Dim x As Double
    x = xi
    Dim v As Double
    v = vi
    Dim vCutoff

    Dim dx As Double
    dx = Derivatives(x, v, rateConst, equilConst, initConc, epsilon, initFlow, dx)

    Dim h As Double
    h = 50
    Dim hi As Double
    hi = h

    Dim xnew As Double
    Dim xerror As Double
    Dim hNext As Double
    Dim xScaled As Double

    Do 'Loop for the calculation of RK4, calling up the subfunction for each iteration, until the difference between the the final value of v and the value before it is less than that of the stepsize
        vCutoff = v + vout
        h = dv 'h which is the step size is equal to the change of v, basically just another term for stepsize'

        If (vCutoff - v < h) Then h = vCutoff - v  'This is basically a trimmming step, this is used if the increment exceeds the solution determining v value (If the program calculates past the point needed)

        If (vCutoff > vf) Then vCutoff = vf

        Call RK_KC(x, v, rateConst, equilConst, initConc, epsilon, initFlow, dx, h, dv, xnew, xerror) 'calls the RK4 values calculated in the subfunction, incorparating them into the loop, and declaring the new value of x, so the next iteration is able to start
        x = xnew

        PrintValue ("x")
        Range("a5").Select
        PrintValues (x.ToString())


        If (v >= vf) Then Exit Do 'If the volume is greater than or equal to the cutoff volume then the program will exit the loop and the final values or solutions can be determined

    Loop

End Function

Public Function RK_KC(ByRef x As Double, ByRef v As Double, ByVal rateConst As Double, ByVal equilConst As Double, ByVal initConc As Double, ByVal epsilon As Double, ByVal initFlow As Double, ByRef dx As Double, ByRef h As Double, ByRef dv As Double, ByVal xnew As Double, ByVal xerror As Double)
    'NOTE: This function is called by the Integrator subfunction, and the calculations in this function are implemented into a loop
    'so therefore iterations for every value of x and vnew are computed, and the subfunctions calling the integrator function will display the results in an excel spreadsheet

    'Implementing the RKF_KC method into the main function by defining the six k values, as well as the temp values of x calculated between each stored value of k2,k3,k4,k5,k6
    Dim k1 As Double, k2 As Double, k3 As Double, k4 As Double, k5 As Double, k6 As Double
    Dim xTemp As Double 'Called xTemp because it is a placeholder

    Dim a2, a3, a4, a5, a6 As Double
    a2 = 0.2
    a3 = 0.3
    a4 = 0.6
    a5 = 1
    a6 = 0.875
    'The parameters of A column in the butcher table for Cash-Karp RKF method

    Dim b21, b31, b32, b41, b42, b43, b51, b52, b53, b54, b61, b62, b63, b64, b65 As Double
    b21 = 0.2
    b31 = 3 / 40
    b32 = 9 / 40
    b41 = 0.3
    b42 = -0.9
    b43 = 1.2
    b51 = -11 / 54
    b52 = 2.5
    b53 = -70 / 27
    b54 = 35 / 27
    b61 = 1631 / 55296
    b62 = 175 / 512
    b63 = 575 / 13824
    b64 = 44275 / 110592
    b65 = 253 / 4096
    'The parameters of B matrix in the butcher table for Cash-Karp RKF method

    Dim c1, c3, c4, c6 As Double
    c1 = 37 / 378
    c3 = 250 / 621
    c4 = 125 / 594
    c6 = 512 / 1771
    'The parameters of C row in the butcher table for Cash-Karp RKF method

    Dim dc1, dc3, dc4, dc5, dc6 As Double 'The parameters of the differences of between the fourth and fifth order Cash-Karp RKF methods
    dc1 = c1 - (2825 / 27648)
    dc3 = c3 - (18575 / 48384)
    dc4 = c4 - (13525 / 55296)
    dc5 = -277 / 14336
    dc6 = c6 - 0.25
    'The parameters of the differences of between the fourth and fifth order Cash-Karp RKF methods, also known as the local truncation error

    xTemp = x + b21 * h * dx 'Calculates the temporary value of x between v;x and v+h;xnew based on k1
    x = xTemp
    v = v + a2 * h
    k2 = Derivatives(x, v, rateConst, equilConst, initConc, epsilon, initFlow, dx)  'Calls the derivative function and calculates the value of k2 using the stepsize, "a" parameters, and the temporary value of x

    xTemp = x + h * (b31 * dx + b32 * k2) 'Calculates the temporary value of x between v;x and v+h;xnew based on k2
    x = xTemp
    v = v + a3 * h
    k3 = Derivatives(x, v, rateConst, equilConst, initConc, epsilon, initFlow, dx) 'Calls the derivative function and calculates the value of k3 using the stepsize, "b" parameters, and the temporary value of x

    xTemp = x + h * (b41 * dx + b42 * k2 + b43 * k3) 'Calculates the temporary value of x between v;x and v+h;xnew based on k3
    x = xTemp
    v = v + a4 * h
    k4 = Derivatives(x, v, rateConst, equilConst, initConc, epsilon, initFlow, dx) 'Calls the derivative function and calculates the value of k4 using the stepsize, "b" parameters, and the temporary value of x

    xTemp = x + h * (b51 * dx + b52 * k2 + b53 * k3 + b54 * k4) 'Calculates the temporary value of x between v;x and v+h;xnew based on k4
    x = xTemp
    v = v + a5 * h
    k5 = Derivatives(x, v, rateConst, equilConst, initConc, epsilon, initFlow, dx) 'Calls the derivative function and calculates the value of k5 using the stepsize, "b" parameters, and the temporary value of x

    xTemp = x + h * (b61 * dx + b62 * k2 + b63 * k3 + b64 * k4 + b65 * k5) 'Calculates the temporary value of x between v;x and v+h;xnew based on k5
    x = xTemp
    v = v + a6 * h
    k6 = Derivatives(x, v, rateConst, equilConst, initConc, epsilon, initFlow, dx) 'Calls the derivative function and calculates the value of k6 using the stepsize, "b" parameters, and the temporary value of x

    xnew = x + h * (c1 * dx + c3 * k3 + c4 * k4 + c6 * k6) 'The new value of x based on the Cash-Karp RKF butcher table

    xerror = h * (dc1 * dx + dc3 * k3 + dc4 * k4 + dc5 * k5 + dc6 * k6) 'the error value of x associated with Cash-Karp RKF

End Function

Public Function Derivatives(ByRef x As Double, ByVal v As Double, ByVal rateConst As Double, ByVal equilConst As Double, ByVal initConc As Double, ByVal epsilon As Double, ByVal initFlow As Double, ByRef dx As Double)
    'NOTE: This function is called by the RK4 subfunction to calculate k values

    'Defined ODE, which is the change of conversion (dx) over the change of volume(dv). Within this ODE there are constant values of
    'Initial flowrate, Initial Concentration, the epsilon value, the equilibrium constant. the reaction rate constant, and finally the value of conversion which is x
    dx = ((rateConst * initConc) / (initFlow)) * ((((1 - x) / (1 + (epsilon * x))) - ((4 * initConc * x ^ (2)) / (equilConst * (1 + (epsilon * x)) ^ (2)))))

End Function


Public Function GetValue(cellToGet As String) As Double
    Dim wb As Workbook
    Dim ws As Worksheet
    Dim TxtRng  As Range
    Dim returnVal As Double

    Set wb = ActiveWorkbook
    Set ws = wb.Sheets("Sheet1")

    ws.Unprotect

    Set TxtRng = ws.Range(cellToGet)

    returnVal = CDbl(val(TxtRng.Value))

    GetValue = returnVal
    Exit Function

End Function

Public Function PrintValue(str As String)
    Dim wb As Workbook
    Dim ws As Worksheet
    Dim TxtRng  As Range

    Set wb = ActiveWorkbook
    Set ws = wb.Sheets("Sheet1")

    ws.Unprotect

    Set TxtRng = ws.Range("A5")
    TxtRng.Value = str

End Function

Answer 1

Here is an example of using Cash-Karp with a different diff eq.

Option Explicit

Const PI As Double = 3.14159265358979

Public Function CalcSlope(ByVal x As Double, ByVal y As Double) As Double
    CalcSlope = Cos(PI * x) / Sqr(1 + y ^ 2)
End Function

Public Function CASHKARP_1(ByVal h As Double, ByVal x0 As Double, ByVal y0 As Double, ByRef y_err) As Double
    Dim K1 As Double, K2 As Double, K3 As Double, K4 As Double, _
    K5 As Double, K6 As Double

    K1 = h * CalcSlope( _
        x0, _
        y0)

    K2 = h * CalcSlope( _
        h * (1# / 5#) + x0, _
        K1 * (1# / 5#) + y0)

    K3 = h * CalcSlope( _
        h * (3# / 10#) + x0, _
        K1 * (3# / 40#) + K2 * (9# / 40#) + y0)

    K4 = h * CalcSlope( _
        h * (3# / 5#) + x0, _
        K1 * (3# / 10#) + K2 * (-9# / 10#) + K3 * (6# / 5#) + y0)

    K5 = h * CalcSlope( _
         h * (1#) + x0, _
         K1 * (-11# / 54#) + K2 * (5# / 2#) + K3 * (-70# / 27#) + K4 * (35# / 27#) + y0)

    K6 = h * CalcSlope( _
         h * (7# / 8#) + x0, _
         K1 * (1631# / 55296#) + K2 * (175# / 512#) + K3 * (575# / 13824#) + K4 * (44275# / 110592#) + K5 * (253# / 4096#) + y0)


    CASHKARP_1 = (37# / 378#) * K1 + (250# / 621#) * K3 _
        + (125# / 594#) * K4 + (512# / 1771#) * K6

    y_err = (2825# / 27648# - 37# / 378#) * K1 + _
            (18575# / 48384# - 250# / 621#) * K3 + _
            (13525# / 55296# - 125 / 594#) * K4 + _
            (277# / 14336) * K5 + (1# / 4# - 512# / 1771#) * K6

End Function

Public Sub SolveRk()
    Dim count As Long

    count = Range(Range("F3"), Range("F3").End(xlDown)).Rows.count
    Range("F3").Resize(count, 3).ClearContents

    Dim x0 As Double, y0 As Double, h0 As Double, xf As Double

    x0 = Range("x0")
    y0 = Range("y0")
    xf = Range("xf")
    h0 = Range("h0")

    Dim x As Double, h As Double, y As Double, e As Double, r As Double
    x = x0: h = h0: y = y0

    Dim rk_results() As Variant
    ReDim rk_results(1 To 10, 1 To 3)

    Const tol As Double = 0.0000001
    count = 1
    rk_results(1, 1) = x0: rk_results(1, 2) = y: rk_results(1, 3) = h
    Do
        x = x0 + h
        y = y0 + CASHKARP_1(h, x0, y0, e)

        If Abs(e) < tol Then
            count = count + 1

            If UBound(rk_results, 1) < count Then
                ExpandMatrix rk_results, 2 * UBound(rk_results, 1)
            End If

            rk_results(count, 1) = x: rk_results(count, 2) = y: rk_results(count, 3) = h

            If x - xf > -tol Then Exit Do
            h = h * 1.15
            If x + h > xf Then
                h = xf - x
            End If
            x0 = x: y0 = y
        Else
            h = h / 1.2
        End If

    Loop

    Range("F3").Resize(count, 3).Value = rk_results

End Sub

Public Sub ExpandMatrix(ByRef values() As Variant, Optional ByVal new_rows As Long = 0, Optional ByVal new_cols As Long = 0)
    If new_rows = 0 Then new_rows = UBound(values, 1)
    If new_cols = 0 Then new_cols = UBound(values, 2)

    Dim temp() As Variant
    ReDim temp(1 To new_rows, 1 To new_cols)

    Dim i As Long, j As Long

    For i = 1 To UBound(values, 1)
        For j = 1 To UBound(values, 2)
            temp(i, j) = values(i, j)
        Next j
    Next i

    values = temp
End Sub