Math 124B March 08, 2012 Viktor Grigoryan
17 Finite dierences for the heat equation
In the example considered last time we used the forward dierence for u
t
and the centered dierence for
u
xx
in the heat equation to arrive at the following dierence equation.
u
n+1
j
u
n
j
t
=
u
n
j+1
2u
n
j
+ u
n
j1
(x)
2
. (1)
Denoting s = t/(x)
2
, this lead to the FTCS scheme,
u
n+1
j
= s(u
n
j+1
+ u
n
j1
) + (1 2s)u
n
j
, (2)
which is an explicit numerical scheme. We also saw last time that choosing s = 1 lead to a wild numerical
solution, while the choice s =
1
2
resulted in a somewhat accurate approximation to the solution.
To better understand the eect of s on the outcome of the scheme, let us consider a separated solution
of the dierence equation (1), u
n
j
= X
j
T
n
. Substituting this solution gives
X
j
T
n+1
= s(X
j+1
T
n
+ X
j1
T
n
) + (1 2s)X
j
T
n
.
Dividing both sides of the above equation by u
n
j
, we obtain
T
n+1
T
n
= (1 2s) + s
X
j+1
+ X
j1
X
j
= , (3)
where is independent of both j and n. We then have for the T component
T
n+1
= T
n
, and hence, T
n
=
n
T
0
.
For the scheme (2) to be stable, we need to have || 1, since otherwise it would lead to solutions that
grow exponentially in time.
The equation for the X component is
1 2s + s
X
j+1
+ X
j1
X
j
= .
Plugging in a discretized Fourier mode X
j
= e
ikjx
into this equation, we obtain
= 1 2s + s
e
ik(j+1)x
+ e
ik(j1)x
e
ikjx
, or = 1 2s + s(e
ikx
+ e
ikx
).
Using Eulers formula, we can rewrite the last identity as
= 1 2s + 2s cos kx = 1 2s(1 cos kx).
Since cos kx 1, we see that 1 for all values of s. On the other hand, the condition 1 is
equivalent to
1 1 2s(1 cos kx).
The worst case happens when the frequency k is such that cos kx 1, then the above inequality would
give 1 1 4s, which implies the following necessary condition for stability
t
(x)
2
= s
1
2
. (4)
1
17.1 Boundary conditions
Last time we saw how one uses the discretized initial data to march forward in time using the explicit
scheme (2). In the presence of Dirichlet boundary conditions, the discretized boundary data is also used
when computing the numerical solution.
For the Neumann boundary conditions,
u
x
(0, t) = g(t), u
x
(l, t) = h(t),
the values of the solution on the boundary points of the grid are not immediately available, so one
needs to use a dierence approximation for the conditions in order to march forward in time. Using
the forward or backward dierences to approximate u
x
on the boundary would introduce a local trun-
cation error of order O(x), which is bigger than the truncation error O(x)
2
, thus contaminating
the numerical solution. To have errors of the same order as the dierence equation, we must use the
centered dierences. For these we introduce ghost points u
n
1
and u
n
J+1
in addition to the grid points
u
n
0
, u
n
1
, . . . , u
n
J
. Then the centered dierence approximation for the Neumann conditions will be
g(nt) = g
n
=
u
n
1
u
n
1
2x
, and h(nt) = h
n
=
u
n
J+1
u
n
J1
2x
. (5)
From the above identities we can compute the values u
n
1
and u
n
J+1
at the ghost points at time level n,
which will then be used to compute u
n+1
0
and u
n+1
J
, thus arriving at numerical values at the boundaries.
17.2 The implicit BTCS scheme
Instead of approximating u
t
in the heat equation u
t
= u
xx
by the forward dierence, which resulted in
the dierence equation (1), we can use the backward dierence, which at time level n + 1 will give the
dierence equation
u
n+1
j
u
n
j
t
=
u
n+1
j+1
2u
n+1
j
+ u
n+1
j1
(x)
2
. (6)
Denoting as before s = t/(x)
2
, and separating values at dierent time levels, we can rewrite the
above equation as
su
n+1
j1
+ (1 + 2s)u
n+1
j
su
n+1
j1
= u
n
j
. (7)
For a xed n, the last equation can be thought of as a system for u
n+1
1
, u
n+1
2
, . . . , u
n+1
J
. Thus, (7) gives
an implicit scheme, which corresponds to the template
s
1 + 2s
s
1 + 2s
1
1 + 2s
To nd the numerical solution using this scheme, one needs to solve the system of linear equations (7).
In the presence of Dirichlet boundary conditions, this system can be written in the following vector form
_
_
_
_
_
_
_
1 + 2s s 0 . . 0
s 1 + 2s s 0 . .
0 s 1 + 2s . .
. . . . 0
. . . s
0 . . 0 s 1 + 2s
_
_
_
_
_
_
_
_
_
u
n+1
1
u
n+1
2
.
.
u
n+1
J1
u
n+1
J
_
_
=
_
_
u
n
1
+ su
n+1
0
u
n
2
.
.
u
n
J1
u
n
J
+ su
n+1
J
_
_
.
We can see that the matrix of the coecients is tridiagonal, which can be inverted by ecient algorithms.
Although (7) is implicit, it requires about twice as many arithmetic operations as the scheme (2).
Let us now check for what values of s the BTCS scheme (7) will be stable. Plugging the separated
2
solution u
n
j
=
n
e
ikjx
into the scheme and dividing both sides by u
n
j
, we have
se
ikx
+ (1 + 2s) se
ikx
= 1.
Factoring on the left hand side, and solving for it gives
=
1
1 + 2s s(e
ikx
+ e
ikx
)
=
1
1 + 2s(1 cos kx)
,
so || 1, since 1 cos kx 0. Thus, the implicit scheme (7) is stable for all values of s, i.e.
unconditionally stable.
17.3 The -scheme
The two schemes for the heat equation considered so far have their advantages and disadvantages. On
the one hand we have the FTCS scheme (2), which is explicit, hence easier to implement, but it has
the stability condition t
1
2
(x)
2
. The last fact requires very small mesh size for the time variable,
which leads one to consider more time steps to reach the values at a certain time.
On the other hand, the implicit scheme BTCS (7) requires more arithmetic operations to nd the
values at a certain time step, but it is unconditionally stable, allowing one to chose a larger mesh size for
the time variable. Notice also that the two schemes use dierent sets of points in the computation of u
n+1
j
Combining the two schemes with dierent weights into a single scheme may emphasize some of the
advantages of these schemes, and also be more accurate, since it would use a larger set of points to
compute the same values.
Taking a parameter 0 1, we form the scheme
u
n+1
j
u
n
j
t
= (1 )
u
n
j+1
2u
n
j
+ u
n
j1
(x)
2
+
u
n+1
j+1
2u
n+1
j
+ u
n+1
j1
(x)
2
. (8)
Observe that when = 0 the above scheme is exactly the BTCS scheme (7), while when = 1 the
scheme become the FTCS scheme (2). Clearly for 0 < 1 this scheme is implicit.
Denoting, as before, s = t/(x)
2
, we analyze the stability of the scheme (8) by substituting in the
separated solution u
n
j
=
n
e
ikjx
. This leads to
1 = (1 )s(e
ikx
2 + e
ikx
) + s(e
ikx
2 + e
ikx
).
Solving for , we obtain
=
1 2(1 )s(1 cos kx)
1 + 2s(1 cos kx)
.
Then the stability condition || 1 is equivalent to
1
1 2(1 )s(1 cos kx)
1 + 2s(1 cos kx)
1. (9)
The inequality on the right is always satised, since the numerator of the fraction is always less than or
equal to one, while the denominator is always more than or equal to one. Thus we need to only check
the inequality on the left, which is equivalent to
1 2s(1 cos kx) 1 2s(1 cos kx) + 2s(1 cos kx),
or, after combining like terms,
2 s(4 2)(1 cos kx).
3
The worst case again corresponds to the frequency, for which cos kx 1, resulting in
2 2s(4 2).
The last inequality will be always satised, if
1
2
, since then the right hand side would be nonnegative.
Thus, the -scheme (8) is unconditionally stable for
1
2
1. On the other hand, when 0 <
1
2
,
the stability condition is
t
(x)
2
= s
1
2 4
.
Notice that for = 0 this condition exactly coincides with the stability condition (4) for the FTCS
scheme, which was expected, since the -scheme reduces to the FTCS scheme for = 0.
In the special case =
1
2
the scheme (8) is called the Crank-Nicholson scheme. It can be written as
the average of (2) and (7),
s
2
u
n+1
j1
+ (1 + s)u
n+1
j
s
2
u
n+1
j+1
=
s
2
u
n
j1
+ (1 s)u
n
j
+
s
2
u
n
j1
, (10)
and will have the template
1
2
s
1 + s
1
2
s
1 + s
1
2
s
1 + s
1 s
1 + s
1
2
s
1 + s
The Crank-Nicholson scheme (10) is more accurate than (2) and (7) for small values of t, however, it
is the most computationally involved.
17.4 Conclusion
Using either the forward dierence approximation or the backward dierence approximation for the
time derivative in the heat equation, we obtained the explicit scheme (2) and the implicit scheme (7)
respectively. The stability analysis of these schemes showed that for stability of the rst scheme we need
t
1
2
(x)
2
, while the second scheme is unconditionally stable.
Combining these two schemes with dierent weights into one gave as the implicit -scheme (8), which
is unconditionally stable for
1
2
. In the special case of =
1
2
it is called the Crank-Nicholson scheme,
which is given by (10).
4
