Трудности с функцией RK4 VBA. Функция возвращает только #REF! ценность - PullRequest
Новая
       66

Трудности с функцией RK4 VBA. Функция возвращает только #REF! ценность

1 голос
/ 19 апреля 2020

Мне было интересно, смогу ли я получить какую-то помощь, чтобы найти, где я ошибаюсь в своей функции. Я пытался написать функцию для использования метода RK4 (Runge-Kutta 4) для оценки ODE. Я попытался изменить несколько вещей, но независимо от того, что я делаю, функция будет возвращать только «#REF!». Я приложил мой код RK4 ниже (с соответствующим кодом для дифференциального уравнения, используемого в функции). Я пытаюсь применить метод Эйлера для первых нескольких точек, а затем перейти к использованию метода RK4. 

 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

Вот моя функция, 

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

1 Ответ

1 голос
/ 19 апреля 2020

Попробуйте: 

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

Основная проблема заключается в том, что имя RK44 относится к ячейке, поэтому Excel запутывается. Переименуйте функцию в SIMRK44, и проблема исчезнет.

Я настроил код, чтобы он стал настоящей симуляцией, в которой результаты каждого шага зависят от предыдущего шага. В вашем коде произошла путаница с yold, которая не нужна.

Также для лучшей читаемости каждого значения, которое должно быть реальным коэффициентом, а не целым числом, я использовал литералы. Так что 2# вместо 2. Внутри компилятор все равно вносит изменения, становится яснее, что 1/6# не является целочисленным делением.

Кроме того, VBA по умолчанию передает аргументы функции ByRef, что может вызвать побочные эффекты. Вы должны поставить ByVal перед каждым аргументом, чтобы показать Excel, что функция не будет изменять аргументы. Только тогда его можно использовать как UDF.

Если бы меня попросили написать интегратор RK4, я бы придерживался стандартного псевдокода, приведенного в учебниках. Проще и проще понять, что это за намерения. 

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. Кроме того, приведенный выше код может обрабатывать переменные временные шаги, так как количество итераций не фиксировано.

...