San José State University
Thayer Watkins
Silicon Valley
& Tornado Alley

The Relationship of Solutions
to Partial Differential Equations and
their Corresponding Finite Difference Schemes

Partial Differential Equations in Operator Format

A typical partial differential equation (PDE)involved in physical processes is of the form:

Ltq = Lx,y,zq

where q is some dependent variable and Lt is a polynomial in the partial derivative operator with respect to time; i.e.; ∂/∂t, which can also be represented more succinctly as ∂t. Likewise Lx,y,z is a polynomial in the partial derivative operators with respect to the spatial variables x, y and z; i.e.; ∂x, ∂y, and ∂z. For example, the heat equation

∂T/∂t = K∇2T

is, in the operator format,

(∂t)T = K(∂x2+∂y2+∂z2)T

The above format presumes that the PDE is homogeneous; i.e.; that each term involves the dependent variable so there is no constant term or forcing term.

Finite Difference Schemes in Operator Format

The finite difference approximation of the partial derivatives can be represented in terms of the forward difference operator Δ. For example, the forward difference approximation of the first derivative is:

∂q/∂x = (qi+1 - qi)/h

where h is the gridlength Δx. This can be represented in terms of the operator Uxq(x)=q(x+h) as

∂q/∂x = ((Ux-1)qi)/h

In terms of this representation the second derivative takes the form

2q/∂x2 = ((Ux2-2Ux+1)qi)/h2

The backward difference approximation to the first derivative has the operator representation

∂q/∂x = ((1 - Ux-1)qi)/h

Each of the above approximations are defective in that they represent a better approximation of the derivatives at points other than at x. In the case of the second derivative the above formula would better represent the derivative at x+h than at x. To get the derivatives at x the formulas should be:

∂q/∂x = ((Ux-Ux-1)qi)/2h
2q/∂x2 = ((Ux -2 + Ux-1)qi)/h2

When such finite differencing schemes are substituted for the derivatives in a PDE the result is a difference equation. The resulting difference equation is typically of a higher order than the PDE so there cannot be a one-to-one correspondence between the independent solutions of the PDE and the independent solutions to the difference equation. The extra solution types of the difference equation compared to the differential equation are usually called computational modes, to contrast them with the ones which correspond to solution types of the differential equation, which are called physical modes.

It is common practice to use different approximation schemes for the different order derivatives in an equation; e.g., a forward difference for the first derivative and a centered difference for the second derivative.

In any case, the difference equation that results from the substitution of the approximations for the derivatives can be considered as an equation in terms of the the forward operator U and the gridlengths or the time step.

Methods of Finding Analytical Solutions
to Partial Differential Equations

While there are existence theorems that establish that there exitst solutions to a given PDE, there is no reason to expect that a given PDE will have solutions which can be expressed in terms of the arbitrary simple functions which have been historically defined. Therefore there is no general method for finding such analytical solutions or even establishing whether or not they exist. However, for the important PDE's involved in physical processes, there are some techniques that are fruitful such as separation of variables or the method of characteristics.

Separation of variables involves assuming the solution to the PDE is a product of functions of the separate independent variables. For example, for the heat equation in one dimension

∂T/∂t = K2(∂2T/∂x2

it is assumed that the temperature T(x,t) = R(x)S(t).

This results in the equation

R(x)dS/dt = K2S(t)d2R/dx2
which after division by RS is
(1/S)dS/dt = K2(1/R)d2R/dx2

The LHS is a function of t alone and the RHS is a function of x alone. This means the common value of the LHS and RHS must be a constant, called the separation constant. For the heat equation in order to satisfy the boundary condition that the temperature approach a finite limit as x goes to infinity the separation constant must be a negative number. This requirement is achieved by making the separation constant equal to -λ2. This reduces the solution of the heat equation to the solution of two ordinary differential equations

dS/dt = -λ2S
d2R/dx2 = -(λ/K)2R

The ordinary differential equations which arise under separation of variables typically do have analytical solutions and these solutions can be found by assuming a solution of the form eyt and the conditions that have to be satisfied by y. This involves finding the roots of a polynomial equation. Likewise the associated difference equation has analytical solutions of the form xi and finding the conditions to be satisfied by x. These values of x can be found by finding the roots of a polynomial equation. Note that because of the different forms of assumed solution the relationship between parameters y and x for the solutions is that eyh=x and hence y=ln(x)/h, where h is the gridlength or timestep.

The polynomials involves in finding the solutions to the differential equations and the difference equation have an interesting relationship. Suppose we have a polynomial equation in y

anyn+ an-1yn-1+... a2y2+ a1y+ a0 = 0.

and we have a polynomial relationship between y and x

y = bmxm+ bm-1xm-1+ ... b2x2+ b1x+ b0

When this polynomial is substituted into the polynomial for y the result is a polynomial in x, of degree nm. The interesting question is the relationship of the solutions to the polynomial in x to the solutions to the polynomial in y.

The polynomial in y can be factored to give:

(y - α1)(y - α2)×...× (y - αn) = 0.

The roots of the polynomial; i.e.; the values of y that make the value of the polynomial zero; are then α1, α2,...,αn. When the polynomial relationship between y and x is substituted into the polynomial in y the result involves terms of the form:

(bmxm+ bm-1xm-1+ ... b2x2+ b1x+ b0 - αi)

Such terms can be factored to

(x - γi1)(x - γi2)...(x - γim)

where γij are the roots of the polynomial

bmxm+ bm-1xm-1+ ... b2x2+ b1x+ b0 - αi = 0.

These are the roots of the polynomial in x. If αi is zero then the γ values would be the same as the roots of the polynomial

bmxm+ bm-1xm-1+ ... b2x2+ b1x+ b0 = 0.

which could be called the β values. The γ values can thus be considered perturbations of the β values based upon the α values, the roots of the polynomial equation in y. For the relevant finite difference approximation the β values are all equal to unity.

Consider the heat equation once again. After separation of variables the equations to be solved are:

dS/dt = -λ2S
d2R/dx2 = -(λ/K)2R

These are to be solved separately and lead to the trivally simple polynomials:

q + λ2 = 0
p2 + (λ/K)2 = 0

The solutions are of course

q = - λ2
p = ±i(λ/K)

where i is the square root of -1.

Now consider a forward finite difference approximation for the time derivative (U - 1)S/h, where h is the time step. This leads to the polynomial

s - 1 + hλ2 = 0
s = 1 - hλ2

Note that the relationship between a root y for the solution for the differential equation and a root x for the solution for the difference equation is y=ln(x)/h. In this case the relation is

ln(s)/h = ln(1±hλ2)/h
which for small h reduces to
ln(s)/h = ±λ2

For the second derivative with respect to x the finite difference approximation might be (U - 2 + U-1)R/k2 where k is the gridlength. The polynomials to be solved for the spatial solution are:

r - 2 + r-1 + k2i(λ/K) = 0
r - 2 + r-1 - k2i(λ/K) = 0

These reduce to the quadratic equations

r2 - 2r + 1 + rk2i(λ/K) = 0
r2 - 2r + 1 - rk2i(λ/K) = 0

which have the solutions

r = 1 + k2i(λ/K)/2 ±((1+k2i(λ/K)/2)2 -1)1/2
r = 1 - k2i(λ/K)/2 ±((1-k2i(λ/K)/2)2 -1)1/2

which reduce to

r = 1 + ik2(λ/K)/2 ±(-k4(λ/K)2/4 + ik2(λ/K))1/2
r = 1 - ik2(λ/K)/2 ±(-k4(λ/K)2/4 - ik2(λ/K))1/2

An ik2(λ/K)/2 term may be factored from the square root expression to give the roots in the form

r = 1 - ik2(λ/K)/2(1 ± (1 - 2i/k2(λ/K))1/2)
r = 1 + ik2(λ/K)/2(1 ± (1 + 2i/k2(λ/K))1/2)

There are thus four roots; two pairs in which each pair is a perturbation about some values. The magnitude of the deviations from these roots depends upon the gridlength k. As k goes to zero both roots of a pair go to the common value.

HOME PAGE OF applet-magic
HOME PAGE OF Thayer Watkins