NUMERICAL METHODS FOR PARABOLIC EQUATIONS As a

NUMERICAL METHODS FOR PARABOLIC EQUATIONS
LONG CHEN
As a model problem of general parabolic equations, we shall mainly consider the following heat equation and study corresponding finite difference methods and finite element
methods

in Ω × (0, T ),
 ut − ∆u = f
u = 0
on ∂Ω × (0, T ),
(1)

u(·, 0) = u0
in Ω.
Here u = u(x, t) is a function of spatial variable x ∈ Ω and time variable t ∈ (0, T ).
The Laplace differential operator ∆ is taking with respect to the spatial variable. For the
simplicity of exposition, we consider only homogenous Dirichlet boundary condition and
comment on the adaptation to Neumann and other type of boundary conditions. Comparing
with elliptic equation, we also need to assign the value at time t = 0 which is called
initial condition. The ending time T could be +∞. For parabolic equations, the boundary
∂Ω × (0, T ) ∪ Ω × {t = 0} is called the parabolic boundary. Therefore the initial condition
can be also thought as a boundary condition.
1. BACKGROUND ON HEAT EQUATION
For the homogenous Dirichlet boundary condition without source term, in the steady
state, i.e., ut = 0, we obtain the Laplace equation
∆u = 0 in Ω and u|∂Ω = 0.
So u = 0 no matter what the initial condition is. Indeed the solution will decay to zero
exponentially. Let us consider the simple 1-D problem
ut = uxx in R1 × (0, T ),
u(·, 0) = u0 .
We apply Fourier transfer in space
Z
u
ˆ(k, t) =
u(x, t)e−ikx dx.
R
2
Then u
cx = (−ik)ˆ
u, u
d
ˆ, and ubt = u
ˆt . So we get the following ODE for each
xx = −k u
Fourier coefficient u
ˆ(k, t)
u
ˆt = −k 2 u
ˆ, u
ˆ(·, 0) = u
ˆ0
2
The solution in frequency domain is u
ˆ(k, t) = u
ˆ0 e−k t . We apply inverse Fourier transform back to (x, t) coordinate and get
Z
(x−y)2
1
e− 4t u0 (y) dy.
u(x, t) = √
4πt R
For a general bounded domain Ω, we cannot apply Fourier transform. Instead we can
use the eigenfunctions of A = −∆0 (Laplace operator with zero Dirichlet boundary condition). Since A is SPD, we know there exists an orthogonormal basis formed by eigen2
functions of
function in such bases
PA, i.e., L = span{φ1 , φ2 , . . . , }. We expand the
u(x, t) = k uk (t)φk (x). The heat equation then becomes ukt (t) = λk uk , uk (0) = uk0
1
2
LONG CHEN
and the solution is uk = uk0 e−λk t for k = 1, 2, . . .. Again each component will exponentially decay to zero since the eigenvalue λk of A is positive. And the larger the eigenvalue
is, the faster the decay rate. In practice, however, it is much harder to finding out all eigenvalue and eigenfunctions than solving the heat equation numerically. We will talk about
finite difference and finite element methods for solving the heat equation.
2. F INITE DIFFERENCE METHODS FOR 1-D HEAT EQUATION
In this section, we consider a simple 1-D heat equation
(2)
(3)
ut = uxx + f
in (0, 1) × (0, T ),
u(0) = u(1) = 0, u(x, 0) = u0 (x).
to illustrate the main issues in the numerical methods for solving parabolic equations.
Let Ω = (0, 1) be decomposed into a uniform grid {0 = x0 < x1 < . . . < xN +1 = 1}
with xi = ih, h = 1/N , and time interval (0, T ) be decomposed into {0 = t0 < t1 <
. . . < tM = T with tn = nδt, δt = T /M . The tensor product of these two grids give a
two dimensional rectangular grid for the domain Ω × (0, T ). We now introduce three finite
difference methods by discretizing the equation (2) on grid points.
2.1. Forward Euler method. We shall approximate the function value u(xi , tn ) by Uin
and uxx by second order central difference
n
n
Ui−1
+ Ui+1
− 2Uin
.
2
h
For the time derivative, we use the forward Euler scheme
uxx (xi , tn ) ≈
(4)
ut (xi , tn ) ≈
Uin+1 − Uin
.
δt
Together with the initial condition and the source Fin = f (xi , tn ), we then end with a
system
(5)
(6)
n
U n + Ui+1
− 2Uin
Uin+1 − Uin
= i−1
+ Fin ,
δt
h2
Ui0 = u0 (xi ),
1 ≤ i ≤ N, 1 ≤ n ≤ M
1 ≤ i ≤ N, n = 0.
To write (5) in a compact form, we introduce the parameter λ = δt/h2 and the vector
n t
) . Then (5) can be written as, for n = 0, · · · , M
U n = (U1n , U2n , . . . , UN
U n+1 = AU n + δtF n ,
where



A = I + λ∆h = 


1 − 2λ
λ
...
0
0
λ
1 − 2λ
...
λ
0
0
λ
...
1 − 2λ
λ
0
0
...
λ
1 − 2λ



.


Starting from t = 0, we can evaluate point values at grid points from the initial condition
and thus obtain U 0 . After that, the unknown at next time step is computed by one matrixvector multiplication and vector addition which can be done very efficiently without storing
the matrix. This method is also called time marching.
NUMERICAL METHODS FOR PARABOLIC EQUATIONS
3
Remark 2.1. Because of the homogenous Dirichlet boundary condition, the boundary
index i = 0, N + 1 is not included. For Neumann boundary condition, we do need to
impose equations on these two boundary nodes and introduce ghost points for accurately
discretize the Neumann boundary condition; See Chapter: Finite difference methods for
elliptic equations.
The first issue is on the stability in time. When f = 0, i.e., heat equation without
source, in the continuous level, the solution should exponential decay. In the discrete level,
we have U n+1 = AU n and want to control the magnitude of U in certain norm.
Theorem 2.2. When the time step δt ≤ h2 /2, the forward Euler method is stable in the
maximum norm in the sense that if U n+1 = AU n then
kU n k∞ ≤ kU 0 k∞ ≤ ku0 k∞ .
Proof. By the definition of the norm of a matrix
kAk∞ =
max
i=1,··· ,N
N
X
|aij | = 2λ + |1 − 2λ|.
j=1
If δt ≤ h2 /2, then kAk∞ = 1 and consequently
kU n k∞ ≤ kAk∞ kU n−1 k∞ ≤ kU n−1 k∞ .
Exercise 2.3. Construct an example to show, numerically or theoretically, that if δt >
h2 /2, then
kU n k∞ > kU 0 k∞ .
Theorem 2.2 and Exercise 2.3 imply that to ensure the stability of the forward Euler
method, we have to choose the time step δt in the size of h2 which is very restrictive, say,
h = 10−3 then t = 10−6 . Although in one step, it is efficient, it will be very expensive
to reach the solution at the ending time T by moving at such tiny time step. For time
dependent equations, we shall consider not only the computational cost for one single time
step but also the total time to arrive a certain stopping time.
The second issue is the accuracy. We first consider the consistency error. Define unI =
(u(x1 , tn ), u(x2 , tn ), · · · , u(xN , tn ))t . We denote a general differential operator as L and
its discritization as Lnh . Note that the action Lu is the operator L acting on a continuous
function u while Lnh is applied to vectors such as U n or unI . The consistent error or
so-called truncation error is to pretend we know the exact function values and see what
the error from the approximation of the differential operator. More precisely, for the heat
equation we define Lu = ut − ∆u and Lnh V = (V n+1 − V n )/δt − ∆h /h2 V n and the
truncation error
τ n = Lnh unI − (Lu)nI
as the truncation error. The error is denoted by
E n = U n − unI .
The estimate of the truncation error is straightforward using Taylor expansion.
Lemma 2.4. Suppose u is smooth enough. For the forward Euler method, we have
kτ n k∞ ≤ C1 δt + C2 h2 .
4
LONG CHEN
The convergence or the estimate of the error E n is a consequence of the stability and
consistency. The method can be written as
Lnh U n = F n = fIn = (Lu)nI .
We then obtain the error equation
Lnh (U n − unI ) = (Lu)nI − Lnh unI ,
(7)
which can be simply written as Lnh E n = −τ n .
Theorem 2.5. For the forward Euler method, when δt ≤ h2 /2 and the solution u is smooth
enough, we have
kU n − unI k∞ ≤ Ctn (δt + h2 ).
Proof. We write out the specific error equation for the forward Euler method
E n+1 = AE n − δt τ n .
Consequently
E n = An E 0 − δt
n−1
X
An−l−1 τ l .
l=1
0
Since E = 0, by the stability and consistency, we obtain
kE n k∞ ≤ δt
n−1
X
kτ l k∞ ≤ Cnδt(δt + h2 ) = Ctn (δt + h2 ).
l=1
2.2. Backward Euler method. Now we study backward Euler to remove the strong constrain of the time step for the stability. The method is simply using backward difference to
approximate the time derivative. We list the system below:
(8)
(9)
n
U n + Ui+1
− 2Uin
Uin − Uin−1
= i−1
+ Fin ,
δt
h2
Ui0 = u0 (xi ),
1 ≤ i ≤ N, 1 ≤ n ≤ M
1 ≤ i ≤ N, n = 0.
In the matrix form (8) reads as
(10)
(I − λ∆h )U n = U n−1 + δtF n .
Starting from U 0 , to compute the value at the next time step, we need to solve an algebraic
equation to obtain
U n = (I − λ∆h )−1 (U n−1 + δtF n ).
The inverse of the matrix, which involves the stiffness matrix of Laplacian operator, is not
easy in high dimensions. For 1-D problem, the matrix is tri-diagonal and can be solved
very efficiently by the Thomas algorithm.
The gain is the unconditional stability.
Theorem 2.6. For the backward Euler method without source term, i.e., (I − λ∆h )U n =
U n−1 , we always have the stability
kU n k∞ ≤ kU n−1 k∞ ≤ ku0 k∞ .
NUMERICAL METHODS FOR PARABOLIC EQUATIONS
5
Proof. We shall rewrite the backward Euler scheme as
n
n
(1 + 2λ)Uin = Uin−1 + λUi−1
+ λUi+1
.
Therefore for any 1 ≤ i ≤ N ,
(1 + 2λ)|Uin | ≤ kU n−1 k∞ + 2λkU n k∞ ,
which implies
(1 + 2λ)kU n k∞ ≤ kU n−1 k∞ + 2λkU n k∞ ,
and the desired result then follows.
The truncation error of the backward Euler method can be obtained similarly.
Lemma 2.7. Suppose u is smooth enough. For the backward Euler method, we have
kτ n k∞ ≤ C1 δt + C2 h2 .
We then use stability and consistency to give error analysis of the backward Euler
method.
Theorem 2.8. For backward Euler method, when the solution u is smooth enough, we
have
kU n − unI k∞ ≤ Ctn (δt + h2 ).
Proof. We write the error equation for the backward Euler method as
n
n
(1 + 2λ)Ein = Ein−1 + λEi−1
+ λEi+1
− δt τin .
Similar to the proof of the stability, we obtain
kE n k∞ ≤ kE n−1 k∞ + δtkτ n k∞
Consequently, we obtain
kE n k∞ ≤ δt
n−1
X
kτ l k∞ ≤ Ctn (δt + h2 ).
l=1
Exercise 2.9. Apply backward Euler to the example you constructed in Exercise 2.3 to
show numerically the scheme is stable.
From the error analysis, to have optimal convergent order, we also need δt ≈ h2 which
is still restrictive. Next we shall give a unconditional stable scheme with second order
truncation error O(δt2 + h2 ).
2.3. Crank-Nicolson method. To improve the truncation error, we need to use central
difference for time discritization. We keep the forward discretization as (Uin+1 − Uin )/δt
but now treat it is an approximation of ut (xi , tn+1/2 ). That is we discretize the equation
n+1/2
at (xi , tn+1/2 ). The value Ui
is taken as average of Uin and Uin+1 . We then end with
the scheme: for 1 ≤ i ≤ N, 1 ≤ n ≤ M
n+1
n+1
n+1
n
n
+ Ui+1
− 2Uin
1 Ui−1
1 Ui−1 + Ui+1 − 2Ui
Uin+1 − Uin
n+1/2
=
+
+ Fi
,
δt
2
h2
2
h2
and for 1 ≤ i ≤ N, n = 0
Ui0 = u0 (xi ).
(11)
6
LONG CHEN
In matrix form, (45) can be written as
1
1
(I − λ∆h )U n+1 = (I + λ∆h )U n + F n+1/2 .
2
2
It can be easily verify that this time the truncation error is improved as
Lemma 2.10. Suppose u is smooth enough. For Crank-Nicolson method, we have
kτ n k∞ ≤ C1 δt2 + C2 h2 .
Exercise 2.11. Prove the maximum norm stability of Crank-Nicolson method with assumption on λ. What is the weakest assumption on λ you need?
In the next section, we shall prove Crank-Nicolson method is unconditionally stable in
the l2 norm and thus obtain the following second order convergence; See Exercise 13.
Theorem 2.12. For Crank-Nicolson method, we have
kE n k ≤ Ctn (C1 δt2 + C2 h2 ).
3. VON N EUMANN ANALYSIS
One way to study the L2 stability of the heat equation is through Fourier analysis and
its discrete version. When Ω = R, we can use Fourier analysis. For Ω = (0, 1), we can
use eigenfunctions of Laplacian operator. To study the matrix equation, we establish the
discrete counter part.
Lemma 3.1. If A = diag(b, a, b) be a N × N matrix, then the eigenvalue of A is
λk = a + 2b cos θk ,
k = 1, · · · , N
and its corresponding eigenvector is
√
φk = 2 (sin θk , sin 2θk , · · · , sin N θk ),
where
θk = kθ =
kπ
.
N +1
Proof. It is can be easily verified by the direct calculation.
Remark 3.2. The eigenvectors do not depend on the values a and b!
We then define a scaling of l2 inner product of two vectors as
(U , V )h = h U t V = h
N
X
Ui Vi ,
i=1
which is a mimic of L2 inner product of two functions if we thought Ui and Vi represent
the values of corresponding functions at xi .
N
It is straightforward to verify the eigenvectors
√ {φk }k=1 forms an orthonormal basis of
N
R with respect to (·, ·)h . Indeed the scaling 2 is introduced for the normalization. For
any vector V ∈ RN , we expand it in this new basis
V =
N
X
vˆk φk .
k=1
The discrete Parseval identity is
kV kh = kˆ
v k,
NUMERICAL METHODS FOR PARABOLIC EQUATIONS
7
where v
ˆ = (ˆ
v1 , · · · , vˆN )t is the vector formed by the coefficients. Theˆin the coefficient
indicates this is a mimic of Fourier transform in the discrete level.
Suppose
U n+1 = AU n ,
with A = diag(b, a, b). Then the stability in k · kh is related to the spectrum radius of A.
That is
kU n+1 kh ≤ kAkh kU n kh ,
and since A is symmetric with respect to (·, ·)h ,
kAkh = ρ(A) = max |λk (A)|.
1≤k≤N
1/2
Note that k · kh = h k · k. The stability in the scaled l2 norm is equivalent to that in
k · k norm. The scaling h is chosen to be consistent with the scaling of discrete Laplace
operator.
Now we analyze stability of the three numerical methods we have discussed. Recall
that λ = δ/h2 .
Forward Euler.
U n+1 = AU n ,
A = diag(λ, 1 − 2λ, λ).
Thus
ρ(A) = max |1 − 2λ + 2λ cos θk |.
1≤k≤N
It is easy to verify that
λ ≤ 1/2 =⇒ ρ(A) ≤ 1.
Let us write A = AN , the uniform stable with respect to the dimension N implies λ ≤ 1/2
sup ρ(AN ) ≤ 1 =⇒ λ ≤ 1/2.
N
Therefore we obtain the same condition as that for the maximum norm stability.
Backward Euler.
U n+1 = A−1 U n ,
Thus
ρ(A−1 ) =
A = diag(−λ, 1 + 2λ, −λ).
1
1
= max
,
ρ(A) 1≤k≤N |1 + 2λ − 2λ cos θk |
and
ρ(A−1 ) ≤ 1 for any λ > 0.
Therefore backward Euler is also unconditionally stable in l2 norm.
Crank-Nicolson.
U n+1 = A−1 BU n , where
1
A = diag(−λ, 2 + 2λ, −λ),
2
1
B = diag(λ, 2 − 2λ, λ)
2
Since A and B have the same eigenvectors, we have
ρ(A−1 B) = max
1≤k≤N
λk (B)
|1 − λ + λ cos θk |
= max
,
λk (A) 1≤k≤N |1 + λ − λ cos θk |
8
LONG CHEN
and
ρ(A−1 B) ≤ 1
for any λ > 0.
Therefore we proved Crank-Nicolson method is unconditionally stable in l2 norm.
Question 3.3. Is Crank-Nicolson method unconditional stable in the maximum norm?
For general schemes, one can write out the matrix form and plug in the coefficients to
do the stability analysis. This methodology is known as von Neumann analysis. Formally
we replace
(12)
Ujn ← ρn eiθj
in the scheme and obtain a formula of the amplification factor ρ = ρ(θ). The stability is
obtained by considering
max ρ(θ).
0≤θ≤π
To fascinate the calculation, one can factor out ρn eiθj and consider only the difference of
indices.
In (12), a more precise notation is θ = θk which represents the frequency of eigenvect
tors φk = eiθk (1:N ) which can take care of more general boundary conditions. The U n
is expanded in this basis with a coefficient vˆkn = ρnk . We skip the index k and take the
maximum of θ over [0, π] while hπ ≤ θk ≤ (1 − h)π.
As an example, the readers are encouraged to do the following exercise. Use as an
example not exercise.
Exercise 3.4 (The θ method). For θ ∈ [0, 1], we use the scheme: for 1 ≤ i ≤ N, 1 ≤ n ≤
M
n+1
n+1
n
Ui−1
+ Ui+1
− 2Uin+1
U n + Ui+1
− 2Uin
Uin+1 − Uin
= (1 − θ) i−1
+
θ
,
δt
h2
h2
and for 1 ≤ i ≤ N, n = 0
Ui0 = u0 (xi ).
(13)
Give a complete error analysis (stability, consistency, and convergence) of the θ method.
Note that θ = 0 is forward Euler, θ = 1 is backward Euler, and θ = 1/2 is Crank-Nicolson
method.
Exercise 3.5 (Leap-frog method). The scheme is an explicit version of Crank-Nicolson.
For 1 ≤ i ≤ N, 1 ≤ n ≤ M
(14)
n
U n + Ui+1
− 2Uin
Uin+1 − Uin−1
= i−1
.
2δt
h2
The computation of U n+1 can be formulated as linear combination of U n and U n−1 .
Namely it is an explicit scheme. To start the computation, we need two ‘initial’ values.
One is the real initial condition
Ui0 = u0 (xi ),
1 ≤ i ≤ N.
The other U 1 , so-called computational initial condition, can be computed using one step
forward or backward Euler method.
Give a complete error analysis (stability, consistency, and convergence) of the leap-frog
method.
NUMERICAL METHODS FOR PARABOLIC EQUATIONS
9
4. VARIATIONAL FORMULATION AND ENERGY ESTIMATE
The first try is to multiply a test function v ∈ H01 (Ω) and apply the integration by part.
We obtain a variational formulation of the heat equation (1): given an f ∈ L2 (Ω) × (0, T ],
for any t > 0, find u(·, t) ∈ H01 (Ω), ut ∈ L2 (Ω) such that
(15)
(ut , v) + a(u, v) = (f, v), for all v ∈ H01 (Ω).
We then refine the weak formulation (15). The right hand side could be generalized to
f ∈ H −1 (Ω). Since ∆ map H01 (Ω) to H −1 (Ω), we can treat ut (·, t) ∈ H −1 (Ω) for a
fixed t. We then introduce the Sobolev space for the time dependent functions
!1/q
Z
T
Lq (0, T ; W k,p (Ω)) := {u(x, t) | kukLq (0,T ;W k,p (Ω)) :=
0
ku(·, t)kqk,p dt
< ∞}.
Our refined weak formulation will be: given f ∈ L2 (0, T ; H −1 (Ω)) and u0 ∈ H01 (Ω),
find u ∈ L2 (0, T ; H01 (Ω)) and ut ∈ L2 (0, T ; H −1 (Ω)) such that
hut , vi + a(u, v) = hf, vi,
∀v ∈ H01 (Ω), and a.e.t ∈ (0, T )
(16)
u(·, 0) = u0
We assume the equation (16) is well posed. For the existence and uniqueness of the solution
(u, ut ), we refer to [1].
Remark 4.1. The topology for the time variable should be also treat in L2 sense. But in
(15) and (16) we still pose the equation point-wise in time. In particular, one has to justify
the point value u(·, 0) does make sense for a L2 type function. This can be justified by the
regularity theory.
To easy the stability analysis, we treat t as a parameter and the function u = u(x, t) as
a mapping
u : [0, T ] → H01 (Ω),
defined as
u(t)(x) := u(x, t)
(x ∈ Ω, 0 ≤ t ≤ T ).
With a slight abuse of the notation, here we still use u(t) to denote the map. The norm
ku(t)k or ku(t)k1 is taken with respect to the spatial variable.
We then introduce the differential operator
L : L2 (0, T ; H01 (Ω)) → L2 (0, T ; H −1 (Ω)) × H01 (Ω)
as
(Lu)(·, t) = ∂t u − ∆u in H −1 (Ω), for t ∈ (0, T ] a.e.
(Lu)(·, 0) = u(·, 0).
Then the equation (16) can be written as
Lu = (f, u0 ).
Here we explicitly include the initial condition. Note that the spatial boundary condition is
build into the space H01 (Ω).
We shall prove the stability result of L. It is known as energy estimate following [3].
10
LONG CHEN
Theorem 4.2. Suppose (u, ut ) is the solution of (16) and ut ∈ L2 (0, T ; L2 (Ω)), then for
t ∈ (0, T ] a.e.
Z t
(17)
ku(t)k ≤ ku0 k +
kf (s)k ds
0
(18)
t
Z
2
|u(s)|21
ku(t)k +
2
|u(t)|21 +
Z
kf (s)k2−1 ds,
ds ≤ ku0 k +
0
(19)
t
Z
0
t
kut (s)k2 ds ≤ |u0 |21 +
Z
0
t
kf (s)k2 ds.
0
Proof. The solution is defined via the action of all test function. The art of the energy
estimate is to choose a correct test function to extract desirable information.
We first choose v = u to obtain
(ut , u) + a(u, u) = (f, u).
We manipulate these three terms as:
Z
1 2
1 d
d
• (ut , u) =
(u )t =
kuk2 = kuk kuk;
2
2
dt
dt
Ω
• a(u, u) = |u|21 ;
1
• |(f, u)| ≤ kf kkuk or |(f, u)| ≤ kf k−1 |u|1 ≤ (kf k2−1 + |u|21 ).
2
The inequality (17) is an easy consequence of the following inequality
kuk
d
kuk ≤ kf kkuk.
dt
From
1 d
1
kuk2 + |u|21 ≤ (kf k2−1 + |u|21 ),
2 dt
2
we get
d
kuk2 + |u|21 ≤ kf k2−1 .
dt
Integrating over (0, t), we obtain (18).
The inequality (19) can be proved similarly by choosing v = ut and left as an exercise.
From (18), we can obtain the ‘right’ stability for the operator
L : L2 (0, T ; H01 (Ω)) → L2 (0, T ; H −1 (Ω)) × H01 (Ω)
as
kukL2 (0,T ;H01 (Ω)) ≤ kuk0 + kf kL2 (0,T ;H −1 (Ω)) .
Since the equation is posed a.e for t, we could obtain maximum type estimate in time. For
example, (17) can be formulated as
kukL∞ (0,T ;L2 (Ω)) ≤ ku0 k + kf kL1 (0,T ;L2 (Ω)) ,
and (19) implies
kukL∞ (0,T ;H01 (Ω)) ≤ |u0 |1 + kf kL2 (0,T ;L2 (Ω)) .
Exercise 4.3. Prove the stability result (19).
NUMERICAL METHODS FOR PARABOLIC EQUATIONS
11
Exercise 4.4. Prove the energy estimate
(20)
kuk ≤ e−λt ku0 k +
Z
t
e−λ(t−s)) kf k−1 ds,
0
where λ = λmin (−∆) > 0. The estimate (20) shows that the effect of the initial data is
exponential decay.
5. F INITE ELEMENT METHODS : SEMI - DISCRETIZATION IN SPACE
Let {Th , h → 0} be a quasi-uniform family of triangulations of Ω. The semi-discretized
finite element method is: given f ∈ V0h × (0, T ], u0,h ∈ Vh , find uh ∈ L2 (0, T ; Vh ) such
that
(∂t uh , vh ) + a(uh , vh ) = hf, vh i,
∀vh ∈ Vh , t ∈ R+ .
(21)
uh (·, 0) = u0,h
The scheme (21) is called semi-discretization since uh is still a continuous function of
t. The initial condition u0 is approximated by u0,h ∈ Vh and the choice of u0,h is not
unique.
PN
We can expand uh = i=1 ui (t)ϕi (x), where ϕi is the standard hat basis at the vertex
xi for i = 1, · · · , N , the number of interior nodes, and the corresponding coefficient ui (t)
now is a function of time t. The solution uh can be computed by solving an ODE system
u˙ + Au = f ,
where u = (u1 , · · · , uN )t , A is the stiffness matrix, and f = (f1 , · · · , fN )t .
We shall apply our abstract error analysis to estimate the error u − uh in certain norm.
The setting is
R
1/2
T
• X = L2 (0, T ; H01 (Ω)), and kukX = 0 |u(t)|21 dt
R
1/2
T
• Y = L2 (0, T ; H −1 (Ω), and kf kY = 0 kf (t)k2−1 dt
R
1/2
T
• Xh = L2 (0, T ; Vh ), and kuh kXh = 0 |uh (t)|21 dt
R
1/2
T
• Yh = L2 (0, T ; V0h ), and kfh kYh = 0 kfh (t)k2−1,h dt
. Recall the dual
norm
hfh , vh i
, for fh ∈ V0h .
kfh k−1,h = sup
vh ∈Vh |vh |1
• Ih = Rh (t) : H01 (Ω) → Vh is the Ritz-Galerkin projection, i.e., Rh u ∈ Vh such
that
a(Rh u, vh ) = a(u, vh ), ∀vh ∈ Vh .
• Πh = Qh (t) : H −1 (Ω) → V0h is the projection
hQh f, vh i = hf, vh i,
∀vh ∈ Vh .
• Ph : Xh → X is the natural inclusion
• L : X → Y × H01 (Ω) is Lu = ∂t u − ∆u, Lu(·, 0) = u(·, 0), and Lh = L|Xh :
Xh → Yh × Vh . The equation we are solving is
Lh uh = Qh f,
uh (·, 0) = u0,h .
t ∈ (0, T ] a.e. in V0h ,
12
LONG CHEN
5.1. Stability. Adapt the proof of the energy estimate for L, we will obtain similar stability results for Lh .
Theorem 5.1. Suppose uh satisfy Lh uh = fh , uh (·, 0) = u0,h , then
Z t
(22)
kfh (s)k ds
kuh (t)k ≤ ku0,h k +
0
Z t
Z t
(23)
kfh (s)k2−1,h ds,
|uh (s)|21 ds ≤ ku0,h k2 +
kuh (t)k2 +
0
0
Z t
Z t
2
2
2
(24)
|uh (t)|1 +
k∂t uh (s)k ds ≤ |u0,h |1 +
kfh (s)k2 ds.
0
0
Note that in the energy estimate (23), the k · k−1 is replaced by a weaker norm k · k−1,h
since we can apply the inequality
hfh , uh i ≤ kfh k−1,h |uh |1 .
The weaker norm k · k−1,h can be estimated by
kf k−1,h ≤ kf k−1 ≤ Ckf k.
5.2. Consistency. Recall that the consistency error is defined as Lh Ih u−Lh uh = Lh Ih u−
Πh Lu. The choice of Ih = Rh simplifies the consistency error analysis.
Lemma 5.2. For the semi-discretization, we have the error equation
Lh (Rh u − uh ) = ∂t (Rh u − u), t > 0, in V0h ,
(25)
(Rh u − uh )(·, 0) = Rh u0 − u0,h .
(26)
Proof. Let A = −∆. By our definition of consistency, the error equation is: for t > 0
Lh (Ih u − uh ) = Lh Ih u − Πh Lu = Lh Rh u − Qh Lu = ∂t (Rh u − u) + (ARh u − Qh Au).
The desired result then follows by noting that ARh = Qh A in V0h .
The error equation (25) holds in V0h which is a weak topology. The motivation we
choose Ih = Rh is that in this weak topology
hARh u, vh i = (∇Rh u, ∇vh ) = (∇u, ∇v) = hAu, vh i = hQh Au, vh i,
or in operator form
ARh = Qh A.
This technique is firstly proposed by Wheeler [4].
Apply the stability to the error equation, we obtain the following estimate on discrete
error.
Theorem 5.3. The solution uh of (21) satisfy the following error estimate
Z t
(27)
kRh u − uh k ≤ ku0,h − Rh u0 k +
k∂t u − ∂t Rh uk ds
0
2
t
Z
|(uh −Rh u)|21
(28) kRh u−uh k +
2
(29)
Z
k∂t u−∂t Rh uk2−1,h ds,
0
t
2
k∂t (Rh u−uh )k ds ≤
0
t
ds ≤ ku0,h −Rh u0 k +
0
|Rh u−uh |21 +
Z
|u0,h −Rh u0 |21 +
Z
0
t
k∂t u−∂t Rh uk2 ds.
NUMERICAL METHODS FOR PARABOLIC EQUATIONS
13
We then estimate the two terms in these error estimates. The first issue is on the choice
of u0,h . An optimal one is obviously u0,h = Rh u0 so that no error from the the initial
condition. However, this choice requires the inversion of a stiffness matrix which is not
cheap. A simple choice would be the nodal interpolation, i.e. u0,h (xi ) = u0 (xi ) or any
other choice with optimal approximation property
ku0,h − Rh u0 k ≤ ku0 − u0,h k + ku0 − Rh u0 k . h2 ku0 k2 ,
(30)
and similarly
|u0,h − Rh u0 |1 . hku0 k2 .
According to (20), the affect of the initial boundary error will be exponentially decay to
zero as t gros. So in practice, we can choose the simple nodal interpolation.
On the estimate of the second term, we assume ut ∈ H 2 (Ω) to get
k∂t u − ∂t Rh uk = k(I − Rh )ut k . h2 kut k2 .
(31)
When the H 2 regularity result holds for ∆ operator (for example, the domain is smooth or
convex), the negative norm is estimated as
k∂t u − ∂t Rh uk−1,h ≤ k∂t u − ∂t Rh uk−1 ≤ Ck∂t u − ∂t Rh uk . h2 kut k2 .
5.3. Convergence. The convergence of the discrete error comes from the stability and
consistency. For functions in finite element space, the H 1 semi-norm can be estimate by a
simple inverse inequality
|Rh u − uh |1 ≤ Ch−1 kRh u − uh k.
Theorem 5.4. Suppose the solution u to (16) satisfying ut ∈ L2 (0, T ; H 2 (Ω)) and the H 2
regularity holds. Let uh be the solution of (21) with u0,h satisfying (30). We then have
Z t
−1
2
(32)
h |Rh u − uh |1 + kRh u − uh k ≤ Ch ku0 k2 +
kut k2 ds .
0
To estimate the true error u − uh , we need the approximation error estimate of the
projection Rh
h−1 ku − Rh uk + |u − Rh u|1 ≤ Chkuk2 .
Theorem 5.5. Suppose the solution u to (16) satisfying ut ∈ L2 (0, T ; H 2 (Ω)). Then
the solution uh of (21) with u0,h having optimal approximation property (30) satisfy the
following optimal order error estimate:
Z t
(33)
h−1 |u − uh |1 + ku − uh k ≤ Ch2 ku0 k2 +
kut k2 ds .
0
5.4. *Superconvergence and maximum norm estimate. Recall that the error eh = Rh u−
uh satisfies the evolution equation (25) with eh (0) = 0
∂t eh + Ah eh = Qh ∂t (Rh u − u).
(34)
Therefore
Z
eh (t) =
t
e−Ah (t−s) τh ds,
with τh = Qh ∂t (Rh u − u).
0
Due to the smoothing effect of the semi-group e−Ah t , we have the following estimate.
Here we follow the work by Garcia-Archilla and Titi [2].
14
LONG CHEN
Lemma 5.6 (Smoothing property of heat kernel). For τh ∈ Vh , we have
Z t
−Ah (t−s)
e
(35)
max A
τ
ds
h h ≤ C| log h| max kτh k.
0≤t≤T
0≤t≤T
0
Proof. Let λm and λM be the minimal and maximal eigenvalue of Ah . Then it is easy to
check

λM e−λM (t−s) if (t − s) ≤ λ−1

M ,
−Ah (t−s) Ah ≤ (t − s)−1
if λ−1
≤
(t
−
s) ≤ λ−1
e
m ,
M


−λm (t−s)
−1
λm e
if (t − s) ≥ λm .
Note that λm = O(1) and λM ≤ Ch−2 . We get
Z t
−Ah (t−s)
e
max A
τ
ds
h h ≤ C| log h| max kτh k.
0≤t≤T
0≤t≤T
0
Theorem 5.7 (Superconvergence in H 1 -norm). Suppose the solution u to (16) satisfying
ut ∈ L∞ (0, T ; H 2 (Ω)). Let uh be the solution of (21) with u0,h = Rh u0 . Then
max |Rh u − uh |1 ≤ C| log h| h2 max kut k2 .
(36)
0≤t≤T
2
0≤t≤T
2
When ut ∈ L (0, T ; H (Ω)), then
(37)
|Rh u − uh |1 ≤ Ch2
Z
t
1/2
.
kut k22 ds
0
Proof. We multiply A1/2 to (34) and apply Lemma 5.6 to get
Z
T
1/2
−Ah (T −s) 1/2
|eh (T )|1 = kAh eh (T )k = e
Ah τh ds
0
Z
T
−1/2
≤
e−Ah (T −s) Ah ds max kAh τh k
0
0≤t≤T
≤ C| log h|h2 max kut k2 .
0≤t≤T
−1/2
In the last step, we have used the fact kAh τh k = kτh k−1 ≤ kτh k.
To get (37), we use the energy estimate (29).
Since the optimal convergent rate for |u − Rh u|1 or |u − uh |1 is only first order, the
estimate (36) and (37)) are called superconvergence.
To control the maximum norm, we use the discrete embedding result (for 2-D only)
kRh u − uh k∞ ≤ C| log h||Rh u − uh |1 ,
and the error estimate of Rh in the maximum norm
ku − Rh uk∞ ≤ C| log h|h2 kuk2,∞ ,
to obtain the following result.
Theorem 5.8 (Maximum norm estimate for linear element in two dimensions). Suppose
the solution u to (16) satisfying u ∈ L∞ (0, T ; W 2,∞ ) and ut ∈ L2 (0, T ; W 2,∞ ) or ut ∈
NUMERICAL METHODS FOR PARABOLIC EQUATIONS
15
L∞ (0, T ; W 2,∞ (Ω)). Let uh be the solution of (21) with u0,h = Rh u0 . Then in two
dimensions
"
Z t
1/2 #
2
2
k(u − uh )(t)k∞ ≤ C| log h|h kuk2,∞ +
kut k2 ds
(38)
,
0
(39)
2 2
max k(u − uh )(t)k∞ ≤ C| log h| h max [ku(t)k2,∞ + kut (t)k2 ] .
0≤t≤T
0≤t≤T
For high order elements, we could get superconvergence in L2 norm too. Let us define
the order of the polynomial as the degree plus 1, which is the optimal order when measuring
the approximation property in Lp norm. For example, the order of the linear polynomial is
2. When the order of polynomial r is bigger than 3 (,i.e., quadratic and above polynomial),
we can prove a stronger estimate for the negative norm
ku − Rh uk−1 ≤ Chr+1 kukr .
(40)
Using the technique in Lemma 5.6 and Theorem 5.7, we have the following estimate on
the L2 norm.
Theorem 5.9 (Superconvergence in L2 -norm). Suppose the solution u to (16) satisfying
ut ∈ L∞ (0, T ; H 2 (Ω)). Let uh be the solution of (21) with u0,h = Rh u0 . Then
(41)
max k(Rh u − uh )(t)k ≤ C| log h|hr+1 max k(ut )(t)k2r .
0≤t≤T
0≤t≤T
When ut ∈ L2 (0, T ; H 2 (Ω)), then
(42)
kRh u − uh k ≤ Chr+1
Z
t
1/2
.
kut k2r ds
0
Proof. From (34), we have
Z
T
−Ah (T −s)
keh (T )k = e
τh ds
0
Z
T
−1/2
−Ah (T −s) −1/2
e
Ah
ds max kAh τh k
≤
0≤t≤T
0
≤ C| log h|kτh k−1 .
In the second step, we have used the estimate
Z
T
−1/2
e−Ah (T −s) Ah
ds ≤ C| log h|,
0
−1/2
which can be proved by the estimate ke−Ah (T −s) Ah k ≤ (t − s)−1
To prove (42), we simply use the energy estimate (28) and (40).
In 1-D, using the inverse inequality kRh u − uh k∞ ≤ h−1/2 kRh u − uh k, we could
obtain the superconvergence in the maximum norm
Z t
1/2
(43)
kRh u − uh k∞ ≤ Chr+1/2
kut k2r ds
.
0
Again using the inverse inequality kuh − Rh uk∞ ≤ Ch−1 kuh − Rh uk,, the superconvergence of L2 norm, and the maximum norm estimate of Ritz-Galerkin projection, for
r ≥ 3 ku − Rh uk∞ ≤ Chr kukr,∞ , we can improve the maximum norm error estimate.
16
LONG CHEN
Theorem 5.10 (Maximum norm estimate for high order elements in two dimensions). Suppose the solution u to (16) satisfying u ∈ L∞ (0, T ; W 2,∞ ) and ut ∈ L2 (0, T ; W 2,∞ ) or
ut ∈ L∞ (0, T ; W 2,∞ (Ω)). Let uh be the solution of (21) with u0,h = Rh u0 . Then in two
dimensions and for r ≥ 3
"
Z t
1/2 #
kut k2r ds
(44)
k(u − uh )(t)k∞ ≤ Chr kukr,∞ +
.
0
6. F INITE ELEMENT METHODS : SEMI - DISCRETIZATION IN TIME
We consider the semi-discretization in time. We first discretize the time interval (0, T )
by a uniform grid with size δt = T /N and denoted by tn = nδt for n = 0, . . . N . A
continuous function in time will be interpolated into a vector by (I n f )(·, tn ) = f (·, tn ).
Recall that A = −∆ : H01 → H −1 . We list three schemes in operator form.
• Forward Euler. u0 = u0
un − un−1
+ Aun−1 = f n−1 .
δt
• Backward Euler. u0 = u0
un − un−1
+ Aun = f n .
δt
• Crank-Nicolson. u0 = u0
un − un−1
+ A(un + un−1 )/2 = f n−1/2 .
δt
Note that these equations hold in H −1 sense. Taking Crank-Nicolson as an example, the
equation reads
1
1 n
(u − un−1 , v) + (∇un + ∇un−1 , ∇v) = (f n−1/2 , v), for all v ∈ H01 .
δt
2
We now study the stability of these schemes. We rewrite the Backward Euler method as
(45)
(I + δtA)un = un−1 + δtf n .
Since A is SPD, λmin (I + δtA) ≥ 1 and consequently, λmax ((I + δtA)−1 ) ≤ 1. This
implies the L2 stability
kun k ≤ kun−1 k + δtkf n k ≤ ku0 k +
(46)
n
X
δtkf k k.
k=1
The stability (47) is the discrete counter part of (17): discretize the integral
a Riemann sum.
Similarly one can derive the L2 stability for the C-N scheme
kun k ≤ ku0 k +
(47)
n
X
R tn
0
kf k ds by
δtkf k−1/2 k.
k=1
The integral
R tn
0
kf k ds is approximated by the middle point rule.
Remark 6.1. For C-N method, the right hand side can be also chosen as (f n + f n−1 )/2.
It corresponds to the trapezoid rule of the integral. For nonlinear problem A(u), it can be
A((un + un−1 )/2) or (A(un ) + A(un−1 ))/2. The choice is problem dependent.
NUMERICAL METHODS FOR PARABOLIC EQUATIONS
17
The energy method can be adapted to the semi-discretization in time too. We chose
v = (un + un−1 )/2 in (45) to get
1 n 2 1 n−1 2
ku k − ku
k + δt|un−1/2 |21 = δt(f n−1/2 , un−1/2 ),
2
2
which implies the counter part of (18)
n
n
X
X
(48)
kun k2 +
δt|uk−1/2 |21 ≤ ku0 k2 +
δtkf k−1/2 k2−1 .
k=1
k=1
Exercise 6.2. Study the stability of Forward Euler method.
We then study the convergence. We use C-N as a typical example. We apply the discrete
operator Ln to the error I n u − un
Ln (I n u − un ) = Ln I n u − I n Lu =
u(·, tn ) − u(·, tn−1 )
− ∂t u(·, tn−1/2 ).
δt
The consistency error is
u(·, tn ) − u(·, tn−1 )
− ∂t u(·, tn−1/2 ) ≤ C(δt)2 .
(49)
δt
By the stability result, we then get
kI n u − un k ≤ Ctn δt2 .
From the consistency error estimate (49), one can easily see the Backward Euler is only
first order in time. And high order schemes can be constructed by using a high order
difference for the time derivative.
7. F INITE E LEMENT M ETHOD OF PARABOLIC EQUATION : FULL DISCRETIZATION
The full discretization is simply a combination of discretization in space and time. We
consider the following three schemes:
Forward Euler. u0h = u0,h
(50)
(
unh − un−1
h
, vh ) + a(uhn−1 , vh ) = hfhn−1 , vh i,
δt
∀vh ∈ Vh , 1 ≤ n ≤ N
Backward Euler. u0h = u0,h
(51)
(
unh − un−1
h
, vh ) + a(unh , vh ) = hfhn , vh i,
δt
∀vh ∈ Vh , 1 ≤ n ≤ N
Crank-Nicolson. u0h = u0,h
unh − un−1
1
n−1/2
h
, vh )+a( (unh +un−1
), vh ) = hfh
, vh i, ∀vh ∈ Vh , 1 ≤ n ≤ N
h
δt
2
We write the semi-discretization and the full discretization as
(52) (
L → Lh → Lnh .
The discrete error is Rh u(tn ) − unh . The setting is
1/2
R
T
• X = L2 (0, T ; H01 (Ω)), and kukX = 0 |u(t)|21 dt
R
1/2
T
• Y = L2 (0, T ; H −1 (Ω), and kf kY = 0 kf (t)k2−1 dt
P
1/2
N
k 2
• Xhδt = Vh × Vh × · · · × Vh , and kuh kXhδt =
k=0 δt|uh |1
18
LONG CHEN
• Yhδt = V0h × V0h × · · · × V0h , and kfh kYhδt =
n
2
P
N
k=0
δtkfhk k2−1
1/2
.
(0, T ; H01 (Ω))
• I Rh : L
→ Vh × Vh × · · · × Vh is the composition of the nodal
interpolation in time and Ritz projection in space
I n Rh u = Rh I n u = (Rh u)(tn ).
• I n Qh : L2 (0, T ; H −1 (Ω)) → Vh × V0h × · · · × V0h is the composition of the nodal
interpolation in time and L2 projection in space
I n Qh f = Qh I n f = (Qh f )(tn ).
• Lnh : Xhδt → Yhδt × Vh is u0h = u0,h and for n ≥ 1

+ n
n

Forward Euler,
Dt uh − ∆uh
n n
− n
Lh uh = Dt uh − ∆unh
Backward Euler,

 − n
n−1
n
Dt uh − ∆(uh + uh )/2 Crank-Nicolson.
where Dt+ is the forward difference and Dt− the backward difference in time. The
solution unh solves the full discrete equation:
Lnh unh = I n+s Qh f,
in V0h , for n ≥ 1,
u0h = u0,h ,
where s = 0 for backward Euler and forward Euler, and s = −1/2 C-N method.
7.1. Stability. We write the discretization using operator form and the stability in L2 norm
is very natural. Let A = −∆|Vh .
Forward Euler.
(53)
unh = (I − δtA)un−1
+ δtfhn−1 .
h
Backward Euler.
(54)
unh = (I + δtA)−1 (un−1
+ δtfhn ).
h
Crank-Nicolson.
(55)
1
1
n−1/2
unh = (I + δtA)−1 (I − δtA)un−1
+
δtf
.
h
h
2
2
Since A = −∆ is symmetric in the L2 inner product, to obtain the stability in L2 norm,
we only need to study the spectral radius of these operators.
Forward Euler.
ρ(I − δtA) = |1 − δtλmax (A)| ≤ 1,
provided
δt ≤
2
λmax (A)
.
Note that λmax (A) = O(h−2 ). We need the time step is in the size of h2 to make the
forward Euler stable.
Backward Euler. For any δt > 0, since A is SPD, i.e., λmin (A) > 0
ρ((I + δtA)−1 ) = (1 + δtλmin (A))−1 ≤ 1.
NUMERICAL METHODS FOR PARABOLIC EQUATIONS
19
Crank-Nicolson. For any δt > 0,
|1 − 12 δt λk (A)|
1
1
ρ((I + δtA)−1 (I − δt A)) = max
≤ 1,
1≤k≤N 1 + 1 δt λk (A)
2
2
2
We summarize the stability as the following theorem which is a discrete version of (17).
Theorem 7.1. For forward Euler, when δt < 1/λmax (A)
kunh k ≤ ku0h k +
(56)
n−1
X
δtkfhk k.
k=0
For backward Euler
kunh k ≤ ku0h k +
(57)
n−1
X
δtkfhk+1 k.
k=0
For Crank-Nicolson method
kunh k ≤ ku0h k +
(58)
n−1
X
k+1/2
δtkfh
k.
k=0
It is interesting to note that the last term in the stability result (56), (57), and (58) are
different Riemann sums of the integral in the continuous level.
Exercise 7.2. Derive the discrete version of (18) and (19).
7.2. Consistency. The error equation will be
Lnh I n Rh u − I n Qh Lu = Dtn I n Rh u − I n ∂t u = (Dtn I n u − I n ∂t u) + Dtn I n (Rh u − u).
Forward Euler.
k(Dtn I n u − I n ∂t u)k = k
u(tn+1 ) − u(tn )
− ∂t u(tn )k ≤ Cδtk∂tt uk,
δt
and
kDtn I n (Rh u
−1
tn+1
Z
− u)k = kδt
−1 2
(Rh ut − ut ) dsk ≤ δt
Z
tn+1
kut k2 ds.
h
tn
tn
Backward Euler.
k(Dtn I n u − I n ∂t u)k = k
u(tn ) − u(tn−1 )
− ∂t u(tn )k ≤ Cδtk∂tt uk
δt
and
kDtn I n (Rh u − u)k = kδt−1
Z
tn
tn−1
(Rh ut − ut ) dsk ≤ δt−1 h2
Z
tn
kut k2 ds.
tn−1
20
LONG CHEN
Crank-Nicolson. For Crank-Nicolson, it is important to note that we are solving
Lnh unh = I n−1/2 Qh f.
So the error equation is
Rh u(tn ) − Rh u(tn−1 )
− ∂t u(tn−1/2 )
δt
1
+ ARh (u(tn ) + u(tn−1 )) − Qh Au(tn−1/2 )
2
= Dt I n (Rh u − u)
Lnh I n Rh u − I n−1/2 Qh Lu =
+ Dt I n u − ∂t u(tn−1/2 )
1
+ Qh A (u(tn ) + u(tn−1 )) − u(tn−1/2 )
2
= I1 + I2 + I3 .
We then estimate each term as
kI1 k = δt−1 h2
Z
tn
kut k2 ds,
tn−1
2
kI2 k ≤ Ck∂ttt ukδt ,
kI3 k ≤ Ck∆∂tt ukδt2 .
7.3. Convergence. The convergence is a consequence of stability, consistency, and the
approximation. For simplicity, we list the result for Crank-Nicolson and assume u0,h =
Rh u0 .
Theorem 7.3. Let u, unh be the solution of (16) and (52), respectively with u0,h = Rh u0 .
Then for n ≥ 0,
Z tn
Z tn
n
2
2
2
ku(tn ) − uh k ≤ Ch
kut k ds + Cδt
(k∂ttt uk + k∂tt ∆uk).
0
0
R EFERENCES
[1] L. C. Evans. Partial Differential Equations. American Mathematical Society, 1998.
[2] B. Garc´Ia-Archilla and E. S. Titi. Postprocessing the galerkin method: The finite-element case. SIAM J.
Numer. Anal., 37(2):470–499, 2000.
[3] V. Thom´ee. Galerkin finite element methods for parabolic problems, volume 25 of Springer Series in Computational Mathematics. Springer-Verlag, Berlin, second edition, 2006.
[4] M. F. Wheeler. A priori L2 error estimates for Galerkin approximations to parabolic partial differential equations. SIAM J. Numer. Anal., 10:723–759, 1973.

Download Report