Usage Notes

Ordinary Differential Equations

An ordinary differential equation is an equation involving one or more dependent variables \(y_i\), an independent variable t, and derivatives of the \(y_i\) with respect to t.

In the initial-value problem (IVP), the initial or starting values of the dependent variables \(y_i\) at a known value \(t = t_0\) are given. Values of \(y_i(t)\) for \(t > t_0\) or \(t < t_0\) are required.

The functions odeRungeKutta and odeAdamsKrogh solve the IVP for ODEs of the form

\[\frac{dy_i}{dt} = y'_i = f_i\left(t,y_1, \ldots y_N\right) \phantom{...} i=1, \ldots N\]

with \(y_i(t = t_0)\) specified. Here \(f_i\) is a user-supplied function that must be evaluated at any set of values \((t, y_1, ..., y_N)\), \(i = 1, ..., N\).

This problem statement is abbreviated by writing it as a system of first-order ODEs, \(y(t) = [y_1(t), ..., y_N(t)]^T\), \(f(t, y) = [f_1(t, y), ..., f_N(t, y)]^T\), so that the problem becomes \(y' = f(t, y)\) with initial values \(y(t_0)\).

The system \(\tfrac{dy}{dt} = y' = f(t,y)\) is said to be stiff if some of the eigenvalues of the Jacobian matrix \(\left\{ \partial y'_i / \partial y_j \right\}\) are large and negative. An alternate definition is based on the disparate integration times using a non-stiff solver compared to an implicit integration solver. Frequently differential equations modeling the behavior of physical systems are stiff, such as chemical reactions proceeding to equilibrium where subspecies effectively complete their reactions in different epochs. An alternate model concerns discharging capacitors such that different parts of the system have widely varying decay rates (or time constants).

Users typically identify stiff systems by the fact that numerical differential equation solvers such as odeRungeKutta are inefficient, or else completely fail. Special methods are often required. The most common inefficiency is that a large number of evaluations of f(t, y) (and hence an excessive amount of computer time) are required to satisfy the accuracy and stability requirements of the software. In such cases, use the IMSL function odeAdamsKrogh. For more discussion about stiff systems, see Gear (1971, Chapter 11) or Shampine and Gear (1979).

The function:doc:/math/differential/odeAdamsKrogh solves systems of ordinary differential equations of order one, order two, or mixed order one and two.

Differential-algebraic Equations

Frequently, it is not possible or not convenient to express the model of a dynamical system as a set of ODEs. Rather, an implicit equation is available in the form

\[g_i\left(t,y,\ldots,y_{\mathrm{N}},y'_1, \ldots, y'_{\mathrm{N}}\right) = 0 \phantom{...} i = 1, \dots, N\]

The \(g_i\) are user-supplied functions. The system is abbreviated as

\(g\left(t, y, y' \right) = \left[ g_1\left(t, y, y' \right), \ldots, g_N\left(t, y, y' \right) \right]^T = 0\)

With initial value \(y(t_0)\). Any system of ODEs can be trivially written as a differential-algebraic system by defining

\[g\left(t,y,y'\right) = f(t,y) - y'\]

The function differentialAlgebraicEqs solves differential-algebraic systems of index 1 or index 0. For a definition of index of a differential-algebraic system, see (Brenan et al. 1989). Also, see Gear and Petzold (1984) for an outline of the computing methods used.

Partial Differential Equations

There is a section Introduction to pde1dMg in this chapter for pde1dMg with greater details. This software is a variable grid-variable order integrator. It solves a problem

\[\begin{split}\begin{array}{c} \displaystyle\sum_{k=1}^{\textit{NPDE}} C_{j,k}\left(x, t, u, u_x\right)\frac{{\partial}u^k}{{\partial}t} = x^{-m} \frac{\partial}{{\partial}x}\left(x^mR_j\left(x, t, u, u_x\right)\right) - Q_j\left(x, t, u, u_x\right), \\ j = 1, \ldots, \textit{NPDE}, \phantom{...} x_L < x < x_R, \phantom{...} t > t_0, \phantom{...} m \in { 0, 1, 2 } \end{array}\end{split}\]

with boundary conditions

\[\begin{split}\begin{array}{l} \beta_j(x,t) R_j\left(x,t,u,u_x\right) = \gamma_j\left(x,t,u,u_x\right), \\ \mathit{atx} = x_L \text{ and } x= x_R, j=1, \ldots \mathit{NPDE} \\ \end{array}\end{split}\]

The function modifiedMethodOfLines solves the IVP problem for systems of the form

\[\frac{\partial u_i}{\partial t} = f_i \left( x,t,u_1, \ldots, u_N, \frac{\partial u_1}{\partial x}, \ldots, \frac{\partial u_N}{\partial x}, \frac{\partial^2 u_1}{\partial x^2}, \ldots, \frac{\partial^2 u_N}{\partial x^2} \right)\]

subject to the boundary conditions

\[\begin{split}\begin{aligned} \alpha_1^{(i)}u_i(a) + \beta_1^{(i)}\frac{\partial u_i}{\partial x}(a) &= \gamma_1(t) \\ \alpha_2^{(i)}u_i(b) + \beta_2^{(i)}\frac{\partial u_i}{\partial x}(b) &= \gamma_2(t) \\ \end{aligned}\end{split}\]

and subject to the initial conditions \(u_i \left( x,t = t_0 \right) = g_i(x)\) ,for \(i = 1, ..., N\). Here, \(f_i\), \(g_i\), \(\alpha_j^{(i)}\), and \(\beta_j^{(i)}\) are user-supplied, \(j = 1, 2\).

The function feynmanKac solves a single equation on a finite interval \(\left[x_{\text{min}}, x_{\text{max}}\right]\) . This equation often arises in applications from financial engineering and that is the primary focus of the document examples. The equation, initial conditions and Feynman-Kac boundary values are given by

\[\begin{split}\begin{array}{c} f_t + \mu(x,t)f_x + \frac{\sigma^2(x,t)}{2}f_{xx} - \kappa(x,t)f = \phi(f,x,t), \\ f(x,T) = p(x), \left\{f_t = \frac{\partial f}{\partial t}, \mathit{etc.}\right\} \end{array}\end{split}\]

\[a(x,t)f+b(x,t)f_x+c(x,t)f_{xx}-d(x,t), \phantom{...} x=x_{\min}x_{\max}\]

The solution is approximated by a piece-wise series of Hermite quintic polynomials on a grid of the interval \(\left[x_{\text{min}}, x_{\text{max}}\right]\) that yields a twice differentiable solution. To assist in the evaluation of the approximate solution and its derivatives there is the function feynmanKacEvaluate.

The function fastPoisson2d solves Laplace’s, Poisson’s, or Helmholtz’s equation in two dimensions. This function uses a fast Poisson method to solve a PDE of the form

\[\frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} + cu = f(x,y)\]

over a rectangle, subject to boundary conditions on each of the four sides. The scalar constant c and the function f are user specified.