3

Numerical Solutions of PDEs

There’s no sense in being precise when you don’t even know what you’re talking
about.- John von Neumann (1903-1957)

Most of the book has dealt with finding exact solutions to some
generic problems. However, most problems of interest cannot be solved ex-
actly. The heat, wave, and Laplace equations are linear partial differential
equations and can be solved using separation of variables in geometries
in which the Laplacian is separable. However, once we introduce nonlin-
earities, or complicated non-constant coefficients intro the equations, some
of these methods do not work. Even when separation of variables or the
method of eigenfunction expansions gave us exact results, the computation
of the resulting series had to be done on a computer and inevitably one
could only use a finite number of terms of the expansion. So, therefore, it
is sometimes useful to be able to solve differential equations numerically.
In this chapter we will introduce the idea of numerical solutions of partial
differential equations. We will introduce finite difference method and the
idea of stability. Other common approaches may be added later.

3.1 The Finite Difference Method

The heat equation can be solved using separation of variables.

However, many partial differential equations cannot be solved exactly and
one needs to turn to numerical solutions. The heat equation is a simple test
case for using numerical methods. Here we will use the simplest method,
finite differences.
Let us consider the heat equation in one dimension,

ut = ku xx .

Boundary conditions and an initial condition will be applied later. The

starting point is figuring out how to approximate the derivatives in this
equation.
84 partial differential equations

Recall that the partial derivative, ut , is defined by

∂u u( x, t + ∆t) − u( x, t)
= lim .
∂t ∆t→∞ ∆t
Therefore, we can use the approximation

∂u u( x, t + ∆t) − u( x, t)
≈ . (3.1)
∂t ∆t
This is called a forward difference approximation.
In order to find an approximation to the second derivative, u xx , we start
with the forward difference
∂u u( x + ∆x, t) − u( x, t)
≈ .
∂x ∆x
Then,
∂u x u x ( x + ∆x, t) − u x ( x, t)
≈ .
∂x ∆x
We need to approximate the terms in the numerator. It is customary to
use a backward difference approximation. This is given by letting ∆x →
−∆x in the forward difference form,
∂u u( x, t) − u( x − ∆x, t)
≈ . (3.2)
∂x ∆t
Applying this to u x evaluated at x = x and x = x + ∆x, we have

u( x, t) − u( x − ∆x, t)
u x ( x, t) ≈ ,
∆x
and
u( x + ∆x, t) − u( x, t)
u x ( x + ∆x, t) ≈ .
∆x
Inserting these expressions into the approximation for u xx , we have

∂2 u ∂u x
=
∂x2 ∂x
u x ( x + ∆x, t) − u x ( x, t)
≈
∆x
u( x +∆x,t)−u( x,t) u( x,t)−u( x −∆x,t)
∆x ∆x
≈ −
∆x ∆x
u( x + ∆x, t) − 2u( x, t) + u( x − ∆x, t)
= . (3.3)
(∆x )2
This approximation for u xx is called the central difference approximation of
u xx .
Combining Equation (3.1) with (3.3) in the heat equation, we have

u( x, t + ∆t) − u( x, t) u( x + ∆x, t) − 2u( x, t) + u( x − ∆x, t)

≈k .
∆t (∆x )2

Solving for u( x, t + ∆t), we find

u( x, t + ∆t) ≈ u( x, t) + α [u( x + ∆x, t) − 2u( x, t) + u( x − ∆x, t)] , (3.4)

 numerical solutions of pdes 85

∆t
where α = k (∆x )2
.
In this equation we have a way to determine the solution at position x and
time t + ∆t given that we know the solution at three positions, x, x + ∆x,
and x + 2∆x at time t.

u( x, t + ∆t) ≈ u( x, t) + α [u( x + ∆x, t) − 2u( x, t) + u( x − ∆x, t)] . (3.5)

A shorthand notation is usually used to write out finite difference schemes.

The domain of the solution is x ∈ [ a, b] and t ≥ 0. We seek approximate val-
ues of u( x, t) at specific positions and times. We first divide the interval
[ a, b] into N subintervals of width ∆x = (b − a)/N. Then, the endpoints of
the subintervals are given by

xi = a + i∆x, i = 0, 1, . . . , N.

Similarly, we take time steps of ∆t, at times

t j = j∆t, j = 0, 1, 2, . . . .

This gives a grid of points ( xi , t j ) in the domain.

At each grid point in the domain we seek an approximate solution to the
heat equation, ui,j ≈ u( xi , t j ). Equation (3.5) becomes
 
ui,j+1 ≈ ui,j + α ui+1,j − 2ui,j + ui−1,j . (3.6)

t Figure 3.1: This stencil indicates the four

types of terms in the finite difference
scheme in Equation (3.7). The black cir-
cles represent the four terms in the equa-
tion, ui,j ui−1,j ui+1,j and ui,j+1 .

j+1

i−1 i i+1 x

Equation (3.7) is the finite difference scheme for solving the heat equation.
This equation is represented by the stencil shown in Figure 3.1. The black
circles represent the four terms in the equation, ui,j ui−1,j ui+1,j and ui,j+1 .
Let’s assume that the initial condition is given by

u( x, 0) = f ( x ).

Then, we have ui,0 = f ( xi ). Knowing these values, denoted by the open

circles in Figure 3.2, we apply the stencil to generate the solution on the
j = 1 row. This is shown in Figure 3.2.
86 partial differential equations

Figure 3.2: Applying the stencil to the t

row of initial values gives the solution at
the next time step.

Figure 3.3: Continuation of the pro- t

cess provides solutions at the indicated
points.

Further rows are generated by successively applying the stencil on each

row, using the known approximations of ui,j at each level. This gives the
values of the solution at the open circles shown in Figure 3.3. We notice that
the solution can only be obtained at a finite number of points on the grid.
In order to obtain the missing values, we need to impose boundary con-
ditions. For example, if we have Dirichlet conditions at x = a,

u( a, t) = 0,

or u0,j = 0 for j = 0, 1, . . . , then we can fill in some of the missing data

points as seen in Figure 3.4.
The process continues until we again go as far as we can. This is shown
in Figure 3.5.
We can fill in the rest of the grid using a boundary condition at x = b.
For Dirichlet conditions at x = b,

u(b, t) = 0,

or u N,j = 0 for j = 0, 1, . . . , then we can fill in the rest of the missing data
points as seen in Figure 3.6.
 numerical solutions of pdes 87

t Figure 3.4: Knowing the values of the so-

lution at x = a, we can fill in more of the
grid.

t Figure 3.5: Knowing the values of the so-

lution at other times, we continue to fill
the grid as far as the stencil can go.

We could also use Neumann conditions. For example, let

u x ( a, t) = 0.

The approximation to the derivative gives

∂u u( a + ∆x, t) − u( a, t)
≈ = 0.
∂x x=a ∆x
Then,
u( a + ∆x, t) − u( a, t),
or u0,j = u1,j , for j = 0, 1, . . . . Thus, we know the values at the boundary
and can generate the solutions at the grid points as before.
We now have to code this using software. We can use MATLAB to do
this. An example of the code is given below. In this example we specify the
length of the rod, L = 1, and the heat constant, k = 1. The code is run for
t ∈ [0, 0.1].
The grid is created using N = 10 subintervals in space and M = 50 time
steps. This gives dx = ∆x and dt = ∆t. Using these values, we find the
numerical scheme constant α = k∆t/(∆x )2 .
88 partial differential equations

Figure 3.6: Using boundary conditions t

and the initial condition, the grid can be
fill in through any time level.

Nest, we define xi = i ∗ dx, i = 0, 1, . . . , N. However, in MATLAB, we

cannot have an index of 0. We need to start with i = 1. Thus, xi = (i − 1) ∗ dx,
i = 1, 2, . . . , N + 1.
Next, we establish the initial condition. We take a simple condition of

u( x, 0) = sin πx.

We have enough information to begin the numerical scheme as developed

earlier. Namely, we cycle through the time steps using the scheme. There is
one loop for each time step. We will generate the new time step from the
last time step in the form
h i
uinew = uiold + α uiold
+1 − 2u old
i + u old
i −1 . (3.7)

This is done suing u0(i ) = uinew and u1(i ) = uiold .

At the end of each time loop we update the boundary points so that the
grid can be filled in as discussed. When done, we can plot the final solution.
If we want to show solutions at intermediate steps, we can plot the solution
earlier.

% Solution of the Heat Equation Using a Forward Difference Scheme

% Initialize Data
% Length of Rod, Time Interval
% Number of Points in Space, Number of Time Steps

L=1;
T=0.1;
k=1;
N=10;
M=50;
dx=L/N;
dt=T/M;
alpha=k*dt/dx^2;
 numerical solutions of pdes 89

% Position

for i=1:N+1
x(i)=(i-1)*dx;
end

% Initial Condition

for i=1:N+1
u0(i)=sin(pi*x(i));
end

% Partial Difference Equation (Numerical Scheme)

for j=1:M
for i=2:N
u1(i)=u0(i)+alpha*(u0(i+1)-2*u0(i)+u0(i-1));
end
u1(1)=0;
u1(N+1)=0;
u0=u1;
end

% Plot solution
plot(x, u1);

3.2 Truncation Error

In the previous section we found a finite difference scheme for

numerically solving the one dimensional heat equation. We have from Equa-
tions (3.5) and (3.7),

u( x, t + ∆t) ≈ u( x, t) + α [u( x + ∆x, t) − 2u( x, t) + u( x − ∆x, t)] .(3.8)

 
ui,j+1 ≈ ui,j + α ui+1,j − 2ui,j + ui−1,j , (3.9)

where α = k∆t/(∆x )2 . For points x ∈ [ a, b] and t ≥ 0, we use the scheme

to find approximate values of u( xi , ti ) = ui,j at positions xi = a + i∆x, i =
0, 1, . . . , N, and times t j = j∆t, j = 0, 1, 2, . . . .
In implementing the scheme we have found that there are errors intro-
duced just like when using Euler’s Method for ordinary differential equa-
tions. These truncations errors can be found by applying Taylor approxi-
mations just like we had for ordinary differential equations. In the schemes
(3.8) and (3.9), we have not use equality. In order to replace the approxi-
mation by an equality, we need t estimate the order of the terms neglected
90 partial differential equations

in a Taylor series approximation of the time and space derivatives we have

approximated.
We begin with the time derivative approximation. We used the forward
difference approximation (3.1),
∂u u( x, t + ∆t) − u( x, t)
≈ . (3.10)
∂t ∆t
This can be derived from the Tayloer series expansion of u( x, t + ∆t) about
∆t = 0,
∂u 1 ∂2 u
u( x, t + ∆t) = u( x, t) + ( x, t)∆t + ( x, t)(∆t)2 + O((∆t)3 ).
∂t 2! ∂t2
∂u
Solving for ∂t ( x, t ), we obtain

∂u u( x, t + ∆t) − u( x, t) 1 ∂2 u
( x, t) = − ( x, t)∆t + O((∆t)2 ).
∂t ∆t 2! ∂t2
We see that we have obtained the forward difference approximation (3.1)
with the added benefit of knowing something about the error terms intro-
duced in the approximation. Namely, when we approximate ut with the
forward difference approximation (3.1), we are making an error of
1 ∂2 u
E( x, t, ∆t) = − ( x, t)∆t + O((∆t)2 ).
2! ∂t2
We have truncated the Taylor series to obtain this approximation and we say
that
∂u u( x, t + ∆t) − u( x, t)
= + O(∆t) (3.11)
∂t ∆t
is a first order approximation in ∆t.
In a similar manor, we can obtain the truncation error for the u x x- term.
However, instead of starting with the approximation we used in Equation
??uxx), we will derive a term using the Taylor series expansion of u( x +
∆x, t) about ∆x = 0. Namely, we begin with the expansion
1 1
u( x + ∆x, t) = u( x, t) + u x ( x, t)∆x + u xx ( x, t)(∆x )2 + u xxx ( x, t)(∆x )3
2! 3!
1
+ u xxxx ( x, t)(∆x )4 + . . . . (3.12)
4!
We want to solve this equation for u xx . However, there are some obstruc-
tions, like needing to know the u x term. So, we seek a way to eliminate
lower order terms. On way is to note that replacing ∆x by −∆x gives
1 1
u( x − ∆x, t) = u( x, t) − u x ( x, t)∆x + u xx ( x, t)(∆x )2 − u xxx ( x, t)(∆x )3
2! 3!
1
+ u xxxx ( x, t)(∆x )4 + . . . . (3.13)
4!
Adding these Taylor series, we have

u( x + ∆x, t) + u( x + ∆x, t) = 2u( x, t) + u xx ( x, t)(∆x )2

2
+ u xxxx ( x, t)(∆x )4 + O((∆x )6 ).
4!
(3.14)
 numerical solutions of pdes 91

We can now solve for u xx to find

u( x + ∆x, t) − 2u( x, t) + u( x + ∆x, t)
u xx ( x, t) =
(∆x )2
2
+ u xxxx ( x, t)(∆x )2 + O((∆x )4 ). (3.15)
4!
Thus, we have that
u( x + ∆x, t) − 2u( x, t) + u( x + ∆x, t)
u xx ( x, t) = + O((∆x )2 )
(∆x )2
is the second order in ∆x approximation of u xx .
Combining these results, we find that the heat equation is approximated
by
u( x, t + ∆t) − u( x, t) u( x + ∆x, t) − 2u( x, t) + u( x + ∆x, t) 
2

= + O ( ∆x ) , ∆t .
∆t (∆x )2
This has local truncation error that is first order in time and second order in
space.

3.3 Stability

Another consideration for numerical schemes for the heat equa-

tion is the stability of the scheme. In implementing the finite difference
scheme,
 
um,j+1 = um,j + α um+1,j − 2um,j + um−1,j , (3.16)
α = k∆t/(∆x )2 , one finds that the solution goes crazy when α is too big.
In other words, if you try to push the individual time steps too far into
the future, then something goes haywire. We can determine the onset of
instability by looking at the solution of this equation for um,j . [Note: We
changed index i to m to avoid confusion later in this section.]
The scheme is actually what is called a partial difference equation for
um,j . We could write it in terms of difference, such as um+1,j − um,j and
um,j+1 − um,j . The furthest apart the time steps are are one unit and the
spatial points are two units apart. We can see this in the stencils in Figure
3.1. So, this is a second order partial difference equation similar to the idea
that the heat equation is a second order partial differential equation. The
heat equation can be solved using the method of separation of variables.
The difference scheme can also be solved in a similar fashion. We will show
how this can lead to product solutions.
We begin by assuming that umj = Xm Tj , a product of functions of the
indices m and j. Inserting this guess into the finite difference equation, we
have
 
um,j+1 = um,j + α um+1,j − 2um,j + um−1,j ,
Xm Tj+1 = Xm Tj + α [ Xm+1 − 2Xm + Xm−1 ] Tj ,
Tj+1 αXm+1 + (1 − 2α) Xm + αXm−1
= . (3.17)
Tj Xm
92 partial differential equations

Noting that we have a function of j equal to a function of m, then we can

set each of these to a constant, λ. Then, we obtain two ordinary difference
equations:

Tj+1 = λTj , (3.18)

αXm+1 + (1 − 2α) Xm + αXm−1 = λXm . (3.19)

The first equation is a simple first order difference equation and can be
solved by iteration:

Tj+1 = λTj ,
= λ(λTj−1 ) = λ2 Tj−1 ,
= λ3 Tj−2 ,
= λ j+1 T0 , (3.20)

The second difference equation can be solved by making a guess in the

same spirit as solving a second order constant coefficient differential equa-
tion.Namely, let Xm = ξ m for some number ξ. This gives

αXm+1 + (1 − 2α) Xm + αXm−1 = λXm ,

h i
ξ m−1 αξ 2 + (1 − 2α)ξ + α = λξ m
αξ 2 + (1 − 2α − λ)ξ + α = 0. (3.21)

This is an equation foe ξ in terms of α and λ. Due to the boundary

conditions, we expect to have oscillatory solutions. So, we can guess that
ξ = |ξ |eiθ , where i here is the imaginary unit. We assume that |ξ | = 1, and
thus ξ = eiθ and Xm = ξ m = eimθ . Since xm = m∆x, we have Xm = eixm θ/∆x .
We define β = theta/∆, to give Xm = eiβxm and ξ = eiβ∆x .
Inserting this value for ξ into the quadratic equation for ξ, we have

0 = αξ 2 + (1 − 2α − λ)ξ + α
= αe2iβ∆x + (1 − 2α − λ)eiβ∆x + α
h i
= eiβ∆x α(eiβ∆x + e−iβ∆x ) + (1 − 2α − λ)
= eiβ∆x [2α cos( β∆x ) + (1 − 2α − λ)]
λ = 2α cos( β∆x ) + 1 − 2α. (3.22)

So, we have found that

umj = Xm Tj = λm ( a cos αxm + b sin αxm ), a2 + b2 = h20 ,

and
λ = 2α cos( β∆x ) + 1 − 2α.

For the solution to remain bounded, or stable, we need |λ| ≤ 1.

Therefore, we have the inequality

−1 ≤ 2α cos( β∆x ) + 1 − 2α ≤ 1.
 numerical solutions of pdes 93

Since cos( β∆x ) ≤ 1, the upper bound is obviously satisfied. Since −1 ≤

cos( β∆x ), the lower bound is satisfied for −1 ≤ −2s + 1 − 2s, or s ≤ 12 .
Therefore, the stability criterion is satisfied when

∆t 1
α=k ≤ .
∆x 2 2