Partial differential equations

NUMERICAL METHODS
These can then be treated in the way indicated in the previous paragraph. The
extension to more than one dependent variable is straightforward.
In practical problems it often happens that boundary conditions applicable
to a higher-order equation consist not of the values of the function and all its
derivatives at one particular point but of, say, the values of the function at two
separate end-points. In these cases a solution cannot be found using an explicit
step-by-step ‘marching’ scheme, in which the solutions at successive values of the
independent variable are calculated using solution values previously found. Other
methods have to be tried.
One obvious method is to treat the problem as a ‘marching one’, but to use a
number of (intelligently guessed) initial values for the derivatives at the starting
point. The aim is then to find, by interpolation or some other form of iteration,
those starting values for the derivatives that will produce the given value of the
function at the finishing point.
In some cases the problem can be reduced by a differencing scheme to a matrix
equation. Such a case is that of a second-order equation for y(x) with constant
coefficients and given values of y at the two end-points. Consider the second-order
equation
y + 2ky + µy = f(x),
(27.85)
with the boundary conditions
y(0) = A,
y(1) = B.
If (27.85) is replaced by a central difference equation,
yi+1 − 2yi + yi−1
yi+1 − yi−1
+ µyi = f(xi ),
+ 2k
h2
2h
we obtain from it the recurrence relation
(1 + kh)yi+1 + (µh2 − 2)yi + (1 − kh)yi−1 = h2 f(xi ).
For h = 1/(N − 1) this is in exactly the form of the N × N tridiagonal matrix
equation (27.30), with
b1 = bN = 1,
ai = 1 − kh,
bi = µh2 − 2,
c1 = aN = 0,
ci = 1 + kh,
i = 2, 3, . . . , N − 1,
and y1 replaced by A, yN by B and yi by h2 f(xi ) for i = 2, 3, . . . , N − 1. The
solutions can be obtained as in (27.31) and (27.32).
27.8 Partial differential equations
The extension of previous methods to partial differential equations, thus involving
two or more independent variables, proceeds in a more or less obvious way. Rather
1030
27.8 PARTIAL DIFFERENTIAL EQUATIONS
than an interval divided into equal steps by the points at which solutions to the
equations are to be found, a mesh of points in two or more dimensions has to be
set up and all the variables given an increased number of subscripts.
Considerations of the stability, accuracy and feasibility of particular calculational schemes are the same as for the one-dimensional case in principle, but in
practice are too complicated to be discussed here.
Rather than note generalities that we are unable to pursue in any quantitative
way, we will conclude this chapter by indicating in outline how two familiar
partial differential equations of physical science can be set up for numerical
solution. The first of these is Laplace’s equation in two dimensions,
∂2 φ ∂2 φ
+ 2 = 0,
∂x2
∂y
(27.86)
the value of φ being given on the perimeter of a closed domain.
A grid with spacings ∆x and ∆y in the two directions is first chosen, so that,
for example, xi stands for the point x0 + i∆x and φi,j for the value φ(xi , yj ). Next,
using a second central difference formula, (27.86) is turned into
φi+1,j − 2φi,j + φi−1,j
φi,j+1 − 2φi,j + φi,j−1
+
= 0,
(∆x)2
(∆y)2
2
(27.87)
2
for i = 0, 1, . . . , N and j = 0, 1, . . . , M. If (∆x) = λ(∆y) then this becomes the
recurrence relationship
φi+1,j + φi−1,j + λ(φi,j+1 + φi,j−1 ) = 2(1 + λ)φi,j .
(27.88)
The boundary conditions in their simplest form (i.e. for a rectangular domain)
mean that
φ0,j ,
φN,j ,
φi,0 ,
φi,M
(27.89)
have predetermined values. Non-rectangular boundaries can be accommodated,
either by more complex boundary-value prescriptions or by using non-Cartesian
coordinates.
To find a set of values satisfying (27.88), an initial guess of a complete set of
values for the φi,j is made, subject to the requirement that the quantities listed
in (27.89) have the given fixed values; those values that are not on the boundary
are then adjusted iteratively in order to try to bring about condition (27.88)
everywhere. Clearly one scheme is to set λ = 1 and recalculate each φi,j as the
mean of the four current values at neighbouring grid-points, using (27.88) directly,
and then to iterate this recalculation until no value of φ changes significantly
after a complete cycle through all values of i and j. This procedure is the simplest
of such ‘relaxation’ methods; for a slightly more sophisticated scheme see exercise
27.26 at the end of this chapter. The reader is referred to specialist books for
fuller accounts of how this approach can be made faster and more accurate.
1031
NUMERICAL METHODS
Our final example is based upon the one-dimensional diffusion equation for
the temperature φ of a system:
∂2 φ
∂φ
=κ 2.
∂t
∂x
(27.90)
If φi,j stands for φ(x0 + i∆x, t0 + j∆t) then a forward difference representation
of the time derivative and a central difference representation for the spatial
derivative lead to the following relationship:
φi+1,j − 2φi,j + φi−1,j
φi,j+1 − φi,j
=κ
.
∆t
(∆x)2
(27.91)
This allows the construction of an explicit scheme for generating the temperature
distribution at later times, given that it is known at some earlier time:
φi,j+1 = α(φi+1,j + φi−1,j ) + (1 − 2α)φi,j ,
(27.92)
where α = κ∆t/(∆x)2 .
Although this scheme is explicit, it is not a good one because of the asymmetric
way in which the differences are formed. However, the effect of this can be
minimised if we study and correct for the errors introduced in the following way.
Taylor’s series for the time variable gives
φi,j+1 = φi,j + ∆t
(∆t)2 ∂2 φi,j
∂φi,j
+
+ ··· ,
∂t
2! ∂t2
(27.93)
using the same notation as previously. Thus the first correction term to the LHS
of (27.91) is
−
∆t ∂2 φi,j
.
2 ∂t2
(27.94)
The first term omitted on the RHS of the same equation is, by a similar argument,
−κ
2(∆x)2 ∂4 φi,j
.
4!
∂x4
(27.95)
But, using the fact that φ satisfies (27.90), we obtain
2 ∂4 φ
∂2 φ
∂
∂2 ∂φ
∂ φ
=
=
κ
κ
= κ2 4 ,
2
2
2
∂t
∂t
∂x
∂x
∂t
∂x
(27.96)
and so, to this accuracy, the two errors (27.94) and (27.95) can be made to cancel
if α is chosen such that
−
2κ(∆x)2
κ2 ∆t
=−
,
2
4!
1032
i.e. α =
1
.
6