Introduction

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

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

The ODE Function (PV-WAVE Advantage) solves the IVP for ODEs of the form:

with yi(t = t0) specified. Here, fi is a user-supplied function that must be evaluated at any set of values (t, y1, ..., yN), i = 0, ..., N.

The previous problem statement is abbreviated by writing it as a system of first-order ODEs, y(t) = [y1(t), ..., yN(t) ]T , f(t, y) = [f1(t, y), ..., fN(t, y)]T, so that the problem becomes y' = f(t, y) with initial values y(t0).

is said to be stiff if some of the eigenvalues of the Jacobian matrix:

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 the Runge-Kutta-Verner fifth-order and sixth-order method, are inefficient or they 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, keyword R_K_V should not be specified when using the ODE Function (PV-WAVE Advantage). For more about stiff systems, see Gear (1971, Chapter 11) or Shampine and Gear (1979).

The ODEBV Function (PV-WAVE Advantage) solves the boundary value problem (BVP) for first order systems of the form:

subject to the boundary conditions:

Both functions f, h are user-supplied. The routine assumes that the user has embedded the problem into a one-parameter family of problems. In this formulation, p is an optional continuation parameter. It can be useful in solving nonlinear problems. When used, p = 0 corresponds to an easy-to-solve problem and p = 1 corresponds to the actual problem.

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:

The gi are user-supplied functions. The system is abbreviated as:

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

The DEA_PETZOLD_GEAR Procedure (PV-WAVE Advantage) 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.

See the Introduction to PDE_1D_MG for greater details about the PDE_1D_MG Procedure (PV-WAVE Advantage). This software is a variable grid-variable order integrator. It solves a problem:

with boundary conditions:

The PDE_MOL Function (PV-WAVE Advantage) solves the IVP problem for systems of the form:

subject to the boundary conditions:

and subject to the initial conditions:

for i = 1, ..., N. Here, fi, gi,:

are user-supplied, j = 1, 2.

The routine FEYNMAN_KAC Function (PV-WAVE Advantage) solves a single equation on a finite interval [xmin, xmax]. This equation often arises in applications from financial engineering and that is the primary focus of the document examples. The equation, initial conditions and boundary values are given by:

The solution is approximated by a piece-wise series of Hermite quintic polynomials on a grid of the interval [xmin, xmax] that yields a twice differentiable solution. To assist in the evaluation of the approximate solution and its derivatives there is the function FEYNMAN_KAC_EVALUATE Function (PV-WAVE Advantage).

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

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

The section describes an algorithm and a corresponding integrator routine PDE_1D_MG Procedure (PV-WAVE Advantage) for solving a system of partial differential equations.

This software is a one-dimensional differential equation solver. It requires the user to provide initial and boundary conditions in addition to a function for the evaluation of ut. The integration method is noteworthy due to the maintenance of grid lines in the space variable, x. Details for choosing new grid lines are given in Blom and Zegeling, (1994). The class of problems solved with PDE_1D_MG Procedure (PV-WAVE Advantage) is expressed by Equation 1 and given in more detail by:

The vector u ≡ [u1, ..., uNPDE]T is the solution. The integer value NPDE ≥ 1 is the number of differential equations. The functions Rj and Qj can be regarded, in special cases, as flux and source terms. The functions u, Cj,k, Rj, Qj are expected to be continuous. Allowed values for the integer m are any of m = 0,1,2. These are respectively for problems in Cartesian, cylindrical or polar, and spherical coordinates. In the two cases with m > 0, the interval [xL, xR] must not contain x = 0 as an interior point.

In the boundary conditions the functions βj and γj are continuous. In the two cases with m > 0, with an endpoint of [xL, xR] at 0, the finite value of the solution at x = 0 must be ensured. This requires the specification of the solution at x = 0, or it implies that

. The initial values satisfy u(x,t0) = u0(x), x ∈ [xL, xR], where u0 is a piece-wise continuous vector function of x with NPDE components.

The user must pose the problem so that mathematical definitions are known for the functions:

These functions are provided to the routine PDE_1D_MG in the form of two user-supplied functions. This form of the usage interface is explained below and illustrated with several examples. u0 can be supplied as the input argument u or by an optional user-supplied function. Users comfortable with the description of this algorithm may skip directly to the Examples for the PDE_1D_MG Procedure (PV-WAVE Advantage).

Equation 1 is approximatedat N = ngrids time-dependent grid values xL = x0 < x1 < ... < xi(t) < ... < xN+1 = xR. Using the total differential

, we can transform the differential equation to the form:

Using central divided differences for the factor ux leads to the system of ordinary differential equations in implicit form:

The terms Ui, Fi respectively represent the approximate solution to the partial differential equation and the value of f(u,x,t) at the point (u,x,t) = (Ui,xi(t),t). The truncation error from this approximation is second-order in the space variable x. The above ordinary differential equations are underdetermined, so additional equations are added for determining the time-dependent grid points. These additional equations contain parameters that can be adjusted by the user. Often it will be necessary to modify these parameters to solve a difficult problem. For this purpose the following quantities are needed:

The values ni are the so-called point concentration of the grid. The parameter κ ≥ 0 denotes a spatial smoothing value. Now the grid points are defined implicitly so that:

The parameter τ ≥ 0 denotes a time-smoothing value. If the value τ is chosen to be large, this results in a fixed spatial grid. Increasing τ from its default value avoids the error condition where grid lines cross. The divisors are defined by:

The value κ determines the level of clustering or spatial smoothing of the grid points. Decreasing κ from its default values also decreases the amount of spatial smoothing. The parameters Mi approximate arc length and help determine the shape of the grid or xi distribution. The parameter τ prevents the grid movement from adjusting immediately to new values of the Mi, thereby avoiding oscillations in the grid that cause large relative errors in the solution. This is important when applied to solutions with steep gradients.

The discrete form of the differential equation and the smoothing equations are combined to yield the implicit system of differential equations:

This is usually a stiff differential-algebraic system. It is solved using the integrator DEA_PETZOLD_GEAR Procedure (PV-WAVE Advantage), documented in this chapter. If DEA_PETZOLD_GEAR is needed during the evaluations of the differential equations or boundary conditions, it must be done in a separate thread to avoid possible problems with the PDE_1D_MG’s internal use of DEA_PETZOLD_GEAR. The only options for DEA_PETZOLD_GEAR set by PDE_1D_MG are the Maximum BDF Order, and the absolute and relative error values, documented as the Max_bdf_order and Atol_rtol_arrays keywords in the FEYNMAN_KAC Function (PV-WAVE Advantage).