Stiff Systems of Ordinary Differential Equations

A "stiff" ODE is a system in which, when expressed in the form y = A x, the matrix A is nearly singular. Under these conditions, rkfixed may oscillate or be unstable.

Radau(init, x1, x2, npoints, D) Uses the RADAU5 method.
Stiffb(init, x1, x2, npoints, D, J) Uses the Bulirsch-Stoer method.
Stiffr(init, x1, x2, npoints, D, J) Uses Rosenbrock method.

Each of these functions returns a matrix in which the first column contains the npoints between x1 and x2 at which the solution to the ODE is evaluated. Subsequent columns contain the corresponding values of the solution(s) to the nth-order ODE(s), and the n − 1 derivatives of the solution(s). Either a single ODE or a system of ODE equations is allowed.

Arguments:

• init is either a vector of n real initial values or a single scalar initial value, in the case of a single ODE.
• x1 and x2 are real, scalar endpoints of the interval over which the solution to the ODE(s) will be evaluated. Initial values in init are the values of the ODE function(s) evaluated at x1.
• npoints is the integer number of points, beyond the initial point, at which the solution is to be approximated. This determines the number of rows (1 + npoints) in the matrix of solutions.
• D(x,y) is an n-element vector-valued function of the independent variable, x, and a vector of functions, y, containing the equations for the first derivatives of all unknown functions in the system of ODEs. To create this vector, cast your equation(s) with a first derivative term by itself on the left-hand-side, with no multipliers, and no higher order derivatives in the equation. For example, a single ODE of the function y(x) which contains a second derivative must be written as a system of equations in y₀(x) and y₁(x), where the first derivative of y₀ is y₁. The single-function ODE

  

  is rewritten for the solver, using vector subscripts, as

  with implied left-hand-side

• J is a function which returns an n × (n + 1) matrix of partial derivatives of the system of equations in D. The first column contains the derivatives with respect to x; the remaining columns form the Jacobian matrix (partials with respect to y_n) for the system of differential equations. For the example above, the matrix J is

  

Notes:

• The functions D and J are provided to the solver without their arguments, and the functions y_n(x) are specified in D and J without their argument. The second argument to D and J must be subscripted as a vector, even if there is only a single entry.
• An advantage to using Radau is that no symbolic Jacobian input is needed, although if J is readily available, using it tends to increase accuracy. Radau uses non-uniform step sizes internally when it solves the differential equation, adding more steps in regions of greater variation of the solution, but returns the solution at the number of equally spaced points specified in npoints.
• If you only require the solution at the last point on the integration interval, use the lowercase versions of these functions.
• You may wish to use ODE Solve Blocks to simplify the entry of your equations and initial conditions.
• Depending on the scale of the problem, and the relative step size used, you may need to reduce the value of TOL to get appropriate solutions.

QuickSheet