dkrnambiar
New Member
This is a VBA routine for Runge Kutta 4th order integration.
Function rhs(x, t, dxdt, p)
'==============================================’
' Input your right hand sides in the form dxdt(1) = expression
' in t, x(1), x(2), etc. The right hand side may also include
' several parameter values: p(1), p(2) etc
'==============================================’
dxdt(1) = -p(1) * x(2) * x(1)
dxdt(2) = p(1) * x(2) * x(1) - p(2) * x(2)
dxdt(3) = p(2) * x(2)
'=============================================’
' Function must end with rhs = 0 statement to return a value
'=============================================’
rhs = 0
End Function
Function rk4(nstep As Integer, h As Double, x As Range, p As Range)
'==================================================
'Applies nstep steps of 4th order RK method to the equations
'specified in module rhs
'h holds the timestep, x(1) holds the value of t, with
'x(2)... holding current x vals.
'p is a range of parameters that can be used in the function rhs.
'==================================================
Dim n, nr, i, step, temp As Integer
n = x.Count 'find how many entries are in array x
nr = x.Rows.Count 'used to check if input range a row or column
neqn = n - 1 'find how many independent variables in equations
Dim k(), xtemp(), ttemp, soln(), dxdt(), xcurr(), xx() As Double
ReDim k(1 To neqn), xtemp(1 To neqn), soln(1 To neqn)
ReDim dxdt(1 To neqn), xcurr(1 To neqn), xx(1 To n)
'==================================================
'Strip off the t value from the x values & store x values in xcurr
'==================================================
t = x(1)
For i = 1 To neqn: xcurr(i) = x(i + 1): Next i
For step = 1 To nstep
For i = 1 To neqn: xtemp(i) = xcurr(i): Next i
temp = rhs(xtemp, t, dxdt, p) 'evaluate rhs with current x value
For i = 1 To neqn
k(i) = h * dxdt(i) 'set up k1
soln(i) = xcurr(i) + k(i) / 6 '& add k1 contrbn to new soln pt
xtemp(i) = xcurr(i) + k(i) / 2 ' set up x for eval of k2
Next i
ttemp = t + h / 2 ' set up t val for eval of k2
temp = rhs(xtemp, ttemp, dxdt, p) 'evaluate rhs to use to eval k2
For i = 1 To neqn
k(i) = h * dxdt(i) 'eval k2
soln(i) = soln(i) + k(i) / 3 'add k2 contrbn to new soln pt
xtemp(i) = xcurr(i) + k(i) / 2 'set up x for eval of k3
Next i
temp = rhs(xtemp, ttemp, dxdt, p) ' evaluate rhs to use to eval k3
For i = 1 To neqn
k(i) = h * dxdt(i) 'evaluate k3
soln(i) = soln(i) + k(i) / 3 'add k3 contrbn to new soln pt
xtemp(i) = xcurr(i) + k(i) 'set up x for eval of k4
Next i
t = t + h 'set up t value for eval of k4
temp = rhs(xtemp, t, dxdt, p) 'evaluate rhs to use to eval k4
For i = 1 To neqn
k(i) = h * dxdt(i) 'eval k4
xcurr(i) = soln(i) + k(i) / 6 'add k4 contrbn to new soln pt
Next i
'==================================================
'End of R-K process for current step.
'xcurr contains current sol for x
'==================================================
Next step
'==================================================
'End of required number of timesteps.
'Build rk4 to hold t value followed by x values.
'==================================================
xx(1) = t: For i = 1 To neqn: xx(i + 1) = xcurr(i): Next i
If nr > 1 Then 'Test to see if x_in in cols, write soln in cols
rk4 = WorksheetFunction.Transpose(xx) 'x_in in cols
Else
rk4 = xx 'x_in in rows
End If
End Function
Spreadsheet set out to use the above VBA routine in an array function call generates an error.
[pre]
[/pre]
Function rhs(x, t, dxdt, p)
'==============================================’
' Input your right hand sides in the form dxdt(1) = expression
' in t, x(1), x(2), etc. The right hand side may also include
' several parameter values: p(1), p(2) etc
'==============================================’
dxdt(1) = -p(1) * x(2) * x(1)
dxdt(2) = p(1) * x(2) * x(1) - p(2) * x(2)
dxdt(3) = p(2) * x(2)
'=============================================’
' Function must end with rhs = 0 statement to return a value
'=============================================’
rhs = 0
End Function
Function rk4(nstep As Integer, h As Double, x As Range, p As Range)
'==================================================
'Applies nstep steps of 4th order RK method to the equations
'specified in module rhs
'h holds the timestep, x(1) holds the value of t, with
'x(2)... holding current x vals.
'p is a range of parameters that can be used in the function rhs.
'==================================================
Dim n, nr, i, step, temp As Integer
n = x.Count 'find how many entries are in array x
nr = x.Rows.Count 'used to check if input range a row or column
neqn = n - 1 'find how many independent variables in equations
Dim k(), xtemp(), ttemp, soln(), dxdt(), xcurr(), xx() As Double
ReDim k(1 To neqn), xtemp(1 To neqn), soln(1 To neqn)
ReDim dxdt(1 To neqn), xcurr(1 To neqn), xx(1 To n)
'==================================================
'Strip off the t value from the x values & store x values in xcurr
'==================================================
t = x(1)
For i = 1 To neqn: xcurr(i) = x(i + 1): Next i
For step = 1 To nstep
For i = 1 To neqn: xtemp(i) = xcurr(i): Next i
temp = rhs(xtemp, t, dxdt, p) 'evaluate rhs with current x value
For i = 1 To neqn
k(i) = h * dxdt(i) 'set up k1
soln(i) = xcurr(i) + k(i) / 6 '& add k1 contrbn to new soln pt
xtemp(i) = xcurr(i) + k(i) / 2 ' set up x for eval of k2
Next i
ttemp = t + h / 2 ' set up t val for eval of k2
temp = rhs(xtemp, ttemp, dxdt, p) 'evaluate rhs to use to eval k2
For i = 1 To neqn
k(i) = h * dxdt(i) 'eval k2
soln(i) = soln(i) + k(i) / 3 'add k2 contrbn to new soln pt
xtemp(i) = xcurr(i) + k(i) / 2 'set up x for eval of k3
Next i
temp = rhs(xtemp, ttemp, dxdt, p) ' evaluate rhs to use to eval k3
For i = 1 To neqn
k(i) = h * dxdt(i) 'evaluate k3
soln(i) = soln(i) + k(i) / 3 'add k3 contrbn to new soln pt
xtemp(i) = xcurr(i) + k(i) 'set up x for eval of k4
Next i
t = t + h 'set up t value for eval of k4
temp = rhs(xtemp, t, dxdt, p) 'evaluate rhs to use to eval k4
For i = 1 To neqn
k(i) = h * dxdt(i) 'eval k4
xcurr(i) = soln(i) + k(i) / 6 'add k4 contrbn to new soln pt
Next i
'==================================================
'End of R-K process for current step.
'xcurr contains current sol for x
'==================================================
Next step
'==================================================
'End of required number of timesteps.
'Build rk4 to hold t value followed by x values.
'==================================================
xx(1) = t: For i = 1 To neqn: xx(i + 1) = xcurr(i): Next i
If nr > 1 Then 'Test to see if x_in in cols, write soln in cols
rk4 = WorksheetFunction.Transpose(xx) 'x_in in cols
Else
rk4 = xx 'x_in in rows
End If
End Function
Spreadsheet set out to use the above VBA routine in an array function call generates an error.
[pre]
Code:
Modelling the spread of an epidemic
Parameters related to the disease Integration parameters
Infectivity 1/Durn. number of timestep
b g steps h
0.002 0.3 5 0.1
Initial values t S I R
0 843 10 20
Integration t S I R
0.00 843.00 10.00 20.00
#REF! #REF! #REF! #REF!