# 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: center{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 center{`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 center{`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 center{`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 center{`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)*of`n`

variables, where`y`

is an`n`

-vector. Argument`y`

is a vector of initial conditions of length`n`

which upon exit contains the solution at`x1 + h`

,`n`

is the number of dependent variables,`x1`

is the initial point,`h`

is the step size, and`derivs`

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 as`rk4(y, n, x1, h, derivs)`

except that you must provide 4 scratch arrays`t1`

-`t4`

of size`n`

. 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)*of`n`

variables, where`y`

is an`n`

-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)*of`n`

variables, where`y`

is an`n`

-vector. Starting with`y`

at`x1`

, this function 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`

. 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)*of`n`

variables, where`y`

is an`n`

-vector using a 4-th order Runge-Kutta method. This function takes a 5-th order Runge-Kutta`step`

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)*of`n`

variables, where`y`

is an`n`

-vector using a 4-th order Runge-Kutta method. This function takes a 5-th order Runge-Kutta`step`

with monitoring of local truncation to ensure accuracy and adjust stepsize. For details, see con{NumericalOrdinaryDifferentialEquations}.