ode_runge_kutta
Solves an initial-value problem for ordinary differential equations using the Runge-Kutta-Verner fifth-order and sixth-order method.
Synopsis
#include <imsl.h>
float imsl_f_ode_runge_kutta_mgr (int task, void **state, …, 0)
void imsl_f_ode_runge_kutta (int neq, float *t, float tend, float y[], void *state, void fcn())
The type double functions are imsl_d_ode_runge_kutta_mgr and imsl_d_ode_runge_kutta.
Required Arguments for imsl_ f_ode_runge_kutta_mgr
int task (Input)
This function must be called with task set to IMSL_ODE_INITIALIZE to set up for solving an ODE system and with task equal to IMSL_ODE_RESET to clean up after it has been solved. These values for task are defined in the include file, imsl.h.
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 imsl_f_ode_runge_kutta
int neq (Input)
Number of differential equations.
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 imsl_f_ode_runge_kutta_mgr. It cannot be directly manipulated.
void fcn (int neq, float t, float *y, float *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
and neq, t, and *y are defined immediately preceding this function.
Synopsis with Optional Arguments
#include <imsl.h>
float imsl_f_ode_runge_kutta_mgr (int task, void **state,
IMSL_TOL, float tol,
IMSL_HINIT, float hinit,
IMSL_HMIN, float hmin,
IMSL_HMAX, float hmax,
IMSL_MAX_NUMBER_STEPS, int max_steps,
IMSL_MAX_NUMBER_FCN_EVALS, int max_fcn_evals,
IMSL_SCALE, float scale,
IMSL_NORM, int norm,
IMSL_FLOOR, float floor,
IMSL_NSTEP, int *nstep,
IMSL_NFCN, int *nfcn,
IMSL_HTRIAL, float *htrial,
IMSL_FCN_W_DATA, void fcn(), void *data,
0)
Optional Arguments
IMSL_TOL, float tol (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*imsl_f_machine(4)
IMSL_HINIT, float hinit (Input)
Initial value for the step size h. Steps are applied in the direction of integration.
Default: hinit = 0.001|tend − t|
IMSL_HMIN, float hmin (Input)
Minimum value for the step size h.
Default: hmin - 0.0.
IMSL_HMAX, float hmax (Input)
Maximum value for the step size h.
Default: hmax = 2.0.
IMSL_MAX_NUMBER_STEPS, int max_steps (Input)
Maximum number of steps allowed.
Default: max_steps = 500.
IMSL_MAX_NUMBER_FCN_EVALS, int max_fcn_evals (Input)
Maximum number of function evaluations allowed.
Default: max_fcn_evals = No enforced limit.
IMSL_SCALE, float scale (Input)
A measure of the scale of the problem, such as an approximation to the Jacobian along the trajectory.
Default: scale = 1.
IMSL_NORM, int norm (Input)
Switch determining the error norm. In the following, ei is the absolute value of the error estimate for yi.
norm | Error norm used |
0 | minimum of the absolute error and the relative error, equals the maximum of ei / max (|yi|, 1) for i = 1, …, neq. |
1 | absolute error, equals maxiei. |
2 | maxi(ei ∕ wi) where wi = max (|yi|, floor). The value of floor is reset using IMSL_FLOOR. |
Default: norm = 0.
IMSL_FLOOR, float floor (Input)
This is used with IMSL_NORM. It provides a positive lower bound for the error norm option with value 2.
Default: floor = 1.0.
IMSL_NSTEP, int *nstep (Output)
Returns the number of steps taken.
IMSL_NFCN, int *nfcn (Output)
Returns the number of function evaluations used.
IMSL_HTRIAL, float *htrial (Output)
Returns the current trial step size.
IMSL_FCN_W_DATA,
void fcn(
int neq,
float t,
float *y,
float *yprime,
void *data),
void *data, (Input)
User-supplied function to evaluate the right-hand side, which also accepts a pointer to data that is supplied by the user. data is a pointer to the data to be passed to the user-supplied function. See
Passing Data to User-Supplied Functions in the introduction to this manual for more details.
Description
The function imsl_f_ode_runge_kutta finds an approximation to the solution of a system of first-order differential equations of the form
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 imsl_f_ode_runge_kutta 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
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 imsl_f_ode_runge_kutta_mgr with IMSL_ODE_INITIALIZE. This is the simplest use of the solver, so none of the default values are changed. The function imsl_f_ode_runge_kutta is then called to integrate from t = 0 to t = 1.
#include <imsl.h>
#include <math.h>
void fcn (int neq, float t, float y[], float yprime[]);
int main()
{
int neq = 1; /* Number of ode’s */
float t = 0.0; /* Initial time */
float tend = 1.0; /* Final time */
float y[1] = {1.0}; /* Initial condition */
void *state;
/* Initialize the ODE solver */
imsl_f_ode_runge_kutta_mgr(IMSL_ODE_INITIALIZE, &state, 0);
/* Integrate from t=0 to tend=1 */
imsl_f_ode_runge_kutta (neq, &t, tend, y, state, fcn);
/* Print the solution and error */
printf("y[%f] = %f\n", t, y[0]);
printf("Error is: %e\n", exp( (double)(-tend) )-y[0]);
}
void fcn (int neq, float t, float y[], float yprime[])
{
yprime[0] = -y[0];
}
Output
y[1.000000] = 0.367879
Error is: -9.149755e-09
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 ≤ t ≤ 10. In the program, y[0] = r and y[1] = f. The ODE solver is initialized by a call to imsl_f_ode_runge_kutta_mgr. The error tolerance is set to 0.0005. Absolute error control is selected by setting IMSL_NORM to the value one. We also request that nstep be set to the current number of steps in the integration. The function imsl_f_ode_runge_kutta is then called in a loop to integrate from t = 0 to t = 10 in steps of δ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 imsl_f_ode_runge_kutta_mgr into the area pointed to by state. The last call to imsl_f_ode_runge_kutta_mgr with IMSL_ODE_RESET releases workspace.
#include <imsl.h>
void fcn(int neq, float t, float y[], float yprime[]);
int main()
{
int neq = 2;
float t = 0.0; /* Initial time */
float tend; /* Final time */
float y[2] = {1.0, 3.0}; /* Initial conditions */
int k;
int nstep;
void *state;
/* Initialize the ODE solver */
imsl_f_ode_runge_kutta_mgr(IMSL_ODE_INITIALIZE, &state,
IMSL_TOL, 0.0005,
IMSL_NSTEP, &nstep,
IMSL_NORM, 1,
0);
printf("\n Start End Density of Density of Number of" );
printf("\n Time Time Rabbits Foxes Steps\n\n");
for (k = 0; k < 10; k++) {
tend = k + 1;
imsl_f_ode_runge_kutta (neq, &t, tend, y, state, fcn);
printf("%3d %12.3f %12.3f %12.3f %12d\n", k, t, y[0], y[1], nstep);
}
imsl_f_ode_runge_kutta_mgr(IMSL_ODE_RESET, &state, 0);
}
void fcn (int neq, float t, float y[], float yprime[])
{
/* Density change rate for Rabbits: */
yprime[0] = 2*y[0]*(1 - y[1]);
/* Density change rate for Foxes: */
yprime[1] = -y[1]*(1 - y[0]);
}
Output
Start End Density of Density of Number of
Time Time Rabbits Foxes Steps
0 1.000 0.078 1.465 4
1 2.000 0.085 0.578 6
2 3.000 0.292 0.250 7
3 4.000 1.449 0.187 8
4 5.000 4.046 1.444 11
5 6.000 0.176 2.256 15
6 7.000 0.066 0.908 18
7 8.000 0.148 0.367 20
8 9.000 0.655 0.188 21
9 10.000 3.157 0.352 23
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 = "#". |