Solves a system of partial differential equations of the form ut = f(x, t, u, ux, uxx) using the method of lines. The solution is represented with cubic Hermite polynomials.

Let M = npdes, N = nx and xi = xbreak(I). The function imsl_f_modified_method_of_lines uses the method of lines to solve the partial differential equation system

with the initial conditions

and the boundary conditions

for k = 1, …, M.

Cubic Hermite polynomials are used in the x variable approximation so that the trial solution is expanded in the series

where φi(x) and ψi(x) are the standard basis functions for the cubic Hermite polynomials with the knots x1 < x2 < … < xN. These are piecewise cubic polynomials with continuous first derivatives. At the breakpoints, they satisfy

According to the collocation method, the coefficients of the approximation are obtained so that the trial solution satisfies the differential equation at the two Gaussian points in each subinterval,

for j = 1, …, N. The collocation approximation to the differential equation is

for k = 1, …, M and j = 1, …, 2(N - 1).

This is a system of 2M(N - 1) ordinary differential equations in 2M N unknown coefficient functions, ak,i and bk,i. This system can be written in the matrix‑vector form as A dc/dt = F (t, c) with c(t0) = c0 where c is a vector of coefficients of length 2M N and c0 holds the initial values of the coefficients. The last 2M equations are obtained from the boundary conditions.

If αk = βk = 0, it is assumed that no boundary condition is desired for the k‑th unknown at the left endpoint. A similar comment holds for the right endpoint. Thus, collocation is done at the endpoint. This is generally a useful feature for systems of first-order partial differential equations.

The input/output array Y contains the values of the ak,i. The initial values of the bk,i are obtained by using the IMSL cubic spline function imsl_f_cub_spline_interp_e_cnd (Interpolation and Approximation) to construct functions

such that

The IMSL function imsl_f_cub_spline_value (Interpolation and Approximation) is used to approximate the values

Optional argument IMSL_INITIAL_VALUE_DERIVATIVE allows the user to provide the initial values of bk,i.

The order of matrix A is 2M N and its maximum bandwidth is 6M ‑ 1. The band structure of the Jacobian of F with respect to c is the same as the band structure of A. This system is solved using a modified version of imsl_f_ode_adams_krogh. Numerical Jacobians are used exclusively. Gear’s BDF method is used as the default because the system is typically stiff. For more details, see Sewell (1982).

Four examples of PDEs are now presented that illustrate how users can interface their problems with IMSL PDE solving software. The examples are small and not indicative of the complexities that most practitioners will face in their applications. A set of seven sample application problems, some of them with more than one equation, is given in Sincovec and Madsen (1975). Two further examples are given in Madsen and Sincovec (1979).

The normalized linear diffusion PDE, ut = uxx, 0 ≤ x ≤ 1, t > t0, is solved. The initial values are t0 = 0, u(x, t0) = u0 = 1. There is a “zero-flux” boundary condition at x = 1, namely ux(1, t) = 0, (t > t0). The boundary value of u(0, t) is abruptly changed from u0 to the value 0, for t >0.

When the boundary conditions are discontinuous, or incompatible with the initial conditions such as in this example, it may be important to use double precision.

Here, Problem C is solved from Sincovec and Madsen (1975). The equation is of diffusion-convection type with discontinuous coefficients. This problem illustrates a simple method for programming the evaluation routine for the derivative, ut. Note that the weak discontinuities at x = 0.5 are not evaluated in the expression for ut. The problem is defined as

In this example, using imsl_f_modified_method_of_lines, the linear normalized diffusion PDE ut = uxx is solved but with an optional use that provides values of the derivatives, ux, of the initial data. Due to errors in the numerical derivatives computed by spline interpolation, more precise derivative values are required when the initial data is u(x, 0) = 1 + cos[(2n - 1)πx], n > 1. The boundary conditions are “zero flux” conditions ux(0, t) = ux(1, t) = 0 for t > 0.

This optional usage signals that the derivative of the initial data is passed by the user. The values u(x, tend) and ux(x, tend) are output at the breakpoints with the optional usage.

In this example, consider the linear normalized hyperbolic PDE, utt = uxx, the “vibrating string” equation. This naturally leads to a system of first order PDEs. Define a new dependent variable ut = v. Then, vt = uxx is the second equation in the system. Take as initial data u(x, 0) = sin(πx) and ut(x, 0) = v(x, 0) = 0. The ends of the string are fixed so u(0, t) = u(1, t) = v(0, t) = v(1, t) = 0. The exact solution to this problem is u(x, t) = sin(πx) cos(πt). Residuals are computed at the output values of t for 0 < t ≤ 2. Output is obtained at 200 steps in increments of 0.01.

Even though the sample code imsl_f_modified_method_of_lines gives satisfactory results for this PDE, users should be aware that for nonlinear problems, “shocks” can develop in the solution. The appearance of shocks may cause the code to fail in unpredictable ways. See Courant and Hilbert (1962), pp 488-490, for an introductory discussion of shocks in hyperbolic systems.