Who's good with math and vba?

[nilodor]

Apparently I'm not anymore.

I've got a series of equations that I'm trying to evaluate at work using vba and it's just not working out for me. Most of the evaluation I can get to work out, it's just the velocity term that is really getting away from me.

The velocity (dv/dt) term is shown below:

The biggest contributions to the velocity term is the Av term and the Can term.

K is a function of dv/dt where K = Ko*sin(theta) and is shown below.


and x can be found through dx/dt = v

I am evaluating the velocity term using a 4th order runge-kutta function with a step of 0.01 as v = f(t,v).

There is a term to check the growth of the error which I am evaluating using Simpson's 1/3 rule, which is shown below:


I have 1 series of data so I've been able to check my Can and Caer terms and they are giving the right numbers. The velocity I obtain from the dv/dt is much larger than it should be. The velocity I obtain from the check of the velocity is slightly higher than it should be. I've also included my code. I'm not really sure if I'm doing the RK method right as I've never evaluated an equation with the d/dt terms on the RHS before and am unsure how to approach them.

I'm sorry if this message is a big jumbled, let me know if there's anything you need to clarify it. I'm getting pretty frustrated as I been working on this off and on for a few days now and am really not sure what I'm doing, let alone what I am doing wrong.

Thanks for any help you can provide!

Sub VelocityEval()

Dim x As Variant
Dim Can As Variant
Dim Pmx As Variant
Dim Lam As Variant
Dim psi As Variant
Dim rn As Variant
Dim kn As Variant
With Sheets("Sheet1")
Can = Array(0, 0, 0) ' Anaerobic Energy
kn = Array(0, 0, 0, 0) ' k for R-K
Pmx = Array(.Cells(3, 18).Value, .Cells(3, 19).Value, .Cells(3, 20).Value) ' Maximum power constants
Lam = Array(.Cells(3, 21).Value, .Cells(3, 22).Value, .Cells(3, 23).Value) ' Lamda term
psi = Array(.Cells(3, 24).Value, .Cells(3, 25).Value, .Cells(3, 26).Value) ' Psi term
rn = Array(.Cells(3, 27).Value, .Cells(3, 28).Value, .Cells(3, 29).Value) ' rn = psi/lam
x = 0 ' horizontal distance
dT = 0.01 ' Time step
'Constants
R = 25
ho = 0.65
hcm = 1
grav = 9.81
Ko = 0.0029
A = 3.96
T = 0.01 ' Starting time
j = 0
out = 10.17 'First output at 10m + 0.17 starting gap
cout = 7
Vw = .Cells(3, 32).Value 'Wind term
intT = 0
hvo = 0
Cano = 0

Do While x < 100.17
' Calculate Aerobic Term
Caer = R * T - R * (1 - Exp(-1 * Lam(2) * T)) / Lam(2)
'Calculate Anerobic Term
For i = 0 To 2 ' Loop through 3 Pmx, psi, lam and rn's
    Temp1 = (Pmx(i) * rn(i) ^ (1 / rn(i))) * ((1 + rn(i)) / rn(i)) ^ ((1 + rn(i)) / rn(i))
    Temp2 = (rn(i) / (Lam(i) + psi(i)))
    Temp3 = Exp(-1 * Lam(i) * T) * (psi(i) + Lam(i) * (1 - Exp(-1 * psi(i) * T))) / (Lam(i) * (Lam(i) + psi(i)))
    Can(i) = Temp1 * (Temp2 - Temp3)
  
Next i

' Sum of the Anaerobic term
CanT = Can(0) + Can(1) + Can(2)

If CanT < 0 Then CanT = 0 'Anaerobic term must always be +ve

'For the first evaluation v = the following
If j = 0 Then
    vel = (2 * (CanT + Caer)) ^ 0.5
    dx = vel * dT ' Change is horizontal distance
    K = Ko
    
Else
    If x < 5 Then hv = hcm - (hcm - ho) * (5 - x) / 5 ' Change in hv
    delT = 0 ' change in T for RK method where v = fn(T,v)
    kmult = 0 ' term to be added to velocity term in RK method
    K = Ko * Sin(Atn(grav / (delV + K * (vi - Vw) ^ 2))) ' Calc K based on change in theta
    
    'RK evaluation for dv/dt
    For d = 0 To 3
        
        kn(d) = (-1 / (vi + kmult)) * (A * (vi + kmult) + K * (vi + kmult) * ((vi + kmult) - Vw) ^ 2 + grav * _
        (hv - ho) - CanT - R * (1 - Exp(-Lam(2) * (T + delT))))
        'Set the variation for the next RK loop for time
        If d = 2 Then delT = dT Else delT = dT / 2
        'and for velocity
        kmult = kn(d) * delT
    Next d
    'Calculate the new velocity from RK
    vel = vi + (1 / 6) * (kn(0) + 2 * (kn(1) + kn(2)) + kn(3)) * dT
    
    'Use simpson's rule to evaluate the error in the velocity and update the velocity term
    int1 = K * vel * ((vel - Vw) ^ 2) * T - grav * (hv - ho)
    int2 = K * vel * ((vel - Vw) ^ 2) * (T + dT / 2) - grav * (hv - ho)
    int3 = K * vel * ((vel - Vw) ^ 2) * (T + dT) - grav * (hv - ho)
    intT = dT * (int1 + 4 * int2 + int3) / 6
    vel = (2 * (CanT + Caer - A * x - intT)) ^ 0.5
    'Change in horizontal distance
    dx = vel * dT
End If

x = x + dx ' New horizontal distance
delV = vel - vi ' Change in velocity for evaluation of new K term
vi = vel ' Velocity for next iteration of RK
 


j = j + 1
'Output every 10m
If x > out Then
    .Cells(3, cout).Value = T
    .Cells(4, cout) = vel
    
    
    For m = 0 To 2
        .Cells(6 + m, cout).Value = Can(m)
    Next m
    .Cells(9, cout).Value = Caer
    cout = cout + 1
    out = out + 10
End If
T = T + dT


Loop

End With


End Sub

[gstelmack]

Originally Posted by nilodor

...using vba...

I'm out.

[nilodor]

Originally Posted by gstelmack

I'm out.

Well maybe you don't need to be good at vba. I think I'm just really screwing up the runge kutta integration for the velocity term. Basically I'm taking v and replacing it with (vi + kh) where vi is the velocity from the previous time step, and kh is the runge kutta k times the step size and coefficient term. Then for time, t is replaced by t+h where h is the coefficient of the time step. I'm really not sure what to do with the dCan/dt or dhv/dt terms and I'm not entirely sure why the other estimation of velocity using simpsons rule is out as well, or why I can't just use that term instead (well if it was giving me the right answer).

[nilodor]

Never mind, apparently I can't take derivatives any more either! I was screwing up the dhv and dCan terms.

Math is a cruel mistress when you forget the rules.