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 ∈ Ω ⊂ Rn and time variable t ∈ (0, T ). The ending time T could be +∞. The Laplace 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. Besides the boundary condition on ∂Ω, we also need to assign the function value at time t = 0 which is called initial condition. 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 simplest 1-D problem ut = uxx in R1 × (0, T ),

u(·, 0) = u0 .

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 the frequency domain is u ˆ(k, t) = u ˆ0 e−k t . We apply the 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 that there exists an orthogonormal basis formed by 1

2

LONG CHEN

eigenfunctions A, i.e., L2 = span{φ1 , φ2 , . . . , }. We expand the function in such bases P of k u(x, t) = k u (t)φk (x). The heat equation then becomes ukt (t) = λk uk , uk (0) = uk0 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, the faster the decay rate. This spectral analysis is mainly for theoretical purpose and the numerical application is restricted to special domains. In practice, for domains of complex geometry, it is much harder to finding out all eigenvalue and eigenfunctions than solving the heat equation numerically. In the following sections we will talk about finite difference and finite element methods. 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 U n + Ui+1 − 2Uin . uxx (xi , tn ) ≈ i−1 h2 For the time derivative, we use the forward Euler scheme Uin+1 − Uin . δt Together with the initial condition and the source Fin = f (xi , tn ), we then end with the system (4)

(5) (6)

ut (xi , tn ) ≈

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 = (U1n , U2n , . . . , UN n

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λ

   .  

NUMERICAL METHODS FOR PARABOLIC EQUATIONS

3

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. Therefore it is also called time marching. 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 the stability in time. When f = 0, i.e., the heat equation without the source, in the continuous level, the solution should be exponential decay. In the discrete level, we have U n+1 = AU n and want to control the magnitude of U in certain norms. 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 /2. Although in one step, it is efficient, it will be very expensive to reach the solution at the ending time T by moving forward with 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 . The error is denoted by E n = U n − unI . Estimate of the truncation error is straightforward using the Taylor expansion.

4

LONG CHEN

Lemma 2.4. Suppose u is smooth enough. For the forward Euler method, we have kτ n k∞ ≤ C1 δt + C2 h2 . 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 n

n

0

E = A E − δt

n−1 X

An−l−1 τ l .

l=1

Since E 0 = 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. Next we introduce the backward Euler method to remove the strong constraint of the time step-size for the stability. The method is simply using the backward difference to approximate the time derivative. We list the resulting linear systems: (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∞ .


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 the 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 the backward Euler method to the example you constructed in Exercise 2.3 to show that 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, (47) 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 the truncation error is improved to second order in time 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 assumptions on λ. What is the weakest assumption on λ you can impose? 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 A popular 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 grid points 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 , with vˆk = (V , φk )h .

k=1

The discrete Parseval identity is kV kh = kˆ v k,


7

where v ˆ = (ˆ v1 , · · · , vˆN )t is the vector formed by the coefficients. Theîn 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

Note that k · kh = h1/2 k · k. The stability in the scaled l2 norm is equivalent to that in k · k norm. Here the scaling h is chosen to be consistent with the scaling of discrete Laplace operator. Now we analyze the l2 -stability of the three numerical methods we have discussed. Recall that λ = δ/h2 . Forward Euler Method. 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 stability (with respect to N ) implies λ ≤ 1/2, i.e., sup ρ(AN ) ≤ 1 =⇒ λ ≤ 1/2. N

Therefore we obtain the same condition as that for the maximum norm stability. Backward Euler Method. U n+1 = A−1 U n ,

A = diag(−λ, 1 + 2λ, −λ).

Thus ρ(A−1 ) =

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.

8

LONG CHEN

Crank-Nicolson Method. 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 |

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? You can google “Crank Nicolson method maxnorm stability” to read more. 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 uniform stability is obtained by considering inequality max ρ(θ).

0≤θ≤π

To fascinate the calculation, one can factor out ρn eiθj and consider only the difference of indices. The eigen-function is changed to eiθ = cos(θ) + i sin(θ) to account for possible dift ferent boundary conditions. The the frequency of discrete eigenvectors φk = eiθk (1:N ) is changed to a continuous variable θ in (12). We skip the index k and take the maximum of θ over [0, π] while the range of the discrete angle θk satisfying hπ ≤ θk ≤ (1 − h)π. As an example, the readers are encouraged to do the following exercise. Use one exercise as an example. Exercise 3.4 (The θ method). For θ ∈ [0, 1], we use the scheme: for 1 ≤ i ≤ N, 1 ≤ n ≤ M n+1 n U n+1 + Ui+1 − 2Uin+1 U n + Ui+1 − 2Uin Uin+1 − Uin (13) = (1 − θ) i−1 + θ i−1 , 2 2 δt h h and for 1 ≤ i ≤ N, n = 0 Ui0 = u0 (xi ). 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 2δt h


9

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.

1

The other U , 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. 4. VARIATIONAL FORMULATION AND ENERGY ESTIMATE We multiply a test function v ∈ H01 (Ω) and apply the integration by part to 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 (Ω).

where a(u, v) = (∇u, ∇v) and (·, ·) denotes the L2 -inner product. 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 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 t ∈ (0, T )a.e. (16) u(·, 0) = u0 where h·, ·i is the duality pair of H −1 and H01 . We assume the equation (16) is well posed. For the existence and uniqueness of the solution [u, ut ], we refer to [1]. In most places, we shall still use the formulation (15) and assume f, ut ∈ L2 . 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 (almost everywhere) in time. In particular, one has to justify the point value u(·, 0) does make sense for a L2 type function which can be proved by the regularity theory of the heat equation. 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 operator L : L2 (0, T ; H01 (Ω)) → L2 (0, T ; H −1 (Ω)) × L2 (Ω)

10

LONG CHEN

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. The spatial boundary condition is build into the space H01 (Ω). In most places, when it is clear from the context, we also use L for the differential operator only. We shall prove several stability results of L which are known as energy estimates in [3]. Theorem 4.2. Suppose [u, ut ] is the solution of (15) 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

ku(t)k2 +

(18)

t

Z

|u(s)|21 ds ≤ ku0 k2 +

0

|u(t)|21

(19)

Z

kf (s)k2−1 ds,

0

t 2

kut (s)k ds ≤

+

t

Z

|u0 |21

Z

t

+

0

kf (s)k2 ds.

0

Proof. The solution is defined via the action of all test functions. The art of the energy estimate is to choose an appropriate 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 d d 1 2 (u )t = kuk2 = kuk kuk; • (ut , u) = 2 dt dt Ω 2 • a(u, u) = |u|21 ;

1 (kf k2−1 + |u|21 ). 2 The inequality (17) is an easy consequence of the following inequality • |(f, u)| ≤ kf kkuk or |(f, u)| ≤ kf k−1 |u|1 ≤

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.


11

From (18), we can obtain a better stability for the operator L : L2 (0, T ; H01 (Ω)) → L2 (0, T ; H −1 (Ω)) × L2 (Ω) as kukL2 (0,T ;H01 (Ω)) ≤ kuk0 + kf kL2 (0,T ;H −1 (Ω)) . Since the equation is posed a.e for t, we could also obtain maximum-norm 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). Exercise 4.4. Prove the energy estimate (20)

kuk2 ≤ e−λt ku0 k2 +

Z

t

e−λ(t−s)) kf k2−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 ⊂ H01 (Ω), 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 (indeed differential) function of t. The initial condition u0 is approximated by u0,h ∈ Vh and choices 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 (22)

M u˙ + Au = f ,

where u = (u1 , · · · , uN )t is the coefficient vector, M , A are the mass matrix and the stiffness matrix, respectively, and f = (f1 , · · · , fN )t . For the linear finite element, one can chose a numerical quadrature to make M diagonal, which is known as the mass lumping, so that the ODE system (22) can be solved efficiently by mature ODE solvers. We shall apply our abstract error analysis developed in Chapter: Unified Error Analysis to estimate the error u − uh in certain norms. 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

12

LONG CHEN

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 kfh k−1,h = sup , for fh ∈ V0h . 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,

in V0h , ∀t ∈ (0, T ] a.e.

uh (·, 0) = u0,h . 5.1. Stability. Adapt the proof of the energy estimate for L, we can obtain similar stability results for Lh . The proof is almost identical and thus skipped here. Theorem 5.1. Suppose uh satisfy Lh uh = fh , uh (·, 0) = u0,h , then Z t kuh (t)k ≤ ku0,h k + (23) kfh (s)k ds 0 Z t Z t kuh (t)k2 + |uh (s)|21 ds ≤ ku0,h k2 + (24) kfh (s)k2−1,h ds, 0 0 Z t Z t |uh (t)|21 + k∂t uh (s)k2 ds ≤ |u0,h |21 + (25) kfh (s)k2 ds. 0

0

Note that in the energy estimate (24), 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 (26) (27)

Lh (Rh u − uh ) = ∂t (Rh u − Qh u), t > 0, in V0h , (Rh u − uh )(·, 0) = Rh u0 − u0,h .

Proof. Let A = −∆. By our definition of consistency, the error equation is: for t > 0 Lh (Rh u − uh ) = Lh Rh u − Qh Lu = ∂t (Rh u − Qh u) + (ARh u − Qh Au). The desired result then follows by noting that ARh = Qh A in V0h . One can also verify (26) by writing out the variational formulation.


13

The error equation (26) holds in V0h which is a weak topology. The motivation to 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 simply 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 the discrete error Rh u − uh . Theorem 5.3. The solution uh of (21) satisfy the following error estimate Z t (28) kRh u − uh k ≤ ku0,h − Rh u0 k + k∂t Qh u − ∂t Rh uk ds 0

(29) t

Z

2

|(uh − Rh u)|21

kRh u − uh k +

2

Z

ds ≤ ku0,h − Rh u0 k +

0

t

k∂t Qh u − ∂t Rh uk2−1,h ds,

0

(30) |Rh u − uh |21 +

Z

t

k∂t (Rh u − uh )k2 ds ≤ |u0,h − Rh u0 |21 +

0

Z

t

k∂t Qh u − ∂t Rh uk2 ds.

0

We then estimate the two terms involved 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 coming from approximation of 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 ,

(31) 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 goes to infinity. So in practice, we can choose the simple nodal interpolation. On the estimate of the second term, we assume ut ∈ H 2 (Ω) and assume H 2 -regularity result hold for Poisson equation (for example, the domain is smooth or convex), then (32)

k∂t Qh u − ∂t Rh uk = kQh (I − Rh )ut k ≤ k(I − Rh )ut k . h2 kut k2 .

The negative norm can be bounded by the L2 -norm as k∂t Qh u − ∂t Rh uk−1,h ≤ k∂t Qh u − ∂t Rh uk−1 ≤ Ck∂t Qh u − ∂t Rh uk . h2 kut k2 . When using quadratic and above polynomial, we can prove a stronger estimate for the negative norm (33)

ku − Rh uk−1 ≤ Chr+1 kukr

where r is the order of the element, then we could obtain a superconvergence in L2 -norm Z t 1/2 r+1 2 kRh u − uh k ≤ Ch kut kr ds . 0

14

LONG CHEN

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 (31). We then have Z t −1 2 (34) 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 (31) satisfy the following optimal order error estimate: Z t −1 2 (35) h |u − uh |1 + ku − uh k ≤ Ch ku0 k2 + kut k2 ds . 0

5.4. *Superconvergence and maximum norm estimate. The error eh = Rh u − uh satisfies the evolution equation (26) with eh (0) = 0 with the chose u0,h = Rh u0 , i.e., ∂t eh + Ah eh = Qh ∂t (Rh u − u).

(36) Therefore t

Z

exp(−Ah (t − s))τh ds,

eh (t) =

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]. Lemma 5.6 (Smoothing property of the heat kernel). For τh ∈ Vh , we have

Z t

(37) max exp(−Ah (t − s))Ah τh ds

≤ 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 (t−s) if (t − s) ≤ λ−1

 λ M e M ,

−Ah (t−s) −1 Ah ≤ (t − s) 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)

max e Ah τh ds

≤ C| log h| max kτh k. 0≤t≤T

0

0≤t≤T

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 (38)

max |Rh u − uh |1 ≤ C| log h| h2 max kut k2 .

0≤t≤T

0≤t≤T


15

When ut ∈ L2 (0, T ; H 2 (Ω)), then 2

Z

|Rh u − uh |1 ≤ Ch

(39)

t

kut k22

1/2 ds .

0

Proof. We multiply A1/2 to (36) and apply Lemma 5.6 to get

Z

T

1/2 −Ah (T −s) 1/2 e Ah τh ds |eh (T )|1 = kAh eh (T )k =

0

Z

T

−1/2 e−Ah (T −s) Ah ds max kAh τh k ≤

0≤t≤T

0 ≤ 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 (39), we use the energy estimate (30).

Since the optimal convergent rate for |u − Rh u|1 or |u − uh |1 is only first order, the estimate (38) and (39)) 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 ∈ 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 (40) , k(u − uh )(t)k∞ ≤ C| log h|h kuk2,∞ + kut k2 ds 0

(41)

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. 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 .

(42)

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 r (Ω)). Let uh be the solution of (21) with u0,h = Rh u0 . Then (43)

max k(Rh u − uh )(t)k ≤ C| log h|hr+1 max k(ut )(t)k2r .

0≤t≤T

0≤t≤T

16

LONG CHEN

When ut ∈ L2 (0, T ; H r (Ω)), then (44)

r+1

Z

kRh u − uh k ≤ Ch

t

kut k2r

1/2 ds .

0

Proof. From (36), we have

Z

T

−Ah (T −s) e τh ds keh (T )k =

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

−Ah (T −s) −1/2 e 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 (44), we simply use the energy estimate (29) and (42).

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 (45) 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. 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 # r 2 (46) . k(u − uh )(t)k∞ ≤ Ch kukr,∞ + kut kr ds 0

6. F INITE ELEMENT METHODS : SEMI - DISCRETIZATION IN TIME In this section, we consider the semi-discretization in time. We first discretize the time interval (0, T ) into 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 Method. u0 = u0 un − un−1 + Aun−1 = f n−1 . δt • Backward Euler Method. u0 = u0 un − un−1 + Aun = f n . δt


17

• Crank-Nicolson Method. 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 as 1 1 n (u − un−1 , v) + (∇un + ∇un−1 , ∇v) = (f n−1/2 , v) for all v ∈ H01 . (47) δt 2 We now study the stability of these schemes. We rewrite the Backward Euler method as (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 +

(48)

n X

δtkf k k.

k=1

The stability (49) 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 +

(49)

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. Which one to chose is problem dependent. The energy method can be adapted to the semi-discretization in time easily. For example, we chose v = (un + un−1 )/2 in (47) 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) (50)

n 2

ku k +

n X

δt|uk−1/2 |21

0 2

≤ ku k +

k=1

n X

δtkf k−1/2 k2−1 .

k=1

Exercise 6.2. Study the stability of the forward and backward 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 ) 2 (51) − ∂ u(·, t ) t n−1/2 ≤ C(δt) . δt By the stability result, we then get kI n u − un k ≤ Ctn δt2 .

18

LONG CHEN

From the consistency error estimate (51), one can easily see the backward and forward Euler methods are only first order in time. High order schemes, e.g., Runge-Kutta method, can be constructed by using a high order approximation of the integral. We first rewrite the differential equation (in time) as an integral equation: Z Z Z 1 tn 1 tn 1 tn ut (·, t) dt + Au(·, t) dt = f (·, t) dt. δt tn−1 δt tn−1 δt tn−1 The first term is simply (u(·, tn ) − u(·, tn−1 )/δt. The right hand side can be approximated as accurately as we want since f is known. The subtle part is the middle one. We can apply high order quadrature using more points but u is only known at tn−1 . So approximation at quadrature points should be computed first. For details, we refer to [?]. Add more on high order (explicit and implicit) discretization in time. 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 Method. u0h = u0,h (52)

1 n (u − un−1 , vh ) + a(un−1 , vh ) = hfhn−1 , vh i, h h δt h

∀vh ∈ Vh , 1 ≤ n ≤ N

Backward Euler Method. u0h = u0,h (53)

1 n (u − un−1 , vh ) + a(unh , vh ) = hfhn , vh i, h δt h

∀vh ∈ Vh , 1 ≤ n ≤ N

Crank-Nicolson Method. u0h = u0,h (54) 1 n n−1/2 (u − un−1 , vh ) + a((unh + un−1 )/2, vh ) = hfh , vh i, ∀vh ∈ Vh , 1 ≤ n ≤ N h h δt h We write the semi-discretization and the full discretization as L → Lh → Lnh . The discrete error is Rh u(tn ) − unh . The setting is R 1/2 T • X = L2 (0, T ; H01 (Ω)), and kukX = 0 |u(t)|21 dt 1/2 R T • Y = L2 (0, T ; H −1 (Ω), and kf kY = 0 kf (t)k2−1 dt 1/2 P N k 2 • Xhδt = Vh × Vh × · · · × Vh , and kuh kXhδt = k=0 δt|uh |1 P 1/2 N k 2 • Yhδt = V0h × V0h × · · · × V0h , and kfh kYhδt = δtkf . k h −1 k=0 • I n Rh : L2 (0, T ; H01 (Ω)) → 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 ).


19

• Lnh : Xhδt → Yhδt × Vh is u0h = u0,h and for n ≥ 1  + n n  Forward Euler, Dt uh − ∆uh n n − n n Lh uh = Dt uh − ∆uh 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 present the stability in L2 norm. Again denoted by Ah = −∆|Vh . Forward Euler Method. unh = (I − δtAh )un−1 + δtfhn−1 . h

(55) Backward Euler Method.

unh = (I + δtAh )−1 (un−1 + δtfhn ). h

(56) Crank-Nicolson Method. (57)

unh

1 1 n−1/2 n−1 −1 (I − δtAh )uh + δtfh . = (I + δtAh ) 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 Method. ρ(I − δtAh ) = |1 − δtλmax (Ah )| ≤ 1, provided δt ≤

2 . λmax (Ah )

Note that λmax (Ah ) = O(h−2 ). We need the time step is in the size of h2 to make the forward Euler stable. Backward Euler Method. For any δt > 0, since Ah is SPD, i.e., λmin (Ah ) > 0 ρ((I + δtAh )−1 ) = (1 + δtλmin (Ah ))−1 ≤ 1. Crank-Nicolson Method. For any δt > 0, |1 − 21 δt λk (Ah )| 1 1 ρ((I + δtAh )−1 (I − δt Ah )) = max ≤ 1, 1≤k≤N 1 + 1 δt λk (Ah ) 2 2 2 We summarize the stability as the following theorem which is a discrete version of (17).

20

LONG CHEN

Theorem 7.1. For forward Euler method, when δt < 1/λmax (Ah ) kunh k

(58)

≤

ku0h k

+

n−1 X

δtkfhk k.

k=0

For backward Euler kunh k ≤ ku0h k +

(59)

n−1 X

δtkfhk+1 k.

k=0

For Crank-Nicolson method kunh k ≤ ku0h k +

(60)

n−1 X

k+1/2

δtkfh

k.

k=0

It is interesting to note that the last term in the stability result (58), (59), and (60) 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

(Rh ut − ut ) dsk ≤ δt−1 h2

tn−1

Z

tn

kut k2 ds. tn−1

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 .


21

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 (54), 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Ía-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ée. 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.