2013-11-23 2 views
0

미분 방정식의 시스템을 해결하기 위해 미분 함수를 정의하려고하지만이 서브 루틴을 호출하는 실제 매크로를 실행할 때 계속 실행 오류 5가 발생합니다. 프로 시저 호출이나 인수가 잘못되었습니다. 이 오류는 x가 1보다 크지 않고 제공된 수식을 사용하여 Qv를 계산할 때 If 문에서 발생합니다. 디버깅을하는 동안 모든 변수에 대한 값이 있지만 오류가 발생하며 이유를 알 수 없습니다. 누군가 제발 도와 드릴까요?런타임 오류 5?

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 * 3.14 * 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)/(3.14 * h0 * Di^2)) 

End Sub 

음이 서브 루틴을 호출하는 매크로는 다음과 같습니다

Sub RungeKutta() 

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

Const n As Integer = 2 
Dim htry As Double 
Const eps As Double = 0.00000001 
Dim x1 As Double 
Dim x2 As Double 
Const nvar As Integer = 2 
Dim h1 As Double 
Const hmin As Double = 0.001 

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

x = 0 

y(0) = 1# 
y(1) = 1# 

Call Derivs(x, y(), dydx()) 

Call odeint(ystart(), nvar, x1, x2, eps, h1, hmin, NOk, NBad) 

' I have a bunch of coding to input the calculations into a spreadsheet that I am omitting 

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 Double, hnext As Double) 

NM1 = n - 1 

Dim ytemp() As Double 
Dim yerr() As Double 
Dim h As Double 
Const Tiny As Double = 10^(-30) 

ReDim ytemp(NM1) 
ReDim yerr(NM1) 

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 NM1 
      yscal(I) = Abs(y(I)) + Abs(h * dydx(I)) + Tiny 
     Next I 

    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 Integer, NBad As Integer) 

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

Dim y() As Double 
Dim yscal() As Double 
Dim dydx() As Double 
Dim x As Double 
Dim h As Double 
Dim hdid As Double 
Dim hnext As Double 
Const n As Integer = 2 

NM1 = n - 1 

nvar = 2 

ReDim y(NM1) 
ReDim dydx(NM1) 
ReDim yscal(NM1) 


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

kmax = 500 

ReDim xp(kmax) 
ReDim yp(NM1, kmax) 

dxsav = (x2 - x1)/kmax 


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 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 * k3(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 

이것은 Runge Kutta 방법입니다.

그래서 세 개의 프로그램을 각각 RKCK로 시작하여 디버깅 한 다음 RKQS로 이동 한 다음 모든 매개 변수가 포함 된 각 테스트 매크로를 작성하여 ODEINT로 보내고 각 프로그램과 함께 계산 된 값을 메시지 상자에 출력했습니다 각 프로그램이 예를 들어 완벽하게 작동

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


dydx(0) = -2 * x * y(0) 

dydx(1) = -3 * y(1) * x^2 

End Sub 

그래서 난 실제 문제 설명 방정식과 각 테스트 매크로를 테스트하기로 결정하고, 방정식의 집합 다음 예제했다. RKCK는 잘 작동 했으므로 RKQS도 성공했습니다. 그런 다음 ODEINT에 도착했을 때 오류 메시지가 나타났습니다.

+0

죄송합니다, 늦었습니다 :) – Sico

답변

2

런타임 오류 5는 "잘못된 프로 시저 호출"오류입니다.

내가 그 라인만큼 Y 배열 오류를 생산할 수있는 방법을 볼 수 없습니다 인덱스 당신은 어떤없이 실행 다음과 유사한이 함수를 호출의 예를 줄 필요 1.

에서 값을 가진다 오류.

Sub test() 
Dim dydx(0 To 1) As Double 
Dim y(0 To 1) As Double 
dydx(0) = 1 
dydx(1) = 2 
y(0) = 1 
y(1) = 2 
Derivs 0.5, y, dydx 

End Sub 

나는 당신의 편집 코드를 실행했는데 오류가

Qv = k * (h0^0.5) * (1 - x) * (y(1) - hstar/h0)^0.5 

발생 때 변수 값은 다음과 같습니다이 의미

y(1) = 0 
hstar = 38.3 
h0 = 80 

:

(y(1) - hstar/h0) = -0.478857734838603 

로를 Jean-François Corbett 언급, squa VBA에서 -ve 번호의 루트를 지원하지 않으므로 런타임 오류 5가 발생합니다.

+1

+1''y 배열에 인덱스 1의 값이 있습니다. '값은'hstar/h0'보다 커야합니다 :) 만약 하나가 있다면 오류가 발생합니다 'y'의 배열 원소 중 하나가 초기화되지 않았거나 그 값이'hstar/h0'보다 '>'가 아니다. –

+0

내가 사용한 코딩과 테스트 방정식을 모두 포함하도록 원래 게시 된 질문을 업데이트했다. 고맙습니다. – user3002894

+0

제 평가가 정확하다는 것을 확인해 주셔서 감사합니다 ... –

0

아마 네거티브의 제곱근을 사용하고있을 것입니다.

x^0.5x이 음수이면 "잘못된 프로 시저 호출 또는 인수가 잘못되었습니다"오류가 표시됩니다.

디버그 모드에서 코드를 확인하여이를 확인하십시오.