Newton Excel Bach, not (just) an Excel Blog

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:

1. Read the name of the ODE evaluation routine and the required input data from the spreadsheet.
2. Convert data in spreadsheet ranges to variant arrays.
3. 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.
4. Evaluate any additional values required by the AlgLib routines
5. Set up any required output arrays
6. Call the appropriate AlgLib routine to initialise the State object
7. 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.
8. Extract the desired results from the State object using the appropriate AlgLib routine
9. If necessary convert the AlgLib output arrays to a form suitable for Excel.
10. Assign the output array to the function return value.
11. 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):

ODE input and output details, click for full size view

Example 1; First Order kinetic process

Example 1 Output

Example 2; Process dependent on X and Y

Example 2 results

Example 3; damped vibration system

Example 3 output

Example 4; Undamped SHM, ball through the centre of the Earth

Example 4 results

Example 5; deflection of a cantilever strut and buckling load

Example 5 results