odeRungeKutta¶

Solves an initial-value problem for ordinary differential equations using the Runge-Kutta-Verner fifth-order and sixth-order method.

Synopsis¶

odeRungeKuttaMgr (task, state)

odeRungeKutta (t, tend, y, state, fcn)

Required Arguments for odeRungeKuttaMgr¶

int task (Input): This function must be called with task set to ODE_INITIALIZE to set up for solving an ODE system and with task equal to ODE_RESET to clean up after it has been solved.
void state (Input/Output): The current state of the ODE solution is held in a structure pointed to by state. It cannot be directly manipulated.

Required Arguments for odeRungeKutta¶

float t (Input/Output)

Independent variable. On input, t is the initial independent variable value. On output, t is replaced by tend, unless error conditions arise.

float tend (Input)

Value of t at which the solution is desired. The value tend may be less than the initial value of t.

float y[] (Input/Output)

Array with neq components containing a vector of dependent variables. On input, y contains the initial values. On output, y contains the approximate solution.

void state (Input/Output)

The current state of the ODE solution is held in a structure pointed to by state. It must be initialized by a call to odeRungeKuttaMgr. It cannot be directly manipulated.

void fcn (neq, t, y, yprime)

User-supplied function to evaluate the right-hand side where float yprime (Output)

Array with neq components containing the vector y′. This function computes

\[\mathtt{yprime} = \frac{dy}{dt} = y' = f(t,y)\]

and neq, t, and y are defined immediately preceding this function.

Optional Arguments¶

tol, float (Input)

Tolerance for error control. An attempt is made to control the norm of the local error such that the global error is proportional to tol.

Default: tol = 100.0\*machine(4)

initialStepsize, float (Input)

Initial value for the step size h. Steps are applied in the direction of integration.

Default: initialStepsize = 0.001|tend − t|

minStepsize, float (Input)

Minimum value for the step size h.

Default: minStepsize - 0.0.

maxStepsize, float (Input)

Maximum value for the step size h.

Default: maxStepsize = 2.0.

maxNumberSteps, int (Input)

Maximum number of steps allowed.

Default: maxNumberSteps = 500.

maxNumberFcnEvals, int (Input)

Maximum number of function evaluations allowed.

Default: maxNumberFcnEvals = No enforced limit.

scale, float (Input)

A measure of the scale of the problem, such as an approximation to the Jacobian along the trajectory.

Default: scale = 1.

norm, int (Input)

Switch determining the error norm. In the following, \(e_i\) is the absolute value of the error estimate for \(y_i\).

`norm`	Error norm used
0	minimum of the absolute error and the relative error, equals the maximum of \(e_i / \max(\|y_i\|, 1)\) for i = 1, …, `neq`.
1	absolute error, equals \(max_ie_i\).
2	\(max_i(e_i ∕ w_i)\) where \(w_i = max(\|y_i\|, \textit{floor})\). The value of floor is reset using `floor`.

Default: norm = 0.

floor, float (Input)

This is used with norm. It provides a positive lower bound for the error norm option with value 2.

Default: floor = 1.0.

nstep, int nstep (Output)

Returns the number of steps taken.

nfcn (Output)

Returns the number of function evaluations used.

htrial (Output)

Returns the current trial step size.

Description¶

The function odeRungeKutta finds an approximation to the solution of a system of first-order differential equations of the form

\[\mathtt{yprime} = \frac{dy}{dt} = y' = f(t,y)\]

with given initial conditions for y at the starting value for t. The function attempts to keep the global error proportional to a user-specified tolerance. The proportionality depends on the differential equation and the range of integration.

The function odeRungeKutta is efficient for nonstiff systems where the evaluations of f(t, y) are not expensive. The code is based on an algorithm designed by Hull et al. (1976, 1978). It uses Runge-Kutta formulas of order five and six developed by J.H. Verner.

Examples¶

Example 1¶

This example solves

\[\frac{dy}{dt} = -y\]

over the interval [0, 1] with the initial condition \(y(0) = 1\). The solution is \(y(t) = e^{-t}\).

The ODE solver is initialized by a call to odeRungeKuttaMgr with ODE_INITIALIZE. This is the simplest use of the solver, so none of the default values are changed. The function odeRungeKutta is then called to integrate from t = 0 to t = 1.

from __future__ import print_function
from numpy import *
from pyimsl.math.odeRungeKutta import odeRungeKutta
from pyimsl.math.odeRungeKuttaMgr import odeRungeKuttaMgr, ODE_INITIALIZE, ODE_RESET


def fcn(neq, t, y, yprime):
    yprime[0] = -y[0]


# Initialize the ODE solver
state = []
odeRungeKuttaMgr(ODE_INITIALIZE, state)

# Integrate from t=0 to tend=1
t = [0.0]
tend = 1.0
y = [1.0]
odeRungeKutta(t, tend, y, state, fcn)

# Print the solution and error
print("y[", t[0], "] = ", y[0])
print("Error is: ", exp(-tend) - y[0])

Output¶

y[ 1.0 ] =  0.3678794411714429
Error is:  -5.551115123125783e-16

Example 2¶

Consider a predator-prey problem with rabbits and foxes. Let r be the density of rabbits, and let f be the density of foxes. In the absence of any predator-prey interaction, the rabbits would increase at a rate proportional to their number, and the foxes would die of starvation at a rate proportional to their number. Mathematically, the model without species interaction is approximated by the equation

\[r' = 2r\]

\[ƒ' = ƒ\]

With species interaction, the rate at which the rabbits are consumed by the foxes is assumed to equal the value 2rf. The rate at which the foxes increase, because they are consuming the rabbits, is equal to rf. Thus, the model differential equations to be solved are

\[r' = 2r − 2rƒ\]

\[ƒ' = -ƒ + rƒ\]

For illustration, the initial conditions are taken to be \(r(0) = 1\) and \(f(0) = 3\). The interval of integration is \(0 \leq t \leq 10\). In the program, \(y[0] = r\) and \(y[1] = f\). The ODE solver is initialized by a call to odeRungeKuttaMgr. The error tolerance is set to 0.0005. Absolute error control is selected by setting norm to the value one. We also request that nstep be set to the current number of steps in the integration. The function odeRungeKutta is then called in a loop to integrate from \(t = 0\) to \(t = 10\) in steps of \(\delta t = 1\). At each step, the solution is printed. Note that nstep is updated even though it is not an argument to this function. Its address has been stored within odeRungeKuttaMgr into the area pointed to by state. The last call to odeRungeKuttaMgr with ODE_RESET releases workspace.

from __future__ import print_function
from numpy import *
from pyimsl.math.odeRungeKutta import odeRungeKutta
from pyimsl.math.odeRungeKuttaMgr import odeRungeKuttaMgr, ODE_INITIALIZE, ODE_RESET


def fcn(neq, t, y, yprime):
    # Density change for Rabbits:
    yprime[0] = 2.0 * y[0] * (1 - y[1])
    # Density change for Foxes:
    yprime[1] = -y[1] * (1 - y[0])


t = [0.0]
y = array([1.0, 3.0])
state = []
nstep = []

# Initialize the ODE solver
odeRungeKuttaMgr(ODE_INITIALIZE, state, tol=0.005, nstep=nstep, norm=1)
print("Start       End      Density of   Density of       Number of")
print("Time        Time       Rabbits      Foxes           Steps\n")
for k in range(0, 10):
    tend = k + 1
    odeRungeKutta(t, tend, y, state, fcn)
    print('%3d %12.3f %12.3f %12.3f %12d' %
          (k, t[0], y[0], y[1], nstep[0].value))
odeRungeKuttaMgr(ODE_RESET, state)

Output¶

Start       End      Density of   Density of       Number of
Time        Time       Rabbits      Foxes           Steps

      1.000        0.082        1.468            3
      2.000        0.089        0.582            5
      3.000        0.304        0.253            6
      4.000        1.496        0.195            7
      5.000        3.901        1.535            9
      6.000        0.167        2.186           12
      7.000        0.070        0.880           14
      8.000        0.162        0.358           15
      9.000        0.721        0.190           16
     10.000        3.361        0.408           17

Fatal Errors¶

`IMSL_ODE_TOO_MANY_EVALS`	Completion of the next step would make the number of function evaluations #, but only # evaluations are allowed.
`IMSL_ODE_TOO_MANY_STEPS`	Maximum number of steps allowed, #, used. The problem may be stiff.
`IMSL_ODE_FAIL`	Unable to satisfy the error requirement. “`tol`” = # may be too small.
`IMSL_STOP_USER_FCN`	Request from user supplied function to stop algorithm. User flag = “#”.