Run time error 5????

I am trying to define a derivative function for solving a system of differential equations, however when I run the actual macro that calls up this sub routine, i keep getting Run-time error 5: Invalid procedure call or argument. This error occurs in the If statement when x is not greater than 1 and Qv is to be calculated using the equation provided. While stepping into to debug, there are values for all of the variables, but it gives me the error and I have no idea why. Can someone help please?


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

Code tags please (that's the # on the menu and paste between the [ CODE ] and [ /CODE ]

Without seeing the calling sub, I would take a wild guess and that you're

a. passing y and dydx as arrays and you don't want to, or
b. that should be y(0), or
c. something else

Code:

---

Option Explicit
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(0) - hstar / h0) ^ 0.5  '  MAYBE y(0)
 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

---

If you put a breakpoint of the line, you can hover the mouse and see which variable is not working as expected

Otherwise, post the calling code and it'd be easier to assist

Paul

Well the code that calls this sub routine is:

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

Which calls this sub routine:

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
which then calls this subroutine:

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

It's the Runge Kutta method.

possibly something in the data?


I did simple driver oprogram and changes the 'x' values and had no run time errors

Code:

---

Option Explicit
Sub drv()
    Dim y() As Double, dydx() As Double
   
    ReDim y(1)
    ReDim dydx(1)
   
    y(0) = 12.3
    y(1) = 23.4
   
    Call Derivs(0, y, dydx)
   
    MsgBox dydx(0)
    MsgBox dydx(1)
   
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 * 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

---

Sorry,

Paul