Partial differential equations

12.1
Taylor Series
In this chapter we discuss finite difference approximations to partial derivatives. The ap-
proximations are based on Taylor series expansions of a function of one or more variables.
Recall that the Taylor series expansion for a function of one variable is given by
f (x + h) = f (x) +
h
1!
f

(x) +
h
2
2!
f

(x) +
· · ·
(12.1.1)
The remainder is given by
f
(n)
(ξ)
h
n
n!
,
ξ(x, x + h).
(12.1.2)
For a function of more than one independent variable we have the derivatives replaced by
partial derivatives. We give here the case of 2 independent variables
f (x + h, y + k)
=
f (x, y) +
h
1!
f
x
(x, y) +
k
1!
f
y
(x, y) +
h
2
2!
f
xx
(x, y)
+
2hk
2!
f
xy
(x, y) +
k
2
2!
f
yy
(x, y) +
h
3
3!
f
xxx
(x, y) +
3h
2
k
3!
f
xxy
(x, y)
+
3hk
2
3!
f
xyy
(x, y) +
k
3
3!
f
yyy
(x, y) +
· · ·
(12.1.3)
The remainder can be written in the form
1
n!

h
∂
∂x
+ k
∂
∂y

n
f (x + θh, y + θk),
0
≤ θ ≤ 1.
(12.1.4)
Here we used a subscript to denote partial differentiation. We will be interested in obtaining
approximation about the point (x
i
, y
j
) and we use a subscript to denote the function values
at the point, i.e. f
i j
= f (x
i
, y
j
).
The Taylor series expansion for f
i+1
about the point x
i
is given by
f
i+1
= f
i
+ hf

i
+
h
2
2!
f

i
+
h
3
3!
f

i
+
· · ·
(12.1.5)
The Taylor series expansion for f
i+1 j+1
about the point (x
i
, y
j
) is given by
f
i+1 j+1
= f
ij
+ (h
x
f
x
+ h
y
f
y
)
i j
+ (
h
2
x
2
f
xx
+ h
x
h
y
f
xy
+
h
2
y
2
f
yy
)
i j
+
· · ·
(12.1.6)
Remark: The expansion for f
i+1 j
about (x
i
, y
j
) proceeds as in the case of a function of one
variable.
277

12.2
Finite Differences
An infinite number of difference representations can be found for the partial derivatives of
f (x, y). Let us use the following operators:
forward difference operator
∆
x
f
i j
= f
i+1 j
− f
i j
(12.2.1)
backward difference operator
∇
x
f
i j
= f
i j
− f
i−1 j
(12.2.2)
centered difference
δ
x
f
i j
= f
i+1 j
− f
i−1 j
(12.2.3)
δ
x
f
i j
= f
i+1/2 j
− f
i−1/2 j
(12.2.4)
averaging operator
µ
x
f
i j
= (f
i+1/2 j
+ f
i−1/2 j
)/2
(12.2.5)
Note that
δ
x
= µ
x
δ
x
.
(12.2.6)
In a similar fashion we can define the corresponding operators in y.
In the following table we collected some of the common approximations for the first
derivative.
Finite Difference
Order (see next chapter)
1
h
x
∆
x
f
i j
O(h
x
)
1
h
x
∇
x
f
i j
O(h
x
)
1
2h
x
δ
x
f
i j
O(h
2
x
)
1
2h
x
(
−3f
i j
+ 4f
i+1 j
− f
i+2 j
) =
1
h
x
(∆
x
−
1
2
∆
2
x
)f
i j
O(h
2
x
)
1
2h
x
(3f
i j
− 4f
i−1 j
+ f
i−2 j
) =
1
h
x
(
∇
x
+
1
2
∇
2
x
)f
i j
O(h
2
x
)
1
h
x
(µ
x
δ
x
−
1
3!
µ
x
δ
3
x
)f
i j
O(h
3
x
)
1
2h
x
δ
x
f
i j
1 +
1
6
δ
2
x
O(h
4
x
)
Table 1: Order of approximations to f
x
The compact fourth order three point scheme deserves some explanation. Let f
x
be v,
then the method is to be interpreted as
(1 +
1
6
δ
2
x
)v
i j
=
1
2h
x
δ
x
f
i j
(12.2.7)
or
1
6
(v
i+1 j
+ 4v
i j
+ v
i−1 j
) =
1
2h
x
δ
x
f
i j
.
(12.2.8)
278

This is an implicit formula for the derivative
∂f
∂x
at (x
i
, y
j
). The v
i j
can be computed from
the f
i j
by solving a tridiagonal system of algebraic equations.
The most common second derivative approximations are
f
xx
|
i j
=
1
h
2
x
(f
i j
− 2f
i+1 j
+ f
i+2 j
) + O(h
x
)
(12.2.9)
f
xx
|
i j
=
1
h
2
x
(f
i j
− 2f
i−1 j
+ f
i−2 j
) + O(h
x
)
(12.2.10)
f
xx
|
i j
=
1
h
2
x
δ
2
x
f
i j
+ O(h
2
x
)
(12.2.11)
f
xx
|
i j
=
1
h
2
x
δ
2
x
f
i j
1 +
1
12
δ
2
x
+ O(h
4
x
)
(12.2.12)
Remarks:
1. The order of a scheme is given for a uniform mesh.
2. Tables for difference approximations using more than three points and approximations
of mixed derivatives are given in Anderson, Tannehill and Pletcher (1984 , p.45).
3. We will use the notation
ˆ
δ
2
x
=
δ
2
x
h
2
x
.
(12.2.13)
The centered difference operator can be written as a product of the forward and backward
operator, i.e.
δ
2
x
f
i j
=
∇
x
∆
x
f
i j
.
(12.2.14)
This is true since on the right we have
∇
x
(f
i+1 j
− f
i j
) = f
i+1 j
− f
i j
− (f
i j
− f
i−1 j
)
which agrees with the right hand side of (12.2.14). This idea is important when one wants
to approximate (p(x)y

(x))

at the point x
i
to a second order. In this case one takes the
forward difference inside and the backward difference outside (or vice versa)
∇
x
p
i
y
i+1
− y
i
∆x

(12.2.15)
and after expanding again
p
i
y
i+1
− y
i
∆x
− p
i−1
y
i
− y
i−1
∆x
∆x
(12.2.16)
or
p
i
y
i+1
− (p
i
+ p
i−1
) y
i
+ p
i−1
y
i−1
(∆x)
2
.
(12.2.17)
Note that if p(x)
≡ 1 then we get the well known centered difference.
279

A
+
∆
 x
O
+
α
 
∆
 x
C
+
B+
∆
 y
β
 
∆
 y
D+
y
x
Figure 61: Irregular mesh near curved boundary
12.3
Irregular Mesh
Clearly it is more convenient to use a uniform mesh and it is more accurate in some cases.
However, in many cases this is not possible due to boundaries which do not coincide with the
mesh or due to the need to refine the mesh in part of the domain to maintain the accuracy.
In the latter case one is advised to use a coordinate transformation.
In the former case several possible cures are given in, e.g. Anderson et al (1984). The
most accurate of these is a development of a finite difference approximation which is valid
even when the mesh is nonuniform. It can be shown that
u
xx
O
∼
=
2
(1 + α)h
x
u
c
− u
O
αh
x
−
u
O
− u
A
h
x

(12.3.1)
Similar formula for u
yy
. Note that for α = 1 one obtains the centered difference approx-
imation.
We now develop a three point second order approximation for
∂f
∂x
on a nonuniform mesh.
∂f
∂x
at point O can be written as a linear combination of values of f at A, O, and B,
∂f
∂x
O
= C
1
f (A) + C
2
f (O) + C
3
f (B) .
(12.3.2)
A
+
∆
 x
O
+
α
 
∆
 x
B
+
x
Figure 62: Nonuniform mesh
We use Taylor series to expand f (A) and f (B) about the point O,
f (A) = f (O
− ∆x) = f(O) − ∆xf

(O) +
∆x
2
2
f

(O)
−
∆x
3
6
f

(O)
± · · ·
(12.3.3)
280

f (B) = f (O + α∆x) = f (O) + α∆xf

(O) +
α
2
∆x
2
2
f

(O) +
α
3
∆x
3
6
f

(O) +
· · · (12.3.4)
Thus
∂f
∂x
O
= (C
1
+ C
2
+ C
3
)f (O) + (αC
3
− C
1
)∆x
∂f
∂x
O
+ (C
1
+ α
2
C
3
)
∆x
2
2
∂
2
f
∂x
2
O
+ (α
3
C
3
− C
1
)
∆x
3
6
∂
3
f
∂x
3
O
+
· · ·
(12.3.5)
This yields the following system of equations
C
1
+ C
2
+ C
3
= 0
(12.3.6)
−C
1
+ αC
3
=
1
∆x
(12.3.7)
C
1
+ α
2
C
3
= 0
(12.3.8)
The solution is
C
1
=
−
α
(α + 1)∆x
,
C
2
=
α
− 1
α∆x
,
C
3
=
1
α(α + 1)∆x
(12.3.9)
and thus
∂f
∂x
=
−α
2
f (A) + (α
2
− 1)f(O) + f(B)
α(α + 1)∆x
+
α
6
∆x
2
∂
3
f
∂x
3
O
+
· · ·
(12.3.10)
Note that if the grid is uniform then α = 1 and this becomes the familiar centered difference.
12.4
Thomas Algorithm
This is an algorithm to solve a tridiagonal system of equations





d
1
a
1
b
2
d
2
a
2
b
3
d
3
a
3
. . .





u =





c
1
c
2
c
3
. . .





(12.4.1)
The first step of Thomas algorithm is to bring the tridiagonal M by M matrix to an upper
triangular form
d
i
← d
i
−
b
i
d
i−1
a
i−1
,
i = 2, 3,
· · · , M
(12.4.2)
c
i
← c
i
−
b
i
d
i−1
c
i−1
,
i = 2, 3,
· · · , M.
(12.4.3)
The second step is to backsolve
u
M
=
c
M
d
M
(12.4.4)
u
j
=
c
j
− a
j
u
j+1
d
j
,
j = M
− 1, · · · , 1.
(12.4.5)
The following subroutine solves a tridiagonal system of equations:
281

subroutine tridg(il,iu,rl,d,ru,r)
c
c
solve a tridiagonal system
c
the rhs vector is destroyed and gives the solution
c
the diagonal vector is destroyed
c
integer il,iu
real rl(1),d(1),ru(1),r(1)
C
C
the equations are
C
rl(i)*u(i-1)+d(i)*u(i)+ru(i)*u(i+1)=r(i)
C
il subscript of first equation
C
iu subscript of last equation
C
ilp=il+1
do 1 i=ilp,iu
g=rl(i)/d(i-1)
d(i)=d(i)-g*ru(i-1)
r(i)=r(i)-g*r(i-1)
1 continue
c
c
Back substitution
c
r(iu)=r(iu)/d(iu)
do 2 i=ilp,iu
j=iu-i+il
r(j)=(r(j)-ru(j)*r(j+1))/d(j)
2 continue
return
end
12.5
Methods for Approximating PDEs
In this section we discuss several methods to approximate PDEs. These are certainly not all
the possibilities.
12.5.1
Undetermined coefficients
In this case, we approximate the required partial derivative by a linear combination of
function values. The weights are chosen so that the approximation is of the appropriate
282

order. For example, we can approximate u
xx
at x
i
, y
j
by taking the three neighboring
points,
u
xx
|
i j
= Au
i+1 j
+ Bu
i j
+ Cu
i−1 j
(12.5.1.1)
Now expand each of the terms on the right in Taylor series and compare coefficients (all
terms are evaluated at i j)
u
xx
=
A

u + hu
x
+
h
2
2
u
xx
+
h
3
6
u
xxx
+
h
4
24
u
xxxx
+
· · ·

+ Bu + C

u
− hu
x
+
h
2
2
u
xx
−
h
3
6
u
xxx
+
h
4
24
u
xxxx
± · · ·

(12.5.1.2)
Upon collecting coefficients, we have
A + B + C = 0
(12.5.1.3)
A
− C = 0
(12.5.1.4)
(A + C)
h
2
2
= 1
(12.5.1.5)
This yields
A = C =
1
h
2
(12.5.1.6)
B =
−2
h
2
(12.5.1.7)
The error term, is the next nonzero term, which is
(A + C)
h
4
24
u
xxxx
=
h
2
12
u
xxxx
.
(12.5.1.8)
We call the method second order, because of the h
2
factor in the error term. This is the
centered difference approximation given by (12.2.11).
12.5.2
Integral Method
The strategy here is to develop an algebraic relationship among the values of the unknowns at
neighboring grid points, by integrating the PDE. We demonstrate this on the heat equation
integrated around the point (x
j
, t
n
). The solution at this point can be related to neighboring
values by integration, e.g.


x
j
+∆x/2
x
j
−∆x/2


t
n
+∆t
t
n
u
t
dt

dx = α


t
n
+∆t
t
n


x
j
+∆x/2
x
j
−∆x/2
u
xx
dx

dt.
(12.5.2.1)
Note the order of integration on both sides.


x
j
+∆x/2
x
j
−∆x/2
( u(x, t
n
+ ∆t)
− u(x, t
n
) ) dx = α


t
n
+∆t
t
n
(u
x
(x
j
+ ∆x/2, t)
− u
x
(x
j
− ∆x/2, t)) dt.
(12.5.2.2)
283

Now use the mean value theorem, choosing x
j
as the intermediate point on the left and
t
n
+ ∆t as the intermediate point on the right,
(u(x
j
, t
n
+ ∆t)
− u(x
j
, t
n
) ) ∆x = α (u
x
(x
j
+ ∆x/2, t
n
+ ∆t)
− u
x
(x
j
− ∆x/2, t
n
+ ∆t)) ∆t.
(12.5.2.3)
Now use a centered difference approximation for the u
x
terms and we get the fully implicit
scheme, i.e.
u
n+1
j
− u
n
j
∆t
= α
u
n+1
j+1
− 2u
n+1
j
+ u
n+1
j−1
(∆x)
2
.
(12.5.2.4)