2. Classification of 2nd-order Quasi-linear PDEs in Two Variables

Most physical systems are governed by second order PDEs. In this course we discuss FDMs for solving such equations. The most reliable algorithms are those which model the physics and so to begin with, we must understand the physical background. We consider the equation for u(x, y)

au_{xx} + bu_{xy} + cu_{yy} = f .  (2.1)

This equation is called quasi-linear provided the functions a, b, c and f do not depend on u_{xx}, u_{xy} or u_{yy}. They may, however depend on x, y, u, u_{x} and u_{y}, so that (2.1) is not necessarily linear.

Suppose we know u, u_{x} and u_{y} along some curve Γ in (x, y)-space. From a point P on Γ we move a small vector displacement (dx, dy) to a new point Q not on Γ. Under what circumstances can we determine uniquely the values of u, u_{x} and u_{y} at Q? We denote the change in these variables by du, d(u_{x}) and d(u_{y}). Then by the chain rule for partial derivatives, du = u_{x}dx + u_{y}dy which is known because u_{x} and u_{y} are known along Γ. Similarly,

\begin{matrix}d(u_{x}) & = u_{xx}dx + u_{xy}dy \\ d(u_{y}) & = u_{xy}dx + u_{yy}dy \end{matrix} \} .  (2.2)

We combine (2.1) and (2.2) in matrix form:

( \begin{matrix}a & b & c \\ dx & dy & 0 \\ 0 & dx & dy \end{matrix} ) ( \begin{matrix}u_{xx} \\ u_{xy} \\ u_{yy} \end{matrix} ) = ( \begin{matrix}f \\ d(u_{x}) \\ d(u_{y}) \end{matrix} ) .  (2.3)

a, b and c are known locally because u, u_{x} and u_{y} are known, and so the 3 × 3 matrix is known. Equation (2.3) will have a unique solution for u_{xx}, u_{xy} and u_{yy} unless the determinant of that matrix vanishes, that is unless

| \begin{matrix}a & b & c \\ dx & dy & 0 \\ 0 & dx & dy \end{matrix} | = a(dy)^{2} − b dxdy + c(dx)^{2} = 0 .  (2.4)

If (2.4) holds, the equation (2.3) will have either no solution or infinitely many. The condition for solutions to exist, which we will use later in the course, is that

a \frac{d(u_{x})}{dx} + c \frac{d(u_{y})}{dy} = f .  (2.5)

This is surprising. If we can choose a direction (dx, dy) which satisfies (2.4) we have the possibility that the second derivatives u_{xx} etc. may not be uniquely defined. In other words, the solution may have discontinuities across the line P Q, with u_{xx} taking different values on each side.

It is very important to know whether our solution can have this property. Equation (2.5) is called the Characteristic equation of (2.1). It is a quadratic in dy/dx with solution

\frac{dy}{dx} = \frac{b ± \sqrt{b^{2} − 4ac}}{2a} .  (2.6)

Equation (2.1) is classified as hyperbolic, parabolic or elliptic according to whether these roots are real. Note the sign of b in the formula.

If \{ \begin{matrix}b^{2} − 4ac & > 0  2 real roots, (2.1) is hyperbolic \\ b^{2} − 4ac & = 0  1 real root, (2.1) is parabolic \\ b^{2} − 4ac & < 0  0 real roots, (2.1) is elliptic. \end{matrix} \} .  (2.7)

For hyperbolic equations, (2.6) is an ODE for y(x) which can be integrated to define two sets of curves (one for the + sign, one for the −), called the characteristics of (2.1).

Example: Significance of Characteristics.

Consider the one-dimensional wave equation for u(x, t) with constant c

u_{xx} = \frac{1}{c^{2}} u_{tt}  for t > 0  with u(x, 0) = F (x)  and u_{t}(x, 0) = G(x) . (2.8)

This equation has the general solution u = f(x − ct) + g(x + ct) for any functions f and g and with the above boundary conditions we have d’Alembert’s solution:

u(x, t) = \frac{1}{2}[F (x + ct) + F (x − ct)] + \frac{1}{2c} ∫_{x−ct}^{x+ct} G(ξ)dξ .  (2.9)

Using (2.6) the characteristics of (2.8) are two families of straight lines

\frac{dt}{dx} = ±\frac{1}{c}  or x ± ct = constant.  (2.10)

From the actual solution, we see that the solution at some point P or (x_{0}, t_{0}) with t_{0} > 0 depends only on some of the initial data, that for which x_{0}−ct_{0} ⩽ x ⩽ x_{0}+ct_{0}. Only points from which characteristics going forwards in time can reach the point P can influence the solution at P. The set of such points is called the domain of dependence of P. In exactly the same way not all points with t > t_{0} can be affected by the solution at P. The collection of such points is called the domain of influence of P.

This behaviour is easy to understand physically, if one interprets characteristics as curves along which information travels at a finite speed. If something happens at P it takes a certain time before news of it reaches another point. For (2.8), the characteristics are parallel lines, but for more general hyperbolic systems they will be curved and may meet. In such cases discontinuities may form. If neighbouring characteristics touch, conflicting information arrives at the same point, leading to the creation of shock waves (such as sonic booms, or pressure fronts.) Such discontinuities then propagate along the characteristics.

It is clear from this example that whether or not characteristics exist is vital for the understanding and therefore the numerical modelling of a problem. They are associated with “time-like” behaviour, and have a characteristic speed associated with them, defining the rate at which information travels. In contrast elliptic problems have no “time-like” variable; x and y behave like space coordinates.

Boundary Conditions for Well-Posed Problems.

A problem involving a PDE is said to be ‘well-posed’ if three conditions hold:

(1) A solution exists;

(2) The solution is unique;

(3) The solution is continuous in the boundary conditions, i.e. small changes in the boundary conditions do not lead to large changes in the solution. If this last condition fails to hold, the problem is non-physical and a disaster for numerical modelling.

Whether or not a problem is well-posed depends critically on whether its boundary conditions are appropriate. Typical boundary conditions are:

(a) boundary value problems (BVP); the PDE holds in some closed region and the solution is constrained all over the boundary;

(b) initial value problems (IVP); One or more constraints are given on some curve (usually t = 0) only partially bounding the region in which the PDE holds.

Hyperbolic Equations

From our discussion of characteristics, it is clear that hyperbolic systems should have initial value conditions. Information spreads out from the initial values at a finite rate.

An example of an ill-posed hyperbolic BVP for u(x, y) with a non-unique solution is

u_{xx} = u_{yy}  in \{ \begin{matrix}0 & < x < 1 \\ 0 & < y < 1 \end{matrix} \}  with \{ \begin{matrix}u(x, 0) = u(x, 1) & = 0 \\ u(0, y) = u(1, y) & = 0 \end{matrix} \} .  (2.11)

This problem has the solution u = A sin nπx sin nπy for any constant A and integer n.

Elliptic Equations

These have no characteristics; no lines along which information travels, which suggests that IVPs are inappropriate. A typical elliptic equation is Laplace’s equation

∇^{2}u ≡ u_{xx} + u_{yy} = 0  in D,  (2.12)

where D is some region of (x, y)-space. It can be shown that this equation together with one boundary condition (say u = f) on the boundary ∂D can give a well-posed problem. In contrast, the IVP

u_{xx} + u_{yy} = 0  in y > 0  with u(x, 0) = u_{y}(x, 0) = 0  (2.13)

is ill-posed, even though the solution u(x, y) = 0 is unique! Suppose we perturb the initial conditions so that u(x, 0) = ε sin nx, and u_{y}(x, 0) = 0, where 0 < ε ≪ 1 and n is arbitrary but large. No matter how big n is, | sin nx| ⩽ 1, so we are not altering the boundary condition by more than ε. The (unique) solution to this new problem is

u = ε sin nx cosh ny ≃ ε sin nx \frac{1}{2}e^{|ny|}  when |ny| ≫ 1 .  (2.14)

So a small distance away from the initial line y = 0, the solution is now exponentially large, whereas for the unperturbed problem it was zero. If a tiny (albeit very wiggly) perturbation to the boundary conditions can lead to a vast difference in the solution the problem is physically meaningless and impossible to model numerically. This example, due to Hadamard, shows that IVPs for elliptic equations are discontinuous in the boundary conditions.

Parabolic Equations A typical example is the diffusion equation for u(x, t) with constant diffusivity K:

u_{t} = Ku_{xx}  with u(x, 0) = f(x) .  (2.15)

As we know u(x, 0), we can calculate u_{t}(x, 0) = Kf′′(x) from the equation, and so we would expect to be able to step away from t = 0. From (2.6), the characteristics are given by the repeated root dt/dx = 0 or t = constant. This corresponds to an infinite speed of propagation of information. Is the solution stable? Once more we consider the boundary condition f(x) = ε sin nx, so that the unique solution of (2.15) is

u(x, t) = ε sin nx e^{−n^{2}Kt} .  (2.16)

When ε = 0 the solution is u = 0. When ε > 0 the solution decays away provided Kt > 0, but if Kt < 0 and n is large, the perturbed solution blows up once more.

Thus, parabolic equations require one initial condition and it is vital that we move “forwards in time.” Physically, parabolic equations describe the smoothing out of an initial configuration towards an equilibrium. Many different initial conditions give rise to almost the same final state. This is why running the process backwards in time is an ill-posed problem. You can’t un-stir a cup of tea!

Summary

Equation type  Appropriate B.C.  Method of solution
Hyperbolic  2 Initial  Step in either direction from initial line
Parabolic  1 Initial  Step in one direction only from initial line
Elliptic  1 Boundary  Must solve everywhere simultaneously

Examples of various equation types

Hyperbolic  Parabolic  Elliptic
(characteristic speed c)  (diffusivity K)
Maxwell’s equations  Heat equation  Electrostatics (Poisson)
Unsteady 1-D compressible  Unsteady incompressible N-S Potential flow
Steady 2-D supersonic  Steady boundary layer  Steady Navier-Stokes

Clearly, an understanding of each type of equation is essential. We shall consider Finite Difference schemes appropriate to each in turn, reflecting the physics, where possible.