User:Griffgruff/sandbox

From Wikipedia, the free encyclopedia

Method of Lines
(with an introduction to
Partial Differential Equations)

The method of lines (MOL) is a general procedure for the solution of partial differential equations (PDEs). Before discussing the MOL, we consider briefly some of the basic features of PDEs.

Contents

[edit] Some PDE Basics

Our physical world is most generally described in scientific and engineering terms with respect to three-dimensional space and time which we abbreviate as spacetime. PDEs provide a mathematical description of physical spacetime, and they are therefore among the most widely used forms of mathematics. As a consequence, methods for the solution of PDEs, such as the MOL, are of broad interest in science and engineering.

As a basic illustrative example of a PDE, we consider

\dfrac{\partial u}{\partial t}=D\dfrac{\partial^{2}u}{\partial x^{2}}\quad (1)

where

u = dependent variable
t = independent variable representing time
x = independent variable representing one
dimension of three-dimensional space
D = constant explained below>


Note that eq. (1) has two independent variables, x and t which is the reason it is classified as a PDE (any differential equation with more than one independent variable is a PDE). A differential equation with only one independent variable is generally termed an ordinary differential equation (ODE); we will consider ODEs later as part of the MOL.


Eq. (1) is termed the diffusion equation. When applied to heat transfer, it is Fourier's second law; the dependent variable u is temperature and D is the thermal diffusivity. When eq. (1) is applied to mass diffusion, it is Fick's second law; u is mass concentration and D is the coefficient of diffusion or the diffusivity.

\scriptstyle{ \frac{\partial u}{\partial t} } is a partial derivative of u with respect to t x is held constant when taking this partial derivative, which is why partial is used to describe this derivative). Eq. (1) is first order in t since the highest order partial dervtaive in t is first order; it is second order in x since the highest order derivative in x is second order. Eq. (1) is linear or first degree since all of the terms are to the first power (note that order and degree can be easily confused).

[edit] Initial and Boundary Conditions

Before we consider a solution to eq. (1), we must specify some auxiliary conditions to complete the statement of the PDE problem, i.e., to have a well posed problem. The number of required auxiliary conditions is determined by the highest order derivative in each independent variable. Since eq. (1) is first order in t and second order in x, it requires one auxiliary condition in t and two auxiliary conditions in x.


t is termed an initial value variable and therefore requires one initial condition (IC). It is an initial value variable since it starts at an initial value, t0, and moves forward over a finite interval \scriptstyle{ t_{0}\leq t\leq t_{f} } or a semi-infinite interval \scriptstyle{ t_{0}\leq t\leq\infty } without any additional conditions being imposed. Typically in a PDE application, the initial value variable is time, as in the case of eq. (1).


x is termed a boundary value variable and therefore requires two boundary conditions (BCs). It is a boundary value variable since it varies over a finite interval \scriptstyle{ x_{0}\leq x\leq x_{f} }, a semi-infinite interval \scriptstyle{ x_{0}\leq x\leq\infty } or a fully infinite interval \scriptstyle{ -\infty\leq x\leq\infty } , and at two different values of x, conditions are imposed on u in eq. (1). Typically, the two values of x correspond to boundaries of a physical system, and hence the name boundary conditions.


As examples of auxiliary conditons for eq. (1),

  • An IC could be
u(x,t=0)=u_{0} \quad (2)



where u0 is a given initial (constant) value of u for all x.
  • Two BCs could be
u(x=x_{0},t)=u_{b} \quad (3a)


\dfrac{\partial u(x=x_{f},t)}{\partial x}=0 \quad (3b)


where ub is a given boundary (constant) value of u for all t.

BCs can be of three types:

    1. If the dependent variable is specified, as in BC (3a), the BC is termed Dirichlet.
    2. If the derivative of the dependent variable is specifed, as in BC (3b), the BC is termed Neumann.
    3. If both the dependent variable and its derivative appear in the BC, it is termed a BC of the third type or a Robin BC.

[edit] Types of PDE Solutions

Eqs. (1), (2) and (3) constitute a complete (well posed) PDE problem and we can now consider what we mean by a solution to this problem. Briefly, the solution of a PDE problem is a function that defines the dependent variable as a function of the independent variables, in this case u(x,t). In other words, we seek a function that when substituted in the PDE and all of its auxiliary conditions satisfies simultaneously all of these equations.


The solution can be of two types:

  • If the solution is an actual mathematical function, it is termed an analytical solution. While analytical solutions are the "gold standard" for PDE solutions in the sense that they are exact, they are also generally difficult to derive mathematically for all but the simplest PDE problems (in much the same way that solutions to nonlinear algebraic equations generally cannot be derived mathematically except for certain classes of nonlinear equations).
  • If the solution is in numerical form, e.g., u(x,t) tabulated numerically as a function of x and t, it is termed a numerical solution. Ideally, the numerical solution is simply a numerical evaluation of the analytical solution. But since an analytical solution is generally unavailable for realistic PDE problems in science and engineering, the numerical solution is an approximation to the analytical solution, and our expectation is that it represents the analytical solution with good accuracy. However, numerical solutions can be computed with modern-day computers for very complex problems, and they will generally have good accuracy (even though this cannot be established directly by comparison with the analytical solution since the latter is usually unknown).

The focus of the MOL is the calculation of accurate numerical solutions.

[edit] PDE Subscript Notation

Before we go on to the general classes of PDEs that the MOL can handle, we briefly discuss an alternative notation for PDEs. Instead of writing the partial derivatives as in eq. (1), we adopt a subscript notation that is easier to state and bears a closer resemblance to the associated computer coding. For example, we can write eq. (1) as

u_{t}=Du_{xx} \quad (4)

where, for example, \scriptstyle{ u_{t}\Rightarrow\frac{\partial u}{\partial t} }. In other words, a partial derivative is represented as the dependent variable with a subscript that defines the independent variable. For a derivative that is of order n, the independent variable is repeated n times, e.g., for eq. (1), \scriptstyle{ u_{xx}\Rightarrow\frac{\partial^{2}u}{\partial x^{2}} }.

[edit] A General PDE System

Using the subscript notation, we can now consider some general PDEs. For example, a general PDE first order in t can be considered

\overline{u}_{t}=\overline{f}(\overline{x},t,\overline{u},\overline{u}_{\overline{x}},\overline{u}_{\overline{xx},\ldots}) \quad (5)

where an overbar (overline) denotes a vector. For example, \overline{u} denotes a vector of n dependent variables

\overline{u}=(u_{1},u_{2},\ldots u_{n})^{T}

i.e., a system of n simultaneous PDEs. Similarly, \scriptstyle{ \overline{f} } denotes an n-vector of derivative functions

\overline{f}=(f_{1},f_{2},\ldots f_{n})^{T}

where T denotes a transpose (here a row vector is transposed to a column vector). Note also that \scriptstyle{ \overline{x} } is a vector of spatial coordinates, so that, for example, in Cartesian coordinates \scriptstyle{ \overline{x}=(x,y,z)^{T} } while in cylindrical coordinates \scriptstyle{ \overline{x}=(r,\theta,z)^{T} }. Thus, eq. (5) can represent PDEs in one, two and three spatial dimensions.


Since eq. (5) is first order in t, it requires one initial condition

\overline{u}(\overline{x},t=0)=\overline{u}_{0}(\overline{x},\overline{u},\overline{u}_{\overline{x}},\overline{u}_{\overline{xx},\ldots}) \quad (6)

where \overline{u}_{0} is an n-vector of initial condition functions

\overline{u}_{0}=(u_{10},u_{20},\ldots u_{n0})^{T}

The derivative vector \scriptstyle{ \overline{f}} of eq. (5) includes functions of various spatial derivatives, (\scriptstyle{ \overline{u},\overline{u}_{\overline{x}},\overline{u}_{\overline{xx},\ldots} }), and therefore we cannot state a priori the required number of BCs. For example, if the highest order derivative in \overline{x} in all of the derivative functions is second order, then we require \scriptstyle{ 2\times n } BCs. We state the general BC requirement of eq. (5) as

\overline{f}_{b}(\overline{x}_{b},\overline{u},\overline{u}_{\overline{x}},\overline{u}_{\overline{xx},\ldots})=0\quad (7)

where the subscript b denotes "boundary". The vector of boundary condition functions, \scriptstyle{ \overline{f}_{b} } has a length (number of functions) determined by the highest order derivative in \overline{x} in each PDE (in eq. (5)) as discussed previously.

[edit] PDE Geometric Classification

Eqs. (5), (6) and (7) constitute a general PDE system to which the MOL can be applied. Before proceeding to the details of how this might be done, we need to discuss the three basic forms of the PDEs as classified geometrically. This geometric classification can be done rigorously if certain mathematical forms of the functions in eqs. (5), (6) and (7) are assumed. However, we will adopt a somewhat more descriptive (less rigorous but more general) form of these functions for the specification of the three geometric classes.


If the derivative functions in eq. (5) contain only first order derivatives in \scriptstyle{ \overline{x} }, the PDEs are classified as first order hyperbolic. As an example, the equation

u_{t}+vu_{x}=0\quad (8)

is generally called the linear advection equation; in physical applications, v is a linear or flow velocity. Although eq. (8) is possibly the simplest PDE, this simplicity is deceptive in the sense that it is can be very difficult to integrate numerically since it propagates discontinuities, a distinctive feature of first order hyperbolic PDEs.


Eq. (8) is termed a conservation law since it expresses conservation of mass, energy or momentum under the conditions for which it is derived, i.e., the assumptions on which the equation is based. Conservation laws are a bedrock of PDE mathematical models in science and engineering, and an extensive literature pertaining to their solution, both analytical and numerical, has evolved over many years.


An example of a first order hyperbolic system (using the notation u_{1}\Rightarrow u,u_{2}\Rightarrow v) is

u_{t}=v_{x} \quad (9a)


v_{t}=u_{x} \quad (9b)

Eqs. (9) constitute a system of two linear, constant coefficient, first order hyperbolic PDEs.


Differentiation and algebraic substitution can occasionally be used to eliminate some dependent variables in systems of PDEs. For example, if eq. (9a) is differentiated with respect to t and eq. (9b) is differentiated with respect to x

u_{tt}=v_{xt} \


v_{tx}=u_{xx} \

we can then eliminate the mixed partial derivative between these two equations (assuming vxt in the first equation equals vtx in the second equation) to obtain

u_{tt}=u_{xx} \quad (10)

Eq. (10) is the second order hyperbolic wave equation.


If the derivative functions in eq. (5) contain only second order derivatives in \overline{x}, the PDEs are classified as parabolic. Eq. (1) is an example of a parabolic PDE.


Finally, if a PDE contains no derivatives in t (e.g., the LHS of eq. (5) is zero) it is classified as elliptic. As an example,

u_{xx}+u_{yy}=0 \quad (11)

is Laplace's equation where x and y are spatial independent variables in Cartesian coordinates. Note that with no derivatives in t, elliptic PDEs require no ICs, i.e., they are entirely boundary value PDEs.


PDEs with mixed geometric characteristics are possible, and in fact, are quite common in applications. For example, the PDE

u_{t}=-u_{x}+u_{xx} \quad (12)

is hyperbolic-parabolic. Since it frequently models convection (hyperbolic) through the term ux and diffusion (parabolic) through the term uxx, it is generally termed a convection-diffusion equation. If additionally, it includes a function of the dependent variable u such as

u_{t}=-u_{x}+u_{xx}+f(u) \quad (13)

then it might be termed a convection-diffusion-reaction equation since f(u) typically models the rate of a chemical reaction. If the function is for only the independent variables, i.e.,

u_{t}=-u_{x}+u_{xx}+g(x,t) \quad (14)

the equation could be labeled an inhomogeneous PDE.


This discussion clearly indicates that PDE problems come in an infinite variety, depending, for example, on linearity, types of coefficients (constant, variable), coordinate system, geometric classification (hyperbolic, elliptic, parabolic), number of dependent variables (number of simultaneous PDEs), number of independent variables (number of dimensions), type of BCs, smoothness of the IC, etc., so it might seem impossible to formulate numerical procedures with any generality that can address a relatively broad spectrum of PDEs. But in fact, the MOL provides a surprising degree of generality, although the success in applying it to a new PDE problem depends to some extent on the experience and inventiveness of the analyst, i.e., MOL is not a single, straightforward, clearly defined approach to PDE problems, but rather, is a general concept (or philosophy) that requires specification of details for each new PDE problem. We now proceed to illustrate the formulation of a MOL numerical algorithm with the caveat that this will not be a general discussion of MOL as it is applied to any conceivable PDE problem.

[edit] Elements of the MOL

The basic idea of the MOL is to replace the spatial (boundary value) derivatives in the PDE with algebraic approximations. Once this is done, the spatial derivatives are no longer stated explicitly in terms of the spatial independent variables. Thus, in effect only the initial value variable, typically time in a physical problem, remains. In other words, with only one remaining independent variable, we have a system of ODEs that approximate the original PDE. The challenge, then, is to formulate the approximating system of ODEs. Once this is done, we can apply any integration algorithm for initial value ODEs to compute an approximate numerical solution to the PDE. Thus, one of the salient features of the MOL is the use of existing, and generally well established, numerical methods for ODEs.


To illustrate this procedure, we consider the MOL solution of eq. (8). First we need to replace the spatial derivative ux with an algebraic approximation. In this case we will use a finite difference (FD) such as

u_{x}\approx\frac{u_{i}-u_{i-1}}{\Delta x} \quad (15)


where i is an index designating a position along a grid in x and Δx is the spacing in x along the grid. Thus, for the left end value of \scriptstyle{ x,\; i=1 }, and for the right end value of \scriptstyle{ x,\; i=M }, i.e., the grid in x has M points. Then the MOL approximation of eq. (8) is

\frac{du_{i}}{dt}=-v\frac{u_{i}-u_{i-1}}{\Delta x},1\leq i\leq M \quad (16)

Note that eq. (16) is written as an ODE since there is now only one independent variable, t. Note also that eq. (16) represents a system of M ODEs.


This transformation of a PDE, eq. (8), to a system of ODEs, eqs. (16), illustrates the essence of the MOL, namely, the replacement of the spatial derivatives, in this case ux, so that a system of ODEs approximates the original PDE. Then, to compute the solution of the PDE, we compute a solution to the approximating system of ODEs. But before considering this integration in t, we have to complete the specification of the PDE problem. Since eq. (8) is first order in t and first order in x, it requires one IC and one BC. These will be taken as

u(x,t=0)=f(x) \quad (17a)

.

u(x=0,t)=g(t) \quad (17b)

.

Since eqs. (16) constitute M initial value ODEs, M initial conditions are required and from eq. (17a), these are

u(x_{i},t=0)=f(x_{i}),1\leq i\leq M\quad (18a)

Also, application of BC (17b) gives for grid point i = 1

u(x_{1},t)=g(t) \quad (18b)

Eqs. (16) and (18) now constitute the complete MOL approximation of eq. (8) subject to eqs. (17). The solution of this ODE system gives the M functions

u_{1}(t),u_{2}(t),\ldots u_{M-1}(t),u_{M}(t) \quad (19)

that is, an approximation to u(x,t) at the grid points \scriptstyle{ i=1,2,\ldots M }.


Before we go on to consider the numerical integration of the approximating ODEs, in this case eqs. (16), we briefly consider further the FD approximation of eq. (15), which can be written as

u_{x}\approx\frac{u_{i}-u_{i-1}}{\Delta x}+O(\Delta x)\quad (20)

where O(Δx) denotes of order Δx, that is, the truncation error of the approximation of eq. (16) is proportional to Δx to the first power (varies linearly with Δx); thus eq. (20) is also termed a first order FD (since Δx is to the first power in the order or truncation error term). The term truncation error reflects the fact that the finite difference of eq. (15) comes from a truncated Taylor series.


Note that the numerator of eq. (15), ui − ui − 1, is a difference in two values of u. Also, the denominator Δx remains finite (nonzero). Hence the name finite difference (and it is an approximation because of the truncated Taylor series, so a more complete description is first order finite difference approximation). In fact, in the limit \scriptstyle{ \Delta x\rightarrow0 } the approximation of eq. (15) becomes exactly the derivative. However, in a practical computed-based calculation, Δx remains finite, so eq. (15) remains an approximation.


Also, eq. (8) typically describes the flow of a physical quantity such as concentration of a chemical species or temperature, represented by u, from left to right with respect to x with velocity v. Then, using the FD approximation of eq. (20) at i involves ui and ui − 1. In a flowing system, ui − 1 is to the left (in x) of ui or is upsteam or upwind of ui (to use a nautical analogy). Thus, eq. (20) is termed a first order upwind FD approximation. Generally, for strongly convective systems such as modeled by eq. (8), some form of upwinding is required in the numerical solution of the descriptive PDEs; we will look at this requirement further in the subsequent discussion.

[edit] ODE Integration within the MOL

We now consider briefly the numerical integration of the M ODEs of eqs. (16). If the derivative \scriptstyle{ \frac{du_{i}}{dt} } is approximated by a first order FD

\frac{du_{i}}{dt}\approx\frac{u_{i}^{n+1}-u_{i}^{n}}{\Delta t}+O(\Delta t) \quad (21)

where n is an index for the variable t (t moves forward in steps denoted or indexed by n), then a FD approximation of eq. (16) is

\frac{u_{i}^{n+1}-u_{i}^{n}}{\Delta t}=-v\frac{u_{i}^{n}-u_{i-1}^{n}}{\Delta x}


or solving for \scriptstyle{ u_{i}^{n+1} },

u_{i}^{n+1}=u_{i}^{n}-(v\Delta t/\Delta x)(u_{i}^{n}-u_{i-1}^{n}),i=1,2,\ldots M \quad (22)


Eq. (22) has the important characteristic that it gives \scriptstyle{ u_{i}^{n+1} } explicitly, that is, we can solve for the solution at the advanced point in \scriptstyle{ t,\; n+1 }, from the solution at the base point n. In other words, explicit numerical integration of eqs. (16) is by the forward FD of eq. (21), and this procedure is generally termed the forward Euler method which is the most basic form of ODE integration.


While the explicit form of eq. (22) is computationally convenient, it has a possible limitation. If the time step Δt is above a critcal value, the calculation becomes unstable, which is manifest by successive changes in the dependent variable, \scriptstyle{ \Delta u=u_{i}^{n+1}-u_{i}^{n} }, becoming larger and eventually unbounded as the calculation moves forward in t (for increasing n). In fact, for the solution of eq. (8) by the method of eq. (22) to remain stable, the dimensionless group (vΔt / Δx), which is called the Courant-Friedricks-Lewy or CFL number, must remain below a critical value, in this case, unity. Note that this stability limit places an upper limit on Δt for a given v and Δx; if one attempts to increase the accuracy of the eq. (22) by using a smaller Δx (larger number of grid points in x by increasing M), a smaller value of Δt is required to keep the CFL number below its critical value. Thus, there is a conflicting requirement of improving accuracy while maintaining stability.


The stability limit of the explicit Euler method as implemented via the forward FD of eq. (21) is to use a backward FD for the derivative in t

\frac{du_{i}}{dt}\approx\frac{u_{i}^{n}-u_{i}^{n-1}}{\Delta t}+O(\Delta t) \quad (23)

so that the FD approximation of eqs. (16) becomes

\frac{u_{i}^{n}-u_{i}^{n-1}}{\Delta t}=-v\frac{u_{i}^{n}-u_{i-1}^{n}}{\Delta x}

or after rearrangement (with (vΔt / Δx) = α)

(1+\alpha)u_{i}^{n}+\alpha u_{i-1}^{n}=u_{i}^{n-1},i=1,2,\ldots M \quad (24)

Note that we cannot now solve eq. (24) explicitly for the solution at the advanced point, \scriptstyle{ u_{i}^{n} }, in terms of the solution at the base point \scriptstyle{ u_{i}^{n-1} }. Rather, eq. (24) is implicit in \scriptstyle{ u_{i}^{n} } because \scriptstyle{ u_{i-1}^{n} } is also unknown; that is, we must solve eq. (24) written for each grid point \scriptstyle{ i=1,2,\dots M } as a simultaneous system of bidiagonal equations (bidiagonal because each of eqs. (24) has two unknowns so that simultaneous solution of the full set of approximating algebraic equations is required to obtain the complete numerical solution \scriptstyle{ u_{1}^{n},u_{2}^{n},\dots u_{M}^{n}) }. Thus, the solution of eqs. (24) is termed the implicit Euler method.


We could then naturally ask why use eqs. (24) when eq. (22) is so much easier to use (explicit calculation of the solution at the next step in t of eq. (22) vs the implicit calculation of eqs. (24)). The answer is that the implicit calculation of eqs. (24) is often worthwhile because the implicit Euler method has no stability limit (is unconditionally stable in comparsion with the explicit method with the stability limit stated in terms of the CFL condition). However, there is a price to pay for the improved stability of the implicit Euler method, that is, we must solve a system of simultaneous algebraic equations; eqs. (24) is an example. Furthermore, if the original ODE system approximating the PDE is nonlinear, we have to solve a system of nonlinear algebraic equations (eqs. (24) are linear, so the solution is much easier). The system of nonlinear equations is typically solved by a variant of Newton's method which can become very demanding computationally if the number of ODEs is large (due to the use of a large number of spatial grid points in the MOL approximation of the PDE, especially when we attempt the solution of two dimensional (2D) and three dimensional (3D) PDEs). If you have had some experience with Newton's method, you may appreciate that the Jacobian matrix of the nonlinear algebraic system can become very large and sparse as the number of spatial grid points increases.


Additionally, although there is no limit for Δt with regard to stability for the implicit method, there is a limit with regard to accuracy. In fact, the first order upwind approximation of ux in eq. (8), eq. (20), and the first order approximation of ut in eq. (8), eq. (21) or (23), taken together limit the accuracy of the resulting FD approximation of eq (8). One way around this accuracy limitation is to use higher order FD approximations for the derivatives in eq. (8).


For example, if we consider the second order approximation for uxat i

u_{x}\approx\frac{u_{i+1}-u_{i-1}}{2\Delta x}+O(\Delta x^{2}) \quad (25)

substitution in eq. (8) gives the MOL approximation of eq. (8)

\frac{du_{i}}{dt}=-v\frac{u_{i+1}-u_{i-1}}{2\Delta x},1\leq i\leq M \quad (26)

We could then reason that if the integration in t is performed by the explicit Euler method, i.e., we use the approximation of eq. (21) for \scriptstyle{ u_{t}=\frac{du_{i}}{dt} }, the resulting numerical solution should be more accurate than the solution from eq. (22). In fact, the MOL approximation based on this idea

u_{i}^{n+1}=u_{i}^{n}-\frac{v\Delta t}{2\Delta x}(u_{i+1}^{n}-u_{i-1}^{n}),i=1,2,\ldots M \quad (27)

is unconditionally unstable; this conclusion can be demonstrated by a von Neumann stability analysis that we will not cover here. This surprising result demonstrates that replacing the derivatives in PDEs with higher order approximations does not necessarily guarantee more accurate solutions, or even stable solutions.

[edit] Numerical Diffusion and Oscillation

Even if the implicit Euler method is used for the integration in t of eqs. (26) to achieve stability (or a more sophisticated explicit integrator in t is used that automatically adjusts Δt to achieve a prescribed accuracy), we would find the solution oscillates unrealistically. This numerical distortion is one of two generally observed forms of numerical error. The other numerical distortion is diffusion which would be manifest in the solution from eq. (22). Briefly, the solution would exhibit excessive smoothing or rounding at points in x where the solution changes rapidly. This overall observation that a first order approximation for ux produces numerical diffusion while higher order approximations for ux produce numerical oscillation is predicted by the Godunov order barrier theorem for the Riemann problem [Wesseling, 1]. To briefly explain, the order barrier is first order and any linear FD approximation above first order will be oscillatory. Eq. (8) is an example of the Riemann problem [Wesseling, 1] if IC eq. (17a) is discontinuous, for example u(x,t = 0) = h(t) where h(t) is the Heaviside unit step function. The (exact) analytical solution is the initial condition function f(x) of eq. (17a) moving left to right with velocity v (from eq. (8)) and without distortion, i.e., u(x,t) = f(x − vt); however, the numerical solution will oscillate if ux in eq. (8) is replaced with a linear approximation of second or higher order.


We should also mention one point of terminology for FD approximations. The RHS of eq. (25) is an example of a centered approximation since the two points at i + 1 and i − 1 are centered around the point i. Eq. (20) is an example of a noncentered, one-sided or upwind approximation since the points i and i − 1 are not centered with respect to i.


Finally, to conclude the discussion of first order hyperbolic PDEs such as eq. (8), since the Godunov theorem indicates that FD approximations above first order will produce numerical oscillations in the solution, the question remains if there are approximations above first order that are nonoscillatory. To answer this question we note first that the Godunov theorem applies to linear approximations; for example, eq. (25) is a linear approximation since u on the RHS is to the first power. If, however, we consider nonlinear approximations for ux, we can in fact develop approximations that are nonoscillatory. The details of such nonlinear approximations are beyond the scope of this discussion, so we will merely mention that they are termed high resolution methods which seek a total variation diminishing (TVD) solution. Such methods, which include flux limiter and weighted essentially nonoscillatory (WENO) methods, seek to avoid non-real oscillations when shocks or discontinuities occur in the solution (such as in the Riemann problem) [Shu, 2].


So far we have considered only the MOL solution of first order PDEs, e.g., eq. (8). We conclude this discussion of the MOL by considering a second order PDE, the parabolic eq. (1). To begin, we need an approximation for the second derivative uxx. A commonly used second order, central approximation is (again, derived from the Taylor series, so the term O(Δx2) represents the truncation error)

u_{xx}\approx\frac{u_{i+1}-2u_{i}+u_{i-1}}{\Delta x^{2}}+O(\Delta x^{2}) \quad (28)

Substitution of eq. (28) in eq. (8) gives a system of approximating ODEs

\frac{du_{i}}{dt}=D\frac{u_{i+1}-2u_{i}+u_{i-1}}{\Delta x^{2}},i=1,2,\ldots M \quad (29)


Eqs. (29) are then integrated subject to IC (2) and BCs (3). This integration in t can be by the explicit Euler method, the implicit Euler method, or any other higher order integrator for initial value ODEs. Generally stability is not as much of a concern as with the previous hyperbolic PDEs (a characteristic of parabolic PDEs which tend to smooth solutions rather than hyperbolic PDEs which tend to propagate nonsmooth conditions). However, stability constraints do exist for explicit methods. For example, for the explicit Euler method with a step Δt in t, the stability constraint is DΔt / Δx2 < constant (so that as Δx is reduced to achieve better spatial accuracy in x, Δt must also be reduced to maintain stability).


Before proceeding with the integration of eqs. (29), we must include BCs (3). The Dirichlet BC at x = x0, eq. (3a), is merely

u_{1}=u_{b} \quad (30)

and therefore the ODE of eqs. (29) for i = 1 is not required and the ODE for i = 2 becomes

\frac{du_{2}}{dt}=D\frac{u_{3}-2u_{2}+u_{b}}{\Delta x^{2}} \quad (31)

[edit] Differential Algebraic Equations

Eq. (30) is algebraic, and therefore in combination with the ODEs of eqs. (29), we have a differential algebraic (DAE) system.


At i = M, we have eqs. (29)

\frac{du_{M}}{dt}=D\frac{u_{M+1}-2u_{M}+u_{M-1}}{\Delta x^{2}}\quad (32)

Note that uM + 1 is outside the grid in x, i.e., M + 1 is a fictitious point. But we must assign a value to uM + 1 if eq. (32) is to be integrated. Since this requirement occurs at the boundary point i = M, we obtain this value by approximating BC (3b) using the centered FD approximation of eq. (25)

u_{x}\approx\frac{u_{M+1}-u_{M-1}}{2\Delta x}=0

or

u_{M+1}=u_{M-1} \quad (33)

We can add eq. (33) as an algebraic equation to our system of equations, i.e., continue to use the DAE format, or we can substitute uM + 1 from eq. (33) into the eq. (32)

\frac{du_{M}}{dt}=D\frac{u_{M-1}-2u_{M}+u_{M-1}}{\Delta x^{2}} \quad (34)

and arrive at an ODE system (eqs. (29) for \scriptstyle{ i=3,\ldots M-1 }, eq. (31) for i = 2 and eq. (34) for i = M). Both approaches, either an ODE system or a DAE system, have been used in MOL studies. Either way, we now have a complete formulation of the MOL ODE or DAE system, including the BCs at i = 1,M in eqs. (29). The integration of these equations then gives the numerical solution \scriptstyle{ u_{1}(t),u_{2}(t),\ldots u_{M}(t) }. The preceding discussion is based on a relatively basic DAE system, but it indicates that integrators designed for DAE systems can play an important role in MOL analysis.


If the implicit Euler method is applied to eqs. (29), we have

\frac{u_{i}^{n+1}-u_{i}^{n}}{\Delta t}=D\frac{u_{i+1}^{n+1}-2u_{i}^{n+1}+u_{i-1}^{n+1}}{\Delta x^{2}},i=1,2,\ldots M

or (with α = DΔt / Δx2)

u_{i+1}^{n+1}-(1/\alpha+2)u_{i}^{n+1}+u_{i-1}^{n+1}=(1/\alpha)u_{i}^{n},i=1,2,\ldots M

which is a tridiagonal system of algebraic equations (three unknowns in each equation). Since such banded systems (the nonzero elements are banded around the main diagonal) are common in the numerical solution of PDE systems, special algorithms have been developed to take advantage of the banded structure, typically by not storing and using the zero elements outside the band. These special algorithms that take advantage of the structure of the problem equations can result in major savings in computation time. In the case of tridiagonal equations, the special algorithm is generally called Thomas' method. If the coefficient matrix of the algebraic system does not have a well defined structure, such as bidiagonal or tridiagonal, but consists of mostly zeros with a relatively small number of nonzero elements, which is often the case in the numerical solution of PDEs, the coefficient matrix is said to be sparse; special algorithms and associated software for sparse systems have been developed that can result in very substantial savings in the storage and numerical manipulation of sparse matrices.


Generally when applying the MOL, the integration of the approximating ODE/DAEs (e.g., eqs. (21) and (29)) is accomplished by using library routines for initial value ODE/DAEs. In other words, the explicit programming of the ODE/DAE integration (such as the explicit or implicit Euler method) is avoided; rather, an established integrator is used. This has the advantage that: (1) the detailed programming of the integration can be avoided, particularly the linear algebra (solution of simultaneous equations) required by an implicit integrator, so that the MOL analysis is substantially simplified, and (2) library routines (usually written by experts) include features that make these routines especially effective (robust) and efficient such as automatic integration step size adjustment and the use of higher order integration methods (beyond the first order accuracy of the Euler methods); also, generally they have been thoroughly tested. Thus, the use of quality ODE/DAE library routines is usually an essential part of MOL analysis. We therefore list at the end of this article some public domain sources of such library routines.

[edit] Higher Dimensions and Different Coordinate Systems

To conclude this discussion of the MOL solution of PDEs, we cover two additional points. First, we have considered PDEs in only Cartesian coordinates, and in fact, just one Cartesian coordinate, x. But MOL analysis can in principle be carried out in any coordinate system. Thus, eq. (1) can be generalized to

\dfrac{\partial u}{\partial t}=D\nabla^{2}u \quad (35)

where \scriptstyle{ \nabla^{2} } is the coordinate independent Laplacian operator which can then be expressed in terms of a particular coordinate system. For example, in cylindrical coordinates eq. (35) is

\frac{\partial u}{\partial t}=D\left(\frac{\partial^{2}u}{\partial r^{2}}+\frac{1}{r}\frac{\partial u}{\partial r}+\frac{1}{r^{2}}\frac{\partial^{2}u}{\partial\theta^{2}}+\frac{\partial^{2}u}{\partial z^{2}}\right) \quad (36)


and in spherical coordinates eq. (35) is


\frac{\partial u}{\partial t}=D\left[\frac{\partial^{2}u}{\partial r^{2}}+\frac{2}{r}\frac{\partial u}{\partial r}+\frac{1}{r^{2}}\left(\frac{\partial^{2}u}{\partial\theta^{2}}+\frac{\cos\theta}{\sin\theta}\frac{\partial u}{\partial\theta}\right)+\frac{1}{r^{2}\sin^{2}\theta}\frac{\partial^{2}u}{\partial\phi^{2}}\right] \quad (37)

.

The challenge then in applying the MOL to PDEs such as eqs. (36) and (37) is the algebraic approximation of the RHS (\scriptstyle{ \nabla^{2}u }) using for example, FDs, finite elements or finite volumes; all of these approximations have been used in MOL analysis, as well as Galerkin, least squares, spectral and other methods. A particularly demanding step is regularization of singularities such as at r=0 (note the number of divisions by r in the RHS of eqs. (36) and (37)) and at θ = 0,π / 2 (note the divisions by sin(θ) in eq. (37)). Thus the application of the MOL typically requires analysis based on the experience and creativity of the analyst (i.e., it is generally not a mechanical procedure from beginning to end).


The complexity of the numerical solution of higher dimensional PDEs in various coordinate systems prompts the question of why a particular coordinate system would be selected over others. The mathematical answer is that the judicious choice of a coordinate system facilitates the implementation of the BCs in the numerical solution. The answer based on physical considerations is that the coordinate system is selected to reflect the geometry of the problem system. For example, if the physical system has the shape of a cylinder, cylindrical coordinates would be used. This choice then facilitates the implementation of the BC at the exterior surface of the physical system (exterior surface of the cylinder). However, this can also lead to complications such as the r = 0 singularities in eq. (36) (due to the variable 1 / r and 1 / r2 coefficients). The resolution of these complications is generally worth the effort rather than the use of a coordinate system that does not naturally conform to the geometry of the physical system. If the physical system is not shaped in accordance with a particular coordinate system, i.e., has an irregular geometry, then an approximation to the physical geometry is used, generally termed body fitted coordinates.


As a concluding point, we might consider the origin of the name "method of lines". If we consider eqs. (29), integration of this ODE system produces the solution \scriptstyle{ u_{2}(t),u_{3}(t),\ldots u_{M}(t) } (note u1(t) = ub, a constant, from BC (3a)). We could then plot these functions in an x − u(x,t) plane as a vertical line at each \scriptstyle{ x \; (i=2,3,\dots M) } with the height of each line equal to u(xi,t). In other words, the plot of the solution would be a set of vertical parallel lines suggesting the name "method of lines" [Schiesser, 3].

[edit] Sources of ODE/DAE Integrators

One of the very useful aspects of MOL is that it enables tried and tested ODE/DAE numerical routines to be used, many of which are in the public domain. The following sources are a good starting point for these routines. For example, the LSODE and VODE series of ODE/DAE integrators [2s], DASSL for DAEs [2s] and the SUNDIALS library [5s] are widely used in MOL analysis; test problems for ODE/DAE routines are also available [6s]. Additionally, routines that can be called from MOL codes are available to perform a variety of complementary computations (e.g., functional approximation by interpolation, evaluation of integrals, maximization and minimization in optimization associated with the solution of PDEs) [1s,3s,4s]


[1s] www.netlib.org/
[2s] www.netlib.org/ode/index.html (emphasis on ODE/DAE software that can be used in MOL analysis)
[3s] gams.nist.gov/
[4s] www.acm.org/toms/
[5s] www.llnl.gov/CASC/software.html
[6s] www.dm.uniba.it/~sim testset/

[edit] References

[1] Wesseling, P. (2001), Principles of Computational Fluid Dynamics, Springer, Berlin

[2] Shu, C-W (1998), Essentially Non-oscillatory and Weighted Essential Non-oscilliatory Schemes for Hyperbolic Conservation Laws. In: Cockburn, B., Johnson, C., Shu, C-W., Tadmor, E. (Eds.), Advanced Numerical Approximation of Nonlinear Hyperbolic Equations, Lecture Notes in Mathematics, vol 1697. Springer, pp325-432.

[3] Schiesser, W. E. (1991), Numerical Method of Lines Integration of Partial Differential Equations, Academic Press, San Diego; see also http://eqworld.ipmnet.ru, http://www.lehigh.edu/~wes1/apci/28apr00.pdf



S. Hamdi (Canada) 
samir.hamdi@utoronto.ca
G. W. Griffiths (UK) 
graham@griffiths1.com
W. E. Schiesser (USA) 
wes1@lehigh.edu