Or to be more correct, the AlgLib Cash-Karp Solver, Cash-Karp being a refinement of the Runge-Kutta method of solving ordinary differential equations.
Solution of differential equations is an iterative process requiring the repeated application of the solver routine followed by evaluation of the differential equation using an external routine. This is performed in AlgLib by a process known as reverse communication, whereby an initialisation routine creates a “State” object which is used to transfer data between the solver routine and the external equation evaluation routine (which is in this case a VBA function). The steps required to wrap this process in an Excel User Defined Function (UDF) are:
- Read the name of the ODE evaluation routine and the required input data from the spreadsheet.
- Convert data in spreadsheet ranges to variant arrays.
- Convert variant arrays (which will always be base 1, 2 dimension arrays) into the form required by the AlgLib routines, which will always be base 0, and may be 1 or 2 dimensional.
- Evaluate any additional values required by the AlgLib routines
- Set up any required output arrays
- Call the appropriate AlgLib routine to initialise the State object
- Loop through the AlgLib solver routine and the external ODE evaluation routine until the solver routine returns “False”, which indicates that it has finished. The VBA function named in “FuncName” is called using the Application.Run method.
- Extract the desired results from the State object using the appropriate AlgLib routine
- If necessary convert the AlgLib output arrays to a form suitable for Excel.
- Assign the output array to the function return value.
- Import the necessary AlgLib routines; in this case the only required modules are ODESolver.bas and ap.bas
This process is illustrated in the code for the UDF ODE() shown below (may take some time to display correctly).
This code and all the related AlgLib and other routines are included in the file ODESolver.zip.
Function ODE(FuncName As String, Initial As Variant, XA As Variant, _ CoeffA As Variant, Optional Eps As Double = 0.000001, _ Optional Step As Double = 0, Optional MaxIts As Long = 100) Dim M As Long, N As Long, State As ODESolverState, _ YA() As Double, XA2() As Double, YA2() As Double, i As Long, _ Rtn As Boolean, Rep As ODESolverReport, NC As Long ' The VBA function arguments are the name of the VBA ' ODE evaluation routine, followed by the required input data. ' Convert data in spreadsheet ranges to variant arrays Initial = GetArray(Initial) XA = GetArray(XA) CoeffA = GetArray(CoeffA) ' Convert variant arrays (which will always be base ' 1, 2 dimension arrays) into the form required by the AlgLib ' routines, which will always be base 0, and may be 1 or 2 dimensional. Rtn = VarAtoDouble1D_0(Initial, YA, N, NC) If N = 1 And NC > 1 Then N = NC Rtn = VarAtoDouble1D_0(XA, XA2, M, NC) ' Evaluate any additional values required by the AlgLib routines MaxIts = MaxIts * M 'Set up any required output arrays ReDim YA2(0 To M - 1, 0 To N - 1) 'Call the appropriate AlgLib routine to initialise the State object Call ODESolverRKCK(YA(), N, XA2, M, Eps, Step, State) ' Loop through the AlgLib solver routine and the external ODE ' evaluation routine until the solver routine returns "False", ' which indicates that it has finished. ' The VBA function named in "FuncName" is called using ' the Application.Run method. Rtn = True i = 0 Do While Rtn = True And i < MaxIts Rtn = ODESolverIteration(State) State.DY = Application.Run(FuncName, State.X, State.Y, CoeffA) i = i + 1 Loop ' Extract the desired results from the State ' object using the appropriate AlgLib routine Call ODESolverResults(State, M, XA2, YA2, Rep) ' If necessary convert the AlgLib output array(s) to ' a form suitable for Excel. In this case YA2 is ' a 2D base 0 array, which may be assigned to the ' function return value without further processing. ' Assign the output array to the function return value ODE = YA2 End Function
The screen shots below show examples of the application of the ODE() function (more details of the final example, buckling of a strut, will be given in the next post):
Pingback: Daily Download 21: Assorted Solvers | Newton Excel Bach, not (just) an Excel Blog
Dear Doug,
thanks for this post, which for single ODE replaces the XNumbers ODE solver.
The question is whether your solver can be applied for system of differential equations…
Yours, Miro
LikeLike
Not as it stands. The link below might be useful:
http://epublications.bond.edu.au/cgi/viewcontent.cgi?article=1130&context=ejsie
LikeLike
Miro – it looks like my previous response was wrong!
See: https://newtonexcelbach.wordpress.com/2014/02/23/rabbits-foxes-and-lorenz-attractors/
and also response to your comment at the XNumbers post.
LikeLike
Pingback: Rabbits, Foxes, and Lorenz Attractors | Newton Excel Bach, not (just) an Excel Blog
Hi
By modifying it Can I use this script for time variable spring and damping coefficient?
Thank you for your help
Bouarfa
LikeLike
Hi Bouarfa,
It is some time since I worked on this spreadsheet, but I think the answer is yes.
If the spring and damping coefficients are just linked to time you can add that into the ODEFunc3 VBA routine.
If the coefficients are linked to displacement see: https://newtonexcelbach.wordpress.com/2014/02/23/rabbits-foxes-and-lorenz-attractors/
Please let us know how you go.
LikeLike
Hi Doug
Thank you for your response (ODEFunc3)
The ode is Q” + (w/2pi)Q’ + (w^2)Q = 0
with w=(2pi/v)v’ where v=v(t) is given by a data Feed.
Bouarfa .
LikeLike