Подкатегория вне диапазона в методе Рунге Кутты

-1

Я программирую метод Runge Kutta с адаптивным размером шага в VBA, и я столкнулся с ошибкой 9 «Subscript out the range». Может кто-то, пожалуйста, помогите мне выяснить, почему и как это исправить?Подкатегория вне диапазона в методе Рунге Кутты

Я прилагаю три отдельные подпрограммы, которые мне нужно было написать, а также начало макроса, который запускает все три программы.

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)

источник

2013-11-18 user3002894

У вас есть опция 'Base 0' и вы сделали, что ваши массивы основаны на '0'? Может ли _ensure_ это с 'ReDim ak2 (от 0 до n-1)' и т. Д. – ja72

Можете ли вы указать, где в коде вы получаете индекс за пределами диапазона? – ja72

Вы пробовали запустить код в режиме отладки? –

Одно из мест, где появляется ошибка, это с Y (I) в:

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

Поскольку Y() не инициализирован данными. Я предполагаю, что вы хотите предоставить некоторые значения, соответствующие значениям x.

источник

2013-11-18 01:02:50 NoChance

В 'Sub RungeKutta()' начальные значения 'y()' должны быть инициализированы. Это то, что вы упоминаете? – ja72

Я только что заметил, что VBA не имеет функции Pi() (называемой функцией Derivs()). Вы должны указать где-то

Public Const PI As Double = 3.14159265358979

и использовать его вместо Pi().

Также размещен Option Base 0 в верхней части кода.

С учетом этих изменений и добавления объявления для h, ytemp(), yerr() в rkqs() я получил его работать с этими значениями:

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()) 
...

источник

2013-11-18 02:17:22 ja72

Подкатегория вне диапазона в методе Рунге Кутты

ответ

Смежные вопросы