NumericalOrdinaryDifferentialEquations¶
numode.spad line 1 [edit on github]
This package is a suite of functions for the numerical integration of an ordinary differential equation of n
variables: centerline{dy/dx = f
(y
, x
) y
is an n
-vector} par All the routines are based on a 4-th order Runge-Kutta kernel. These routines generally have as arguments: n
, the number of dependent variables; x1
, the initial point; h
, the step size; y
, a vector of initial conditions of length n
which upon exit contains the solution at x1 + h
; derivs
, a function which computes the right hand side of the ordinary differential equation: derivs(dydx, y, x)
computes dydx
, a vector which contains the derivative information. par In order of increasing complexity: begin{items} item rk4(y, n, x1, h, derivs)
advances the solution vector to x1 + h
and return the values in y
. item rk4(y, n, x1, h, derivs, t1, t2, t3, t4)
is the same as rk4(y, n, x1, h, derivs)
except that you must provide 4 scratch arrays t1
-t4
of size n
. item Starting with y
at x1
, rk4f(y, n, x1, x2, ns, derivs)
uses ns
fixed steps of a 4-th order Runge-Kutta integrator to advance the solution vector to x2
and return the values in y
. Argument x2
, is the final point, and ns
, the number of steps to take. item rk4qc(y, n, x1, step, eps, yscal, derivs)
takes a 5-th order Runge-Kutta step with monitoring of local truncation to ensure accuracy and adjust stepsize. The function takes two half steps and one full step and scales the difference in solutions at the final point. If the error is within eps
, the step is taken and the result is returned. If the error is not within eps
, the stepsize if decreased and the procedure is tried again until the desired accuracy is reached. Upon input, an trial step size must be given and upon return, an estimate of the next step size to use is returned as well as the step size which produced the desired accuracy. The scaled error is computed as centerline{error = MAX(ABS((y2steps(i) - y1step(i))/yscal(i)))
} and this is compared against eps
. If this is greater than eps
, the step size is reduced accordingly to centerline{hnew = 0.9 * hdid * (error/eps)^(-1/4)
} If the error criterion is satisfied, then we check if the step size was too fine and return a more efficient one. If error
> eps
* (6.0E-04) then the next step size should be centerline{hnext
= 0.9 * hdid * (error
/eps
)^(-1/5
)} Otherwise hnext = 4.0 * hdid
is returned. A more detailed discussion of this and related topics can be found in the book “Numerical Recipies” by W
.Press, B
.P
. Flannery, S
.A. Teukolsky, W
.T
. Vetterling published by Cambridge University Press. Argument step
is a record of 3 floating point numbers (to_try , did , next)
, eps
is the required accuracy, yscal
is the scaling vector for the difference in solutions. On input, step.to_try
should be the guess at a step size to achieve the accuracy. On output, step.did
contains the step size which achieved the accuracy and step.next
is the next step size to use. item rk4qc(y, n, x1, step, eps, yscal, derivs, t1, t2, t3, t4, t5, t6, t7)
is the same as rk4qc(y, n, x1, step, eps, yscal, derivs)
except that the user must provide the 7 scratch arrays t1-t7
of size n
. item rk4a(y, n, x1, x2, eps, h, ns, derivs)
is a driver program which uses rk4qc
to integrate n
ordinary differential equations starting at x1
to x2
, keeping the local truncation error to within eps
by changing the local step size. The scaling vector is defined as centerline{yscal(i) = abs(y(i)) + abs(h*dydx(i)) + tiny
} where y(i)
is the solution at location x
, dydx
is the ordinary differential equation's
right hand side, h
is the current step size and tiny
is 10 times the smallest positive number representable. The user must supply an estimate for a trial step size and the maximum number of calls to rk4qc
to use. Argument x2
is the final point, eps
is local truncation, ns
is the maximum number of call to rk4qc
to use. end{items}
- rk4: (Vector Float, Integer, Float, Float, (Vector Float, Vector Float, Float) -> Void) -> Void
rk4(y, n, x1, h, derivs)
uses a 4-th order Runge-Kutta method to numerically integrate the ordinary differential equation dy/dx = f(y, x) ofn
variables, wherey
is ann
-vector. Argumenty
is a vector of initial conditions of lengthn
which upon exit contains the solution atx1 + h
,n
is the number of dependent variables,x1
is the initial point,h
is the step size, andderivs
is a function which computes the right hand side of the ordinary differential equation. For details, see NumericalOrdinaryDifferentialEquations.
- rk4: (Vector Float, Integer, Float, Float, (Vector Float, Vector Float, Float) -> Void, Vector Float, Vector Float, Vector Float, Vector Float) -> Void
rk4(y, n, x1, h, derivs, t1, t2, t3, t4)
is the same asrk4(y, n, x1, h, derivs)
except that you must provide 4 scratch arrayst1
-t4
of sizen
. For details, see con{NumericalOrdinaryDifferentialEquations}.
- rk4a: (Vector Float, Integer, Float, Float, Float, Float, Integer, (Vector Float, Vector Float, Float) -> Void) -> Void
rk4a(y, n, x1, x2, eps, h, ns, derivs)
is a driver function for the numerical integration of an ordinary differential equation dy/dx = f(y, x) ofn
variables, wherey
is ann
-vector using a 4-th order Runge-Kutta method. For details, see con{NumericalOrdinaryDifferentialEquations}.
- rk4f: (Vector Float, Integer, Float, Float, Integer, (Vector Float, Vector Float, Float) -> Void) -> Void
rk4f(y, n, x1, x2, ns, derivs)
uses a 4-th order Runge-Kutta method to numerically integrate the ordinary differential equation dy/dx = f(y, x) ofn
variables, wherey
is ann
-vector. Starting withy
atx1
, this function usesns
fixed steps of a 4-th order Runge-Kutta integrator to advance the solution vector tox2
and return the values iny
. For details, see con{NumericalOrdinaryDifferentialEquations}.
- rk4qc: (Vector Float, Integer, Float, Record(to_try: Float, did: Float, next: Float), Float, Vector Float, (Vector Float, Vector Float, Float) -> Void) -> Void
rk4qc(y, n, x1, step, eps, yscal, derivs)
is a subfunction for the numerical integration of an ordinary differential equation dy/dx = f(y, x) ofn
variables, wherey
is ann
-vector using a 4-th order Runge-Kutta method. This function takes a 5-th order Runge-Kuttastep
with monitoring of local truncation to ensure accuracy and adjust stepsize. For details, see con{NumericalOrdinaryDifferentialEquations}.
- rk4qc: (Vector Float, Integer, Float, Record(to_try: Float, did: Float, next: Float), Float, Vector Float, (Vector Float, Vector Float, Float) -> Void, Vector Float, Vector Float, Vector Float, Vector Float, Vector Float, Vector Float, Vector Float) -> Void
rk4qc(y, n, x1, step, eps, yscal, derivs, t1, t2, t3, t4, t5, t6, t7)
is a subfunction for the numerical integration of an ordinary differential equation dy/dx = f(y, x) ofn
variables, wherey
is ann
-vector using a 4-th order Runge-Kutta method. This function takes a 5-th order Runge-Kuttastep
with monitoring of local truncation to ensure accuracy and adjust stepsize. For details, see con{NumericalOrdinaryDifferentialEquations}.