10

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. However, we will first begin with a discussion of the
solution of ordinary differential equations in order to get a feel for some
common problems in the solution of differential equations and the notion
of convergence rates of numerical schemes. Then, we turn to the finite
difference method and the ideas of stability. Other common approaches
may be added later.

10.1 Ordinary Differential Equations

10.1.1 Euler’s Method

In this section we will look at the simplest method for solving

first order equations, Euler’s Method. While it is not the most efficient
method, it does provide us with a picture of how one proceeds and can be
improved by introducing better techniques, which are typically covered in
a numerical analysis text.
Let’s consider the class of first order initial value problems of the form
dy
= f ( x, y), y ( x0 ) = y0 . (10.1)
dx
432 partial differential equations

We are interested in finding the solution y( x ) of this equation which passes

through the initial point ( x0 , y0 ) in the xy-plane for values of x in the interval
[ a, b], where a = x0 . We will seek approximations of the solution at N
points, labeled xn for n = 1, . . . , N. For equally spaced points we have
∆x = x1 − x0 = x2 − x1 , etc. We can write these as

xn = x0 + n∆x.

In Figure 10.1 we show three such points on the x-axis.

Figure 10.1: The basics of Euler’s y

Method are shown. An interval of the
x axis is broken into N subintervals.
The approximations to the solutions are
( x2 , y( x2 ))
found using the slope of the tangent to
the solution, given by f ( x, y). Knowing
previous approximations at ( xn−1 , yn−1 ), y2
one can determine the next approxima-
tion, yn . ( x1 , y( x1 ))
y1

y0
( x0 , y0 )

x
x0 x1 x2

The first step of Euler’s Method is to use the initial condition. We repre-
sent this as a point on the solution curve, ( x0 , y( x0 )) = ( x0 , y0 ), as shown in
Figure 10.1. The next step is to develop a method for obtaining approxima-
tions to the solution for the other xn ’s.
We first note that the differential equation gives the slope of the tangent
line at ( x, y( x )) of the solution curve since the slope is the derivative, y0 ( x )0
From the differential equation the slope is f ( x, y( x )). Referring to Figure
10.1, we see the tangent line drawn at ( x0 , y0 ). We look now at x = x1 . The
vertical line x = x1 intersects both the solution curve and the tangent line
passing through ( x0 , y0 ). This is shown by a heavy dashed line.
While we do not know the solution at x = x1 , we can determine the
tangent line and find the intersection point that it makes with the vertical.
As seen in the figure, this intersection point is in theory close to the point
on the solution curve. So, we will designate y1 as the approximation of the
solution y( x1 ). We just need to determine y1 .
The idea is simple. We approximate the derivative in the differential
equation by its difference quotient:
dy y − y0 y − y0
≈ 1 = 1 . (10.2)
dx x1 − x0 ∆x
Since the slope of the tangent to the curve at ( x0 , y0 ) is y0 ( x0 ) = f ( x0 , y0 ),
 numerical solutions of pdes 433

we can write
y1 − y0
≈ f ( x0 , y0 ). (10.3)
∆x
Solving this equation for y1 , we obtain

y1 = y0 + ∆x f ( x0 , y0 ). (10.4)

This gives y1 in terms of quantities that we know.

We now proceed to approximate y( x2 ). Referring to Figure 10.1, we see
that this can be done by using the slope of the solution curve at ( x1 , y1 ).
The corresponding tangent line is shown passing though ( x1 , y1 ) and we
can then get the value of y2 from the intersection of the tangent line with a
vertical line, x = x2 . Following the previous arguments, we find that

y2 = y1 + ∆x f ( x1 , y1 ). (10.5)

Continuing this procedure for all xn , n = 1, . . . N, we arrive at the fol-

lowing scheme for determining a numerical solution to the initial value
problem:

y0 = y ( x0 ),
yn = yn−1 + ∆x f ( xn−1 , yn−1 ), n = 1, . . . , N. (10.6)

This is referred to as Euler’s Method.

dy
Example 10.1. Use Euler’s Method to solve the initial value problem dx = x +
y, y(0) = 1 and obtain an approximation for y(1).
First, we will do this by hand. We break up the interval [0, 1], since we want the
solution at x = 1 and the initial value is at x = 0. Let ∆x = 0.50. Then, x0 = 0,
x1 = 0.5 and x2 = 1.0. Note that there are N = b∆x − a = 2 subintervals and thus

three points.
We next carry out Euler’s Method systematically by setting up a table for the
needed values. Such a table is shown in Table 10.1. Note how the table is set up.
There is a column for each xn and yn . The first row is the initial condition. We also
made use of the function f ( x, y) in computing the yn ’s from (10.6). This sometimes
makes the computation easier. As a result, we find that the desired approximation
is given as y2 = 2.5.

n xn yn = yn−1 + ∆x f ( xn−1 , yn−1 = 0.5xn−1 + 1.5yn−1 Table 10.1: Application of Euler’s

Method for y0 = x + y, y(0) = 1 and
0 0 1 ∆x = 0.5.
1 0.5 0.5(0) + 1.5(1.0) = 1.5
2 1.0 0.5(0.5) + 1.5(1.5) = 2.5

Is this a good result? Well, we could make the spatial increments smaller. Let’s
repeat the procedure for ∆x = 0.2, or N = 5. The results are in Table 10.2.
Now we see that the approximation is y1 = 2.97664. So, it looks like the value
is near 3, but we cannot say much more. Decreasing ∆x more shows that we are
beginning to converge to a solution. We see this in Table 10.3.
434 partial differential equations

Table 10.2: Application of Euler’s n xn yn = 0.2xn−1 + 1.2yn−1

Method for y0 = x + y, y(0) = 1 and
∆x = 0.2.
0 0 1
1 0.2 0.2(0) + 1.2(1.0) = 1.2
2 0.4 0.2(0.2) + 1.2(1.2) = 1.48
3 0.6 0.2(0.4) + 1.2(1.48) = 1.856
4 0.8 0.2(0.6) + 1.2(1.856) = 2.3472
5 1.0 0.2(0.8) + 1.2(2.3472) = 2.97664

Table 10.3: Results of Euler’s Method for ∆x y N ≈ y (1)

y0 = x + y, y(0) = 1 and varying ∆x
0.5 2.5
0.2 2.97664
0.1 3.187484920
0.01 3.409627659
0.001 3.433847864
0.0001 3.436291854

Of course, these values were not done by hand. The last computation
would have taken 1000 lines in the table, or at least 40 pages! One could
use a computer to do this. A simple code in Maple would look like the
following:

> restart:
> f:=(x,y)->y+x;
> a:=0: b:=1: N:=100: h:=(b-a)/N;
> x[0]:=0: y[0]:=1:
for i from 1 to N do
y[i]:=y[i-1]+h*f(x[i-1],y[i-1]):
x[i]:=x[0]+h*(i):
od:
evalf(y[N]);

In this case we could simply use the exact solution. The exact solution is
easily found as
y( x ) = 2e x − x − 1.
(The reader can verify this.) So, the value we are seeking is

y(1) = 2e − 2 = 3.4365636 . . . .

Thus, even the last numerical solution was off by about 0.00027.
Adding a few extra lines for plotting, we can visually see how well the
approximations compare to the exact solution. The Maple code for doing
such a plot is given below.

> with(plots):
Figure 10.2: A comparison of the results
> Data:=[seq([x[i],y[i]],i=0..N)]:
Euler’s Method to the exact solution for
y0 = x + y, y(0) = 1 and N = 10. > P1:=pointplot(Data,symbol=DIAMOND):
> Sol:=t->-t-1+2*exp(t);
 numerical solutions of pdes 435

> P2:=plot(Sol(t),t=a..b,Sol=0..Sol(b)):
> display({P1,P2});

We show in Figures 10.2-10.3 the results for N = 10 and N = 100. In

Figure 10.2 we can see how quickly the numerical solution diverges from
the exact solution. In Figure 10.3 we can see that visually the solutions
agree, but we note that from Table 10.3 that for ∆x = 0.01, the solution is
still off in the second decimal place with a relative error of about 0.8%.

Why would we use a numerical method when we have the exact solu-
tion? Exact solutions can serve as test cases for our methods. We can make
sure our code works before applying them to problems whose solution is
not known.
Figure 10.3: A comparison of the results
There are many other methods for solving first order equations. One Euler’s Method to the exact solution for
commonly used method is the fourth order Runge-Kutta method. This y0 = x + y, y(0) = 1 and N = 100.

method has smaller errors at each step as compared to Euler’s Method.

It is well suited for programming and comes built-in in many packages like
Maple and MATLAB. Typically, it is set up to handle systems of first order
equations.
In fact, it is well known that nth order equations can be written as a sys-
tem of n first order equations. Consider the simple second order equation

y00 = f ( x, y).

This is a larger class of equations than the second order constant coefficient
equation. We can turn this into a system of two first order differential equa-
tions by letting u = y and v = y0 = u0 . Then, v0 = y00 = f ( x, u). So, we have
the first order system

u0 = v,
0
v = f ( x, u). (10.7)

We will not go further into the Runge-Kutta Method here. You can find
more about it in a numerical analysis text. However, we will see that systems
of differential equations do arise naturally in physics. Such systems are
often coupled equations and lead to interesting behaviors.

10.1.2 Higher Order Taylor Methods

Euler’s Method for solving differential equations is easy to un-

derstand but is not efficient in the sense that it is what is called a first order
method. The error at each step, the local truncation error, is of order ∆x,
for x the independent variable. The accumulation of the local truncation er-
rors results in what is called the global error. In order to generalize Euler’s
Method, we need to rederive it. Also, since these methods are typically used
for initial value problems, we will cast the problem to be solved as
dy
= f (t, y), y ( a ) = y0 , t ∈ [ a, b]. (10.8)
dt
436 partial differential equations

The first step towards obtaining a numerical approximation to the solu-

tion of this problem is to divide the t-interval, [ a, b], into N subintervals,

ti = a + ih, i = 0, 1, . . . , N, t0 = a, t N = b,
y
where
b−a
. h=
N
We then seek the numerical solutions
(ti , y(ti ))
ỹi ≈ y(ti ), i = 1, 2, . . . , N,

with ỹ0 = y(t0 ) = y0 . Figure 10.4 graphically shows how these quantities
are related.
( a, y0 ) Euler’s Method can be derived using the Taylor series expansion of of the
(ti , ỹi )
solution y(ti + h) about t = ti for i = 1, 2, . . . , N. This is given by
t
t0 ti tN y ( t i +1 ) = y ( ti + h )
Figure 10.4: The interval [ a, b] is divided h2 00
into N equally spaced subintervals. The = y ( ti ) + y 0 ( ti ) h + y ( ξ i ), ξ i ∈ ( t i , t i +1 ). (10.9)
2
exact solution y(ti ) is shown with the
numerical solution, ỹi with ti = a + ih, 2
Here the term h2 y00 (ξ i ) captures all of the higher order terms and represents
i = 0, 1, . . . , N.
the error made using a linear approximation to y(ti + h).
Dropping the remainder term, noting that y0 (t) = f (t, y), and defining
the resulting numerical approximations by ỹi ≈ y(ti ), we have

ỹi+1 = ỹi + h f (ti , ỹi ), i = 0, 1, . . . , N − 1,

ỹ0 = y ( a ) = y0 . (10.10)

This is Euler’s Method.

Euler’s Method is not used in practice since the error is of order h. How-
ever, it is simple enough for understanding the idea of solving differential
equations numerically. Also, it is easy to study the numerical error, which
we will show next.
The error that results for a single step of the method is called the local
truncation error, which is defined by

y(ti+1 ) − ỹi
τi+1 (h) = − f ( t i , y i ).
h
A simple computation gives

h 00
τi+1 (h) = y ( ξ i ), ξ i ∈ ( t i , t i +1 ).
2
Since the local truncation error is of order h, this scheme is said to be of
order one. More generally, for a numerical scheme of the form

ỹi+1 = ỹi + hF (ti , ỹi ), i = 0, 1, . . . , N − 1,

ỹ0 = y ( a ) = y0 , (10.11)

The local truncation error. the local truncation error is defined by

 numerical solutions of pdes 437

y(ti+1 ) − ỹi
τi+1 (h) = − F ( t i , y i ).
h
The accumulation of these errors leads to the global error. In fact, one
can show that if f is continuous, satisfies the Lipschitz condition,

| f (t, y2 ) − f (t, y1 )| ≤ L|y2 − y1 |

for a particular domain D ⊂ R2 , and

|y00 (t)| ≤ M, t ∈ [ a, b],

then
hM  L(ti − a) 
|y(ti ) − ỹ| ≤ e − 1 , i = 0, 1, . . . , N.
2L
Furthermore, if one introduces round-off errors, bounded by δ, in both the
initial condition and at each step, the global error is modified as
 
1 hM δ  L ( ti − a ) 
|y(ti ) − ỹ| ≤ + e − 1 + |δ0 |e L(ti −a) , i = 0, 1, . . . , N.
L 2 h

Then for small enough steps h, there is a point when the round-off error
will dominate the error. [See Burden and Faires, Numerical Analysis for the
details.]
Can we improve upon Euler’s Method? The natural next step towards
finding a better scheme would be to keep more terms in the Taylor series
expansion. This leads to Taylor series methods of order n.
Taylor series methods of order n take the form

ỹi+1 = ỹi + hT (n) (ti , ỹi ), i = 0, 1, . . . , N − 1,

ỹ0 = y0 , (10.12)

where we have defined

h h ( n −1) ( n )
T (n) (t, y) = y0 (t) + y00 (t) + · · · + y ( t ).
2 n!
However, since y0 (t) = f (t, y), we can write

h 0 h ( n −1) ( n −1)
T (n) (t, y) = f (t, y) + f (t, y) + · · · + f (t, y).
2 n!
We note that for n = 1, we retrieve Euler’s Method as a special case. We
demonstrate a third order Taylor’s Method in the next example.

Example 10.2. Apply the third order Taylor’s Method to

dy
= t + y, y (0) = 1
dt
and obtain an approximation for y(1) for h = 0.1.
The third order Taylor’s Method takes the form

ỹi+1 = ỹi + hT (3) (ti , ỹi ), i = 0, 1, . . . , N − 1,

ỹ0 = y0 , (10.13)
438 partial differential equations

where
h 0 h2 00
T (3) (t, y) = f (t, y) + f (t, y) + f (t, y)
2 3!
and f (t, y) = t + y(t).
In order to set up the scheme, we need the first and second derivative of f (t, y) :

d
f 0 (t, y) = (t + y)
dt
= 1 + y0
= 1+t+y (10.14)
d
f 00 (t, y) = (1 + t + y )
dt
= 1 + y0
= 1+t+y (10.15)

Inserting these expressions into the scheme, we have

h2
 
h
ỹi+1 = ỹi + h (ti + yi ) + (1 + ti + yi ) + (1 + ti + yi ) ,
2 3!
1 h
= ỹi + h(ti + yi ) + h2 ( + )(1 + ti + yi ),
2 6
ỹ0 = y0 , (10.16)

for i = 0, 1, . . . , N − 1.
In Figure 10.2 we show the results comparing Euler’s Method, the 3rd Order
Taylor’s Method, and the exact solution for N = 10. In Table 10.4 we provide are
the numerical values. The relative error in Euler’s method is about 7% and that
of the 3rd Order Taylor’s Method is about 0.006%. Thus, the 3rd Order Taylor’s
Method is significantly better than Euler’s Method.

Table 10.4: Numerical values for Euler’s Euler Taylor Exact

Method, 3rd Order Taylor’s Method, and
exact solution for solving Example 10.2
1.0000 1.0000 1.0000
with N = 10.. 1.1000 1.1103 1.1103
1.2200 1.2428 1.2428
1.3620 1.3997 1.3997
1.5282 1.5836 1.5836
1.7210 1.7974 1.7974
1.9431 2.0442 2.0442
2.1974 2.3274 2.3275
2.4872 2.6509 2.6511
2.8159 3.0190 3.0192
3.1875 3.4364 3.4366

In the last section we provided some Maple code for performing Euler’s
method. A similar code in MATLAB looks like the following:

a=0;
b=1;
 numerical solutions of pdes 439

N=10;
h=(b-a)/N;

% Slope function
f = inline(’t+y’,’t’,’y’);
sol = inline(’2*exp(t)-t-1’,’t’);

% Initial Condition
t(1)=0;
y(1)=1;

y
% Euler’s Method
for i=2:N+1 4

y(i)=y(i-1)+h*f(t(i-1),y(i-1));
t(i)=t(i-1)+h; 3

end
2

1
A simple modification can be made for the 3rd Order Taylor’s Method by
replacing the Euler’s method part of the preceding code by
t
0 .2 .4 .5 .8 1
Figure 10.5: Numerical results for Eu-
ler’s Method (filled circle) and 3rd Order
% Taylor’s Method, Order 3 Taylor’s Method (open circle) for solving
Example 10.2 as compared to exact solu-
y(1)=1;
tion (solid line).
h3 = h^2*(1/2+h/6);
for i=2:N+1
y(i)=y(i-1)+h*f(t(i-1),y(i-1))+h3*(1+t(i-1)+y(i-1));
t(i)=t(i-1)+h;
end

While the accuracy in the last example seemed sufficient, we have to re-
member that we only stopped at one unit of time. How can we be confident
that the scheme would work as well if we carried out the computation for
much longer times. For example, if the time unit were only a second, then
one would need 86,400 times longer to predict a day forward. Of course, the
scale matters. But, often we need to carry out numerical schemes for long
times and we hope that the scheme not only converges to a solution, but
that it converges to the solution to the given problem. Also, the previous
example was relatively easy to program because we could provide a rela-
tively simple form for T (3) (t, y) with a quick computation of the derivatives
of f (t, y). This is not always the case and higher order Taylor methods in
this form are not typically used. Instead, one can approximate T (n) (t, y) by
evaluating the known function f (t, y) at selected values of t and y, leading
to Runge-Kutta methods.
440 partial differential equations

10.1.3 Runge-Kutta Methods

As we had seen in the last section, we can use higher order Taylor
methods to derive numerical schemes for solving

dy
= f (t, y), y ( a ) = y0 , t ∈ [ a, b], (10.17)
dt
using a scheme of the form

ỹi+1 = ỹi + hT (n) (ti , ỹi ), i = 0, 1, . . . , N − 1,

ỹ0 = y0 , (10.18)

where we have defined

h h ( n −1) ( n )
T (n) (t, y) = y0 (t) + y00 (t) + · · · + y ( t ).
2 n!

In this section we will find approximations of T (n) (t, y) which avoid the
need for computing the derivatives.
For example, we could approximate

h
T (2) (t, y) = f (t, y) + f racd f dt(t, y)
2
by
T (2) (t, y) ≈ a f (t + α, y + β)
for selected values of a, α, and β. This requires use of a generalization of
Taylor’s series to functions of two variables. In particular, for small α and β
we have

∂f ∂f
a f (t + α, y + β) = a f (t, y) + (t, y)α + (t, y) β
∂t ∂y
 2 2 ∂2 f

1 ∂ f 2 ∂ f 2
+ (t, y)α + 2 (t, y)αβ + 2 (t, y) β
2 ∂t2 ∂t∂y ∂y
+ higher order terms. (10.19)
df
Furthermore, we need dt (t, y). Since y = y(t), this can be found using a
generalization of the Chain Rule from Calculus III:

df ∂f ∂ f dy
(t, y) = + .
dt ∂t ∂y dt

Thus,  
h ∂f ∂ f dy
T (2) (t, y) = f (t, y) + + .
2 ∂t ∂y dt
Comparing this expression to the linear (Taylor series) approximation of
a f (t + α, y + β), we have

T (2) ≈ a f (t + α, y + β)
h ∂f h ∂f ∂f ∂f
f+ + f ≈ a f + aα + β . (10.20)
2 ∂t 2 ∂y ∂t ∂y
 numerical solutions of pdes 441

We see that we can choose

h h
a = 1, α= , β= f.
2 2
This leads to the numerical scheme
 
h h
ỹi+1 = ỹi + h f ti + , ỹi + f (ti , ỹi ) , i = 0, 1, . . . , N − 1,
2 2
ỹ0 = y0 , (10.21)

This Runge-Kutta scheme is called the Midpoint Method, or Second Order

Runge-Kutta Method, and it has order 2 if all second order derivatives of
f (t, y) are bounded. The Midpoint or Second Order Runge-
Often, in implementing Runge-Kutta schemes, one computes the argu- Kutta Method.

ments separately as shown in the following MATLAB code snippet. (This

code snippet could replace the Euler’s Method section in the code in the last
section.)

% Midpoint Method
y(1)=1;
for i=2:N+1
k1=h/2*f(t(i-1),y(i-1));
k2=h*f(t(i-1)+h/2,y(i-1)+k1);
y(i)=y(i-1)+k2;
t(i)=t(i-1)+h;
end

Example 10.3. Compare the Midpoint Method with the 2nd Order Taylor’s Method
for the problem
y0 = t2 + y, y(0) = 1, t ∈ [0, 1]. (10.22)
The solution to this problem is y(t) = 3et − 2 − 2t − t2 . In order to implement
the 2nd Order Taylor’s Method, we need
h 0
T (2) = f (t, y) + f (t, y)
2
h
= t2 + y + (2t + t2 + y). (10.23)
2
The results of the implementation are shown in Table 10.3.
There are other way to approximate higher order Taylor polynomials. For
example, we can approximate T (3) (t, y) using four parameters by

T (3) (t, y) ≈ a f (t, y) + b f (t + α, y + β f (t, y).

Expanding this approximation and using

h df h2 d f
T (3) (t, y) ≈ f (t, y) + (t, y) + (t, y),
2 dt 6 dt
we find that we cannot get rid of O(h2 ) terms. Thus, the best we can do is
derive second order schemes. In fact, following a procedure similar to the
derivation of the Midpoint Method, we find that
h
a + b = 1, , αb = , β = α.
2
442 partial differential equations

Table 10.5: Numerical values for 2nd Or- Exact Taylor Error Midpoint Error
der Taylor’s Method, Midpoint Method,
exact solution, and errors for solving Ex-
1.0000 1.0000 0.0000 1.0000 0.0000
ample 10.3 with N = 10.. 1.1055 1.1050 0.0005 1.1053 0.0003
1.2242 1.2231 0.0011 1.2236 0.0006
1.3596 1.3577 0.0019 1.3585 0.0010
1.5155 1.5127 0.0028 1.5139 0.0016
1.6962 1.6923 0.0038 1.6939 0.0023
1.9064 1.9013 0.0051 1.9032 0.0031
2.1513 2.1447 0.0065 2.1471 0.0041
2.4366 2.4284 0.0083 2.4313 0.0053
2.7688 2.7585 0.0103 2.7620 0.0068
3.1548 3.1422 0.0126 3.1463 0.0085

There are three equations and four unknowns. Therefore there are many
second order methods. Two classic methods are given by the modified Euler
method (a = b = 12 , α = β = h) and Huen’s method (a = 14 , b = 34 ,
The Fourth Order Runge-Kutta. α = β = 32 h).
The Fourth Order Runge-Kutta Method, which is most often used, is
given by the scheme

ỹ0 = y0 ,
k1 = h f (ti , ỹi ),
h 1
k2 = h f (ti + , ỹi + k1 ),
2 2
h 1
k3 = h f (ti + , ỹi + k2 ),
2 2
k4 = h f (ti + h, ỹi + k3 ),
1
ỹi+1 = ỹi + (k1 + 2k2 + 2k3 + k4 ), i = 0, 1, . . . , N − 1. (10.24)
6
Again, we can test this on Example 10.3 with N = 10. The MATLAB
implementation is given by

% Runge-Kutta 4th Order to solve dy/dt = f(t,y), y(a)=y0, on [a,b]

clear

a=0;
b=1;
N=10;
h=(b-a)/N;

% Slope function
f = inline(’t^2+y’,’t’,’y’);
sol = inline(’-2-2*t-t^2+3*exp(t)’,’t’);

% Initial Condition
t(1)=0;
 numerical solutions of pdes 443

y(1)=1;

% RK4 Method
y1(1)=1;
for i=2:N+1
k1=h*f(t(i-1),y1(i-1));
k2=h*f(t(i-1)+h/2,y1(i-1)+k1/2);
k3=h*f(t(i-1)+h/2,y1(i-1)+k2/2);
k4=h*f(t(i-1)+h,y1(i-1)+k3);
y1(i)=y1(i-1)+(k1+2*k2+2*k3+k4)/6;
t(i)=t(i-1)+h;
end
MATLAB has built-in ODE solvers, as do
MATLAB has built-in ODE solvers, such as ode45 for a fourth order other software packages, like Maple and
Mathematica. You should also note that
Runge-Kutta method. Its implementation is given by there are currently open source pack-
ages, such as Python based NumPy and
[t,y]=ode45(f,[0 1],1); Matplotlib, or Octave, of which some
packages are contained within the Sage
Project.
In this case f is given by an inline function like in the above RK4 code.
The time interval is entered as [0, 1] and the 1 is the initial condition, y(0) =
1.
However, ode45 is not a straight forward RK4 implementation. It is a
hybrid method in which a combination of 4th and 5th order methods are
combined allowing for adaptive methods to handled subintervals of the in-
tegration region which need more care. In this case, it implements a fourth
order Runge-Kutta-Fehlberg method. Running this code for the above ex-
ample actually results in values for N = 41 and not N = 10. If we wanted
to have the routine output numerical solutions at specific times, then one
could use the following form

tspan=0:h:1;
[t,y]=ode45(f,tspan,1);

In Table 10.6 we show the solutions which results for Example 10.3 com-
paring the RK4 snippet above with ode45. As you can see RK4 is much
better than the previous implementation of the second order RK (Midpoint)
Method. However, the MATLAB routine is two orders of magnitude better
that RK4.
There are many ODE solvers in MATLAB. These are typically useful if
RK4 is having difficulty solving particular problems. For the most part, one
is fine using RK4, especially as a starting point. For example, there is ode23,
which is similar to ode45 but combining a second and third order scheme.
Applying the results to Example 10.3 we obtain the results in Table 10.6. We
compare these to the second order Runge-Kutta method. The code snippets
are shown below.

% Second Order RK Method

y1(1)=1;
444 partial differential equations

Table 10.6: Numerical values for Fourth Exact Taylor Error Midpoint Error
Order Runge-Kutta Method, rk45, exact
solution, and errors for solving Example
1.0000 1.0000 0.0000 1.0000 0.0000
10.3 with N = 10. 1.1055 1.1055 4.5894e-08 1.1055 -2.5083e-10
1.2242 1.2242 1.2335e-07 1.2242 -6.0935e-10
1.3596 1.3596 2.3850e-07 1.3596 -1.0954e-09
1.5155 1.5155 3.9843e-07 1.5155 -1.7319e-09
1.6962 1.6962 6.1126e-07 1.6962 -2.5451e-09
1.9064 1.9064 8.8636e-07 1.9064 -3.5651e-09
2.1513 2.1513 1.2345e-06 2.1513 -4.8265e-09
2.4366 2.4366 1.6679e-06 2.4366 -6.3686e-09
2.7688 2.7688 2.2008e-06 2.7688 -8.2366e-09
3.1548 3.1548 2.8492e-06 3.1548 -1.0482e-08

for i=2:N+1
k1=h*f(t(i-1),y1(i-1));
k2=h*f(t(i-1)+h/2,y1(i-1)+k1/2);
y1(i)=y1(i-1)+k2;
t(i)=t(i-1)+h;
end

tspan=0:h:1;
[t,y]=ode23(f,tspan,1);

Table 10.7: Numerical values for Second Exact Taylor Error Midpoint Error
Order Runge-Kutta Method, rk23, exact
solution, and errors for solving Example
1.0000 1.0000 0.0000 1.0000 0.0000
10.3 with N = 10. 1.1055 1.1053 0.0003 1.1055 2.7409e-06
1.2242 1.2236 0.0006 1.2242 8.7114e-06
1.3596 1.3585 0.0010 1.3596 1.6792e-05
1.5155 1.5139 0.0016 1.5154 2.7361e-05
1.6962 1.6939 0.0023 1.6961 4.0853e-05
1.9064 1.9032 0.0031 1.9063 5.7764e-05
2.1513 2.1471 0.0041 2.1512 7.8665e-05
2.4366 2.4313 0.0053 2.4365 0.0001
2.7688 2.7620 0.0068 2.7687 0.0001
3.1548 3.1463 0.0085 3.1547 0.0002

We have seen several numerical schemes for solving initial value prob-
lems. There are other methods, or combinations of methods, which aim
to refine the numerical approximations efficiently as if the step size in the
current methods were taken to be much smaller. Some methods extrapolate
solutions to obtain information outside of the solution interval. Others use
one scheme to get a guess to the solution while refining, or correcting, this
to obtain better solutions as the iteration through time proceeds. Such meth-
ods are described in courses in numerical analysis and in the literature. At
this point we will apply these methods to several physics problems before
continuing with analytical solutions.
 numerical solutions of pdes 445

10.2 The Heat Equation

10.2.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.
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)
≈ . (10.25)
∂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)
≈ . (10.26)
∂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
446 partial differential equations

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)
= . (10.27)
(∆x )2

This approximation for u xx is called the central difference approximation of

u xx .
Combining Equation (10.25) with (10.27) 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)] , (10.28)

∆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)] . (10.29)

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 (10.29) becomes
 
ui,j+1 ≈ ui,j + α ui+1,j − 2ui,j + ui−1,j . (10.30)

Equation (10.31) is the finite difference scheme for solving the heat equa-
tion. This equation is represented by the stencil shown in Figure 10.6. The
black circles represent the four terms in the equation, ui,j ui−1,j ui+1,j and
ui,j+1 .
 numerical solutions of pdes 447

t Figure 10.6: This stencil indicates the

four types of terms in the finite differ-
ence scheme in Equation (10.31). The
black circles represent the four terms in
the equation, ui,j ui−1,j ui+1,j and ui,j+1 .

j+1

i−1 i i+1 x

t Figure 10.7: Applying the stencil to the

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

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 10.7, we apply the stencil to generate the solution on the
j = 1 row. This is shown in Figure 10.7.
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 10.8. 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 10.9.
The process continues until we again go as far as we can. This is shown
in Figure 10.10.
448 partial differential equations

Figure 10.8: Continuation of the pro- t

cess provides solutions at the indicated
points.

Figure 10.9: Knowing the values of the t

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

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 10.11.
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.
 numerical solutions of pdes 449

t Figure 10.10: Knowing the values of the

solution at x = a, we can fill in more of
the grid until we stop.

t Figure 10.11: Using boundary conditions

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

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 .
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 old old
+1 − 2ui + ui −1 . (10.31)
450 partial differential equations

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;

% 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);
 numerical solutions of pdes 451

10.3 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 (10.29) and (10.31),

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

(10.32)
.
 
ui,j+1 ≈ ui,j + α ui+1,j − 2ui,j + ui−1,j , (10.33)

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
(10.32) and (10.33), we have not use equality. In order to replace the approx-
imation by an equality, we need t estimate the order of the terms neglected
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 (10.25),

∂u u( x, t + ∆t) − u( x, t)
≈ . (10.34)
∂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 (10.25)
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 (10.25), 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) (10.35)
∂t ∆t
is a first order approximation in ∆t.
452 partial differential equations

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 + . . . . (10.36)
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 + . . . . (10.37)
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!
(10.38)

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 ). (10.39)
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.

10.4 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 , (10.40)
 numerical solutions of pdes 453

α = 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
10.6. 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
= . (10.41)
Tj Xm

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 , (10.42)

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

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 , (10.44)

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. (10.45)

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
454 partial differential equations

ξ = |ξ |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α. (10.46)

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.

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

10.5 Matrix Formulation