nilodor
09-04-2012, 12:03 PM
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:
http://i1216.photobucket.com/albums/dd375/iaingidley/velterm.png
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.
http://i1216.photobucket.com/albums/dd375/iaingidley/kterm.png
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:
http://i1216.photobucket.com/albums/dd375/iaingidley/checkv.png
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
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:
http://i1216.photobucket.com/albums/dd375/iaingidley/velterm.png
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.
http://i1216.photobucket.com/albums/dd375/iaingidley/kterm.png
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:
http://i1216.photobucket.com/albums/dd375/iaingidley/checkv.png
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