Runge Kutta 方法中的下标超出范围答案

【问题标题】：Subscript out of range in Runge Kutta methodRunge Kutta 方法中的下标超出范围
【发布时间】：2013-11-18 00:48:50
【问题描述】：

我在 VBA 中使用自适应步长对 Runge Kutta 方法进行编程，我遇到了错误 9“下标超出范围”。有人可以帮我弄清楚为什么以及如何解决它吗？

我附上了我需要编写的三个单独的子例程以及运行所有三个程序的宏的开头。

Sub rkck(y() As Double, dydx() As Double, n As Integer, x As Double, h As Double, yout() As Double, yerr() As Double)

Dim NM1 As Integer
Dim I As Integer
Dim ak2() As Double
Dim ak3() As Double
Dim ak4() As Double
Dim ak5() As Double
Dim ak6() As Double
Dim ytemp() As Double

NM1 = n - 1

ReDim ak2(NM1)
ReDim ak3(NM1)
ReDim ak4(NM1)
ReDim ak5(NM1)
ReDim ak6(NM1)
ReDim ytemp(NM1)

Const A2 As Double = 1# / 5#
Const A3 As Double = 3# / 10#
Const A4 As Double = 3# / 5#
Const A5 As Double = 1#
Const A6 As Double = 7# / 8#
Const B21 As Double = 1# / 5#
Const B31 As Double = 3# / 40#
Const B32 As Double = 9# / 40#
Const B41 As Double = 3# / 10#
Const B42 As Double = -9# / 10#
Const B43 As Double = 6# / 5#
Const B51 As Double = -11# / 54#
Const B52 As Double = 5# / 2#
Const B53 As Double = -70# / 27#
Const B54 As Double = 35# / 27#
Const B61 As Double = 1631# / 55296#
Const B62 As Double = 175# / 512#
Const B63 As Double = 575# / 13824#
Const B64 As Double = 44275# / 110592#
Const B65 As Double = 253# / 4096#
Const C1 As Double = 37# / 378#
Const C3 As Double = 250# / 621#
Const C4 As Double = 125# / 594#
Const C6 As Double = 512# / 1771#
Const DC1 As Double = C1 - 2825# / 27648#
Const DC3 As Double = C3 - 18575# / 48384#
Const DC4 As Double = C4 - 13525# / 55296#
Const DC5 As Double = -277# / 14336#
Const DC6 As Double = C6 - 1# / 4#

'First Step
For I = 0 To n - 1
    ytemp(I) = y(I) + B21 * h * dydx(I)
Next I

'Second Step
Call Derivs(x + A2 * h, ytemp(), ak2())

For I = 0 To n - 1
    ytemp(I) = y(I) + h * (B31 * dydx(I) + B32 * ak2(I))
Next I

'Third Step
Call Derivs(x + A3 * h, ytemp(), ak3())

For I = 0 To n - 1
    ytemp(I) = y(I) + h * (B41 * dydx(I) + B42 * ak2(I) + B43 * ak3(I))
Next I

'Fourth Step
Call Derivs(x + A4 * h, ytemp(), ak4())

For I = 0 To n - 1
    ytemp(I) = y(I) + h * (B51 * dydx(I) + B52 * ak2(I) + B53 * ak3(I) + B54 * ak4(I))
Next I

'Fifth Step
Call Derivs(x + A5 * h, ytemp(), ak5())

For I = 0 To n - 1
    ytemp(I) = y(I) + h * (B61 * dydx(I) + B62 * ak2(I) + B63 * ak3(I) + B64 * ak4(I) + B65 * ak5(I))
Next I

'Sixth Step
Call Derivs(x + A6 * h, ytemp(), ak6())

For I = 0 To n - 1
    yout(I) = y(I) + h * (C1 * dydx(I) + C3 * ak3(I) + C4 * ak4(I) + C6 * ak6(I))
Next I

For I = 0 To n - 1
    yerr(I) = h * (DC1 * dydx(I) + DC3 * ak3(I) + DC4 * ak4(I) + DC5 * ak5(I) + DC6 * ak6(I))
Next I    
End Sub


Sub rkqs(y() As Double, dydx() As Double, n As Integer, x As Double, htry As Double, eps As Double, yscal() As Double, hdid As Integer, hnext As Integer)

Const Safety As Double = 0.9
Const PGrow As Double = -0.2
Const PShrink As Double = -0.25
Const ErrCon As Double = (5# / Safety) ^ (1# / PGrow)

h = htry
Do
    Call rkck(y(), dydx(), n, x, h, ytemp(), yerr())

    ErrMax = 0
    For I = 0 To n - 1
        If Abs(yerr(I) / yscal(I)) > ErrMax Then ErrMax = Abs(yerr(I) / yscal(I))
    Next I

    ErrMax = ErrMax / eps

    If ErrMax > 1 Then
        dummy = h
        h = Safety * h * (ErrMax ^ PShrink)

        If h < 0.1 * dummy Then
            h = 0.1 * dummy
        End If

        xNew = x + h

        If xNew = x Then MsgBox "Stepsize underflow in rkqsl", vbExclamation
        ContLoop = True

    Else
        If ErrMax > ErrCon Then
            hnext = Safety * h * (ErrMax ^ PGrow)
        Else
            hnext = 5 * h
        End If

        hdid = h

        x = x + h

        For I = 0 To n - 1
            y(I) = ytemp(I)
        Next I

        ContLoop = False
    End If

Loop While ContLoop
End Sub


Sub odeint(ystart() As Double, nvar As Integer, x1 As Double, x2 As Double, eps As Double, h1 As Double, hmin As Double, NOk As Double, NBad As Double)

Const MaxStp As Double = 10000
Const Tiny As Double = 10 ^ (-30)

x = x1
h = Sgn(x2 - x1) * Abs(h1)
NOk = 0
NBad = 0
kount = -1

For I = 0 To nvar - 1
    y(I) = ystart(I)
Next I

If kmax > 0 Then xsav = x - 2 * dxsav
    For nstp = 1 To MaxStp
        Call Derivs(x, y(), dydx())

        For I = 0 To nvar - 1
            yscal(I) = Abs(y(I)) + Abs(h * dydx(I)) + Tiny
        Next I

        If kmax > 0 Then

            If Abs(x - xsav) > Abs(dxsav) Then

                If kount < kmax - 1 Then
                    kount = kount + 1
                    xp(kount) = x

                    For I = 0 To nvar - 1
                        yp(I, kount) = y(I)
                    Next I

                    xsav = x
                End If
            End If
        End If
        If (x + h - x2) * (x + h - x1) > 0 Then h = x2 - x
            Call rkqs(y(), dydx(), nvar, x, h, eps, yscal(), hdid, hnext)
            If hdid = h Then
                NOk = NOk + 1
            Else
                NBad = NBad + 1
            End If

            If (x - x2) * (x2 - x1) >= 0 Then
                For I = 0 To nvar - 1
                    ystart(I) = y(I)
                Next I

                If Not kmax = 0 Then
                    kount = kount + 1
                    xp(kount) = x

                    For I = 0 To nvar - 1
                        yp(I, kount) = y(I)
                    Next I
                End If
                Exit Sub
            End If

            If Abs(hnext) < hmin Then MsgBox "Stepsize smaller than minimum in odeint!", vbExclamation

            h = hnext
    Next nstp

MsgBox "Too many steps in odeint", vbExclamation          
End Sub


Sub Derivs(x As Double, y() As Double, dydx() As Double)

Const g As Double = 32.1740485564
Const Hr As Double = 100
Const h0 As Double = 80
Const fm As Double = 0.024
Const L As Double = 1500
Const dp As Double = 2
Const tc As Double = 5
Const k As Double = 25.7
Const Di As Double = 5

Dim u0 As Double
Dim Qv As Double
Dim Qv0 As Double
Dim hstar As Double

u0 = ((g * h0 / ((1 / 2) * fm * (L / dp))) * ((Hr / h0) - 1)) ^ (1 / 2)

Qv0 = (u0 * Pi() * Di ^ 2) / 4

hstar = h0 - (Qv0 / k) ^ 2

If x >= 1 Then Qv = 0
Else
    Qv = k * (h0 ^ 0.5) * (1 - x) * (y(1) - hstar / h0) ^ 0.5
End If

dydx(0) = ((tc * g * h0) / (L * u0)) * (((Hr / h0) - y(1)) - ((Hr / h0) - 1) * y(0) * Abs(y(0)))

dydx(1) = ((dp / Di) ^ 2) * (u0 * tc / h0) * y(0) - ((4 * Qv * tc) / (Pi() * h0 * Di ^ 2))
End Sub


Sub RungeKutta()

Dim y() As Double
Dim dydx() As Double
Dim yout() As Double
Dim yerr() As Double
Dim yscal() As Double
Dim hdid As Integer
Dim hnext As Integer
Dim ystart() As Double
Dim hmin As Double
Dim NOk As Double
Dim NBad As Double


Call rkck(y(), dydx(), 2, 0, 2, yout(), yerr())

Call rkqs(y(), dydx(), 2, 0, 2, 0.000001, yscal(), hdid, hnext)

Call odeint(ystart(), 2, 0, 100, 0.000001, 2, hmin, NOk, NBad)

【问题讨论】：

你有Option Base 0 并且你确定你的阵列是基于0 的吗？可以通过ReDim ak2(0 to n-1) 等确保这一点。
您能否说明代码中的下标超出范围的位置？
你试过在调试模式下运行你的代码吗？

标签： vba runge-kutta

【解决方案1】：

出现错误的地方之一是 Y(I) 在：

'First Step 
For I = 0 To n - 1
    ytemp(I) = y(I) + B21 * h * dydx(I)
Next I

因为 Y() 没有用数据初始化。我假设您想提供一些与 x 值相对应的值。

【讨论】：

Sub RungeKutta() 中的初始值y() 需要初始化。这是你说的吗？

【解决方案2】：

我刚刚注意到VBA 没有Pi() 函数（从您的Derivs() 函数调用）。您需要指定某处

Public Const PI As Double = 3.14159265358979

并用它代替Pi()。

还要将Option Base 0 放在代码顶部。

通过这些更改，并在rkqs() 中添加h、ytemp()、yerr() 的声明，我让它以这些值运行：

Sub RungeKutta()

Dim y() As Double
Dim dydx() As Double
Dim yout() As Double
Dim yerr() As Double
Dim yscal() As Double
Dim hdid As Integer
Dim hnext As Integer
Dim ystart() As Double
Dim hmin As Double
Dim NOk As Double
Dim NBad As Double

ReDim y(2), yout(2), yerr(2), dydx(2)
y(0) = 0: y(1) = 100

Call Derivs(0, y, dydx)

Call rkck(y(), dydx(), 2, 0, 2, yout(), yerr())
...

【讨论】：