Finite Difference Methods for Solving Differential Equations

(1)

FINITE DIFFERENCE METHODS

FOR

SOLVING DIFFERENTIAL EQUATIONS

I-Liang Chern

Department of Mathematics National Taiwan University

(2)

(3)

Introduction

The goal of this course is to provide numerical analysis background for finite difference methods for solving partial differential equations. The focuses are the stability and convergence theory. The partial differential equations to be discussed include

• parabolic equations, • elliptic equations,

• hyperbolic conservation laws.

1.1 Finite Difference Approximation

Our goal is to appriximate differential operators by finite difference operators. How to perform approximation? What is the error so produced? We shall assume the underlying functionu :_{R → R} is smooth. Let us define the following finite difference operators:

• Forward difference: D+u(x) := u(x+h)−u(x)_h ,

• Backward difference: D−u(x) := u(x)−u(x−h)h ,

• Centered difference: D0u(x) := u(x+h)−u(x−h)_2h .

Here,h is called the mesh size. By Taylor expansion, we can get • u′(x) = D+u(x) + O(h),

• u′(x) = D₋u(x) + O(h), • u′(x) = D0u(x) + O(h2).

We can also approximate u′_{(x) by other finite difference operators with higher order errors. For}

example,

u′(x) = D3u(x) + O(h3),

(8)

where

D3u(x) = 1

6h(2u(x + h) + 3u(x) − 6u(x − h) + u(x − 2h)) .

These formulae can be derived by performing Taylor expansion ofu at x. For instance, we expand

u(x + h) = u(x) + u′(x)h +h 2 2 u ′′_{(x) +} h3 3!u ′′′_{(x) + · · ·} u(x − h) = u(x) − u′(x)h +h 2 2 u ′′_{(x) −} h3 3!u ′′′_{(x) + · · · .}

Substracting these two equations yields

u(x + h) − u(x − h) = 2u′(x)h +2h 3 3! u ′′′_{(x) + · · · .} This gives u′(x) = D0u(x) −h 2 3!u ′′′_{(x) + · · · = D} 0u(x) + O(h2).

In general, we can derive finite difference approximation foru(k)at specific pointx by the values ofu at some nearby stencil points xj, j = 0, ..., n with n ≥ k. That is,

u(k)(x) =

n

X

j=0

cju(xj) + O(hp−k+1)

for some_{p as larger as possible. Here, the mesh size h denotes max{|x}i− xj|}. As we shall see

later that we can choosep = n. To find the coefficients cj, j = 0, ..., n, we make Taylor expansion

ofu(xj) about the point x:

u(xj) = p X i=0 1 i!(xj− x) i_u(i)_{(x) + O(h}p+1_).

We plug this expansion formula of u(xj) into the finite difference approximation formula for

u(k)_(x): u(k)(x) = n X j=0 cj p X i=0 1 i!(xj − x) i_u(i)_{(x) + O(h}p−k+1_).

Comparing both sides, we obtain

n X j=0 1 i! (xj − x)i hi cj = 1 ifi = k 0 otherwise , for i = 0, ..., p

There arep + 1 equations here, it is natural to choose p = n to match the n + 1 unknowns. This is a n×n Vandermonde system. It is nonsingular if xiare different. The matlab code fdcoeffV(k,xbar,x)

(9)

1.2. BASIC NUMERICAL METHODS FOR ORDINARY DIFFERENTIAL EQUATIONS 5

In the case of uniform grid, using central finite differencing, we can get high order approxima-tion. For instance, letxj = jh, where j ∈ Z. Let uj = u(xj). Then

u′(0) = u1− u−1 h + O(h 2₎ u′′(0) = u1− 2u0+ u−1 h2 + O(h 2₎ u(3) = 1 2h3(u2− 2u1+ 2u0− 2u−1+ u−2) + O(h 2_). Homeworks.

1. Consider xi = ih, i = 0, ..., n. Let ¯x = xm. Find the coefficientsci foru(k)(¯x) and the

coefficient of the leading truncation error for the following cases: • k = 1, n = 2, 3, m = 0, 1, 2, 3.

• k = 2, n = 2, m = 0, 1, 2.

1.2 Basic Numerical Methods for Ordinary Differential Equations

The basic assumption to design numerical algorithm for ordinary differential equations is the smooth-ness of the solutions, which is in general valid provided the coefficients are also smooth. Basic designning techniques include numerical interpolation, numerical integration, and finite difference approximation.

Euler method

Euler method is the simplest numerical integrator for ODEs. The ODE

y′ = f (t, y) (2.1)

is discretized by

yn+1= yn+ kf (tn, yn). (2.2) Here,k is time step size of the discretization. This method is called the forward Euler method. It simply replacedy/dt(tn) by the forward finite difference (yn+1_{− y}n)/k. By Taylor expansion, the local truncation error is

τn:= y′(tn_{) −} y(t

n+1_{) − y(t}n₎

k = O(k)

Leten:= yn_{− y(t}n) be the true error.

Theorem 2.1. Assume_{f ∈ C}1 and suppose the solution y′ = f (t, y) with y(0) = y0 exists on

[0, T ]. Then the Euler method converges at any t ∈ [0, T ]. In fact, the true error en has the following estimate:

|en| ≤ e

λt

λ O(k) → 0, as n → ∞. (2.3)

(10)

Proof. From the regularity of the solution, we have_{y ∈ C}2[0, T ] and

y(tn+1) = y(tn) + kf (tn, y(tn)) + kτn. (2.4)

Taking difference of (2.2) and (2.4), we obtain

|en+1| ≤ |en| + k|f (tn, yn_{) − f (t}n, y(tn_{))| + k|τ}n_| ≤ (1 + kλ)|en| + k|τn|.

where

|f (t, x) − f (t, y)| ≤ λ|x − y|.

The finite difference inequality has a fundamental solution Gn = (1 + λk)n, which is positive providedk is small. Multiplying above equation by (1 + λk)−n−1_{, we obtain}

|em+1|G−m−1 ≤ |em|G−m+ kG−m−1_|τm_|. Summing in_{m from m = 0 to n − 1, we get}

|en| ≤ n−1 X m=0 Gn−m−1_k|τm_{| ≤} n−1_X m=0 GmO(k2) = G n_{− 1} G − 1 O(k 2 ) ≤ G n λ O(k) ≤ eλt λ O(k),

wheret = nk and we have used (1 + λk)n_{≤ e}λt. Remarks.

1. The theorem states that the numerical method converges in[0, T ] as long as the solutions of the ODE exists.

2. The error isO(k) if the solution is in C2[0, T ].

3. The proof above relies on the existence and smoothness of the solution. However, one can also use this approach to prove the local existence theorem by showning the approximate solutions generated by the Euler method form a Cauchy sequence.

Backward Euler method

In many applications, the system is relaxed to a stable solution in a very short period of time. For instance, consider

y′ = y − y¯ τ .

The corresponding solution_{y(t) → ¯y as t ∼ O(τ). In the above forward Euler method, practically,} we should require

(11)

1.2. BASIC NUMERICAL METHODS FOR ORDINARY DIFFERENTIAL EQUATIONS 7

in order to haveGn_{remain bounded. Here,}_{λ is the Lipschitz constant. In the present case, λ = 1/τ .}

Ifτ is very small, the the above forward Euler method will require very small k and lead to inefficient computation. In general, forward Euler method is inefficient (require smallk) if

max ∂f (t, y) ∂y >> 1.

In the case∂f /∂y >> 1, we have no choice to resolve details. We have to take a very small k. However, if∂f /∂y < 0, say for example, y′_{= −λy with λ >> 1. then the backward Euler method} is recommended:

yn+1= yn+ kf (tn+1, yn+1). The error satisfies

|en+1| ≤ en− λk|en+1| + O(k2)

The corresponding fundamental solution isGn := (1 + λk)−n. Notice that the error satisfies

|en| ≤ n−1 X m=0 (1 + λk)−mO(k2) ≤ (1 + λk)_λk−n+1O(k2) ≤ e−λT_λ O(k).

There is no restriction on the size ofk.

Leap frog method

We integratey′ _{= f (t, y) from t}n−1_to_tn+1_:

y(tn+1_{) − y(t}n−1) = Z tn+1

tn−1

f (τ, y(τ )) dτ.

We apply the midpoint rule for numerical integration, we then get y(tn+1_{) − y(t}n−1) = 2kf (tn, y(tn)) + O(k3). The midpoint method (or called leapfrog method) is

yn+1_{− y}n−1= 2kf (tn, yn). (2.5) This is a two-step explicit method.

Homeworks.

1. Prove the convergence theorem for the backward Euler method.

Hint: show thaten+1_{≤ e}n+ (1 + kλ)en+1+ kτn, whereλ is the Lipschitz constant of f . 2. Prove the convergence theorem for the leap-frog method.

Hint: consider the systemyn

(12)

1.3 Runge-Kutta methods

The Runge-Kutta method (RK) is a strategy to integrateR_ttnn+1f dτ by some quadrature method.

RK2 For instance, a second order RK, denoted by RK2, is based on the trapezoidal rule of nu-merical integration. First, we integrate the ODEy′_{= f (t, y) to get}

y(tn+1_{) − y}n= Z tn+1

tn

f (τ, y(τ )) dτ.

Next, this integration is approximated by

1/2(f (tn, yn) + f (tn, yn+1))k.

The latter term involvesyn+1. An explicit Runge-Kutta method approximateyn+1byyn+kf (tn, yn). Thus, RK2 reads ξ1 = f (tn, yn) yn+1 = yn+k 2(f (t n_{, y}n_{) + f (t}n+1_{, y}n_{+ kξ} 1)).

Another kind of RK2 is based on the midpoint rule of integration. It reads

ξ1 = f (tn, yn)

yn+1 = yn+ kf (tn+1/2, yn+ k 2ξ1) The truncation error of RK2 is

y(tn+1_{) − y(t}n) = yn+1_{− y(t}n) + O(k3). RK4 The 4th order Runge-Kutta method has the form

yn+1 = yn+k 6(ξ1+ 2ξ2+ 2ξ3+ ξ4) ξ1 = f (tn, yn) ξ2 = f (tn+ 1 2k, y n₊k 2ξ1) ξ3 = f (tn+1 2k, y n₊k 2ξ2) ξ4 = f (tn+ k, yn+ kξ3)

The truncation error of RK4 is

(13)

1.3. RUNGE-KUTTA METHODS 9

General explicit Runge-Kutta methods The method takes the following general form

yn+1= yn+ k s X i=1 biξi, where ξ1 = f (tn, yn), ξ2 = f (tn+ c2k, yn+ ka21ξ1), ξ3 = f (tn+ c3k, yn+ ka31ξ1+ ka32ξ2), .. . ξs = f (tn+ csk, yn+ k(as1ξ1+ · · · + as,s−1ξs−1)).

We need to specifys (the number of stages), the coefficients aij(1 ≤ j < i ≤ s), bi(i = 1, ..., s)

andci(i = 2, ..., s). We list them in the following Butcher table.

There are_{s(s−1)/2+s+(s−1) unknowns to be determined for a specific scheme. We require the} 0 c2 a21 c3 a31 a32 .. . ... . .. cs as1 as2 · · · as,s−1 b1 b2 · · · bs−1 bs

truncation error to beO(kp+1). To find these coefficients, we need to expand the truncation error formula

y(tn+1_{) − y}n= yn+1_{− y}n+ O(kp+1)

about(tn, yn_{) in terms of derivatives of y(·) at t}n. Then we can get linear equations for the coeffi-cients.

Convergence proof, an example Let us see the proof of the convergence of the two stage Runge-Kutta method. The scheme can be expressed as

yn+1= yn+ kΨ(yn, tn, k) (3.6)

where

Ψ(yn, tn, k) := f (y + 1

2kf (y)). (3.7)

Suppose_{y(·) is a true solution, the corresponding truncation error} τn:= y(t

n+1_{) − y(t}n₎

k − Ψ(y(t

(14)

Thus, the true solution satisfies

y(tn+1_{) − y(t}n) = kΨ(y(tn), tn, k) + kτn The true erroren:= yn_{− y(t}n) satisfies

en+1= en+ k(Ψ(yn, tn_{, k) − Ψ(y(t}n), tn_{, k)) − kτ}n. This implies

|en+1| ≤ |en| + kλ′|en| + k|τn|,

whereλ′is the Lipschitz constant ofΨ(y, t, k) with respect to y. Hence, we get

|en| ≤ (1 + kλ′)n_|e0_{| + k} n−1 X m=0 (1 + kλ′)n−1−m_|τm_| ≤ eλ′t|e0| +e λ′_t λ′ maxm |τ m | Reference:

• Lloyd N. Trefethen, Finite Difference and Sprectral Methods for Ordinary and Partial Differ-ential Equations,

• Randy LeVeque,

• You may also google Runge-Kutta method to get more references.

1.4 Multistep methods

The idea of multi-step method is to derive a relation between, for instance, yn+1_{, y}n_{, y}n−1 _and

y′n andy′n−1 so that the corresponding truncation is small. The simplest multistep method is the midpoint method. Supposeynandyn−1is given. The new stateyn+1is defined by

yn+1_{− y}n−1= 2ky′n= 2kf (tn, yn) The truncation error is

τn:= 1 k y(t

n+1

) − y(tn−1) − 2ky′(tn) = O(k2). Thus, the method is second order.

We can also design a method which involvesyn+1, yn, yn−1andy′n, y′n−1. For instance,

yn+1= yn+k 2(3f (y

n

(15)

1.4. MULTISTEP METHODS 11

The truncation error

τn:= 1 k yn+1_{− y}n+k 2(3f (y n_{) − f (y}n−1₎₎ = O(k2).

A generalr-step multistep method can be written in the form

r X m=0 amyn+1−r+m = k r X m=0 bmy′n+1−r+m= k r X m=0 bmfn+1−r+m. (4.8)

We will always assumear6= 0. When br= 0 the method is explicit; otherwise it is implicit. Ror a

smooth solution of (2.1), we define the truncation errorτnto be

τn:= 1 k r X m=0 amy(tn+1−r+m) − k r X m=0 bmy′(tn+1−r+m) !

Definition 4.1. A multispep method is called of orderp if τn _{= O(k}p_{) uniformly in n. It is called}

consistent ifτn_{(k) → 0 uniformly in n as k → 0.}

Remark. Whenf is smooth, the solution of ODE y′ = f (y) is also smooth. Then the truncation is a smooth function of_{k. In this case, τ (k) → 0 is equivalent to τ(k) = O(k) as k → 0. Let us set} am = 0, bm = 0 for m > r for notational convinience. Taking Taylor expansion about tn+1−r, we

get kτn = r X m=0 am ∞ X j=0 1 j!y (j)_(mk)j_{− k} r X m=0 bm ∞ X j=1 1 (j − 1)!y (j)_(mk)j−1 = r X m=0 am ! y(0)+ ∞ X j=1 1 j! r X m=0 mjam− jmj−1bmkjy(j) = r X m=0 am ! y(0)+ ∞ X j=1 1 j! r X m=0 mj−1(mam− jbm) kjy(j) = ∞ X j=0 1 j! X m=0 mj−1(mam− jbm) kjy(j) = ∞ X j=0 Cjy(j).

(16)

b: C0 = a0+ · · · + ar C1 = (a1+ 2a2+ · · · rar) − (b0+ · · · + br) C2 = 1 2 (a1+ 2 2_a 2+ · · · r2ar) − 2(b1+ · · · + rbr) .. . Cp = r X m=0 mp p! am− r X m=1 mp−1 (p − 1)!bm To obtain orderp scheme, we need to require

Cj = 0, for j = 0, ..., p.

There are_{2(r + 1) unknowns for the coefficients {a}m}rm=0,{bm}rm=0. In principle, we can choose

p = 2r + 1 to have the same number of equations. Unfortunately, there is some limitation from stability criterion which we shall be explained in the next section. The order of accuracyp satisfies

p ≤    r + 2 if r is even, r + 1 if r is odd, r if it is an explicit scheme.

This is the first Dahlquist stability barrier. We shall not discuss here. See Trefethen. Let us see some concrete examples below.

Explicit Adams-Bashforth schemes Whenbr= 0, the method is explicit. Here are some

exam-ples of the explicit schemes called Adams-Bashforth schemes, wherear = 1:

• 1-step: yn+1_{= y}n_{+ kf (y}n₎

• 2-step: yn+1_{= y}n₊k

2(3f (yn) − f (yn−1))

• 3 step: yn+1 = yn+₁₂k(23f (yn_{) − 16f (y}n−1) + 5f (yn−2)) The step size isr and the order is p = s.

Implicit Adams-Moulton schemes Another examples are the Adams-Moulton schemes, where br6= 0 and and the step size r = p

• 1-step: yn+1_{= y}n₊k

2(f (yn+1) + f (yn))

• 2-step: yn+1_{= y}n₊ k

12(5f (yn+1) + 8f (yn) − f (yn−1))

(17)

1.4. MULTISTEP METHODS 13

Sometimes, we can use an explicit scheme to guessyn+1 _{as a predictor in an implicit scheme.}

Such a method is called a predictor-corrector method. A standard one is the following Adams-Bashforth-Moulton scheme: Its predictor part is the Adams-Bashforth scheme:

ˆ

yn+1= yn+ k

12(23f (y

n

) − 16f (yn−1) + 5f (yn−2)) The corrector is the Adams-Moulton scheme:

yn+1= yn+ k 24(9f (ˆy

n+1_{) + 19f (y}n_{) − 5f (y}n−1_{) + f (y}n−2₎₎

The predictor-corrector is still an explicit scheme. However, for stiff problem, we should use im-plicit scheme instead.

Formal algebra Let us introduce the shift operatorZyn= yn+1, or in continuous sense,Zy(t) = y(t + k). Let D be the differential operator. The Taylor expansion

y(t + k) = y(t) + ky′(t) + 1 2!k

2_D2

y(t) + · · · can be expressed formally as

Zy = 1 + (kD) + 1 2!(kD) 2_{+ · · ·} y = ekDy.

The multistep method can be expressed as

Ly := (a(Z) − kb(Z)D) y =a(ekD_{) − kDb(e}kD)y = (C0+ C1(kD) + · · ·) y.

Here, a(Z) = r X m=0 amZm, b(Z) = r X m=0 bmZm

are the generating functions of_{am} and {bm}. A multistep method is of order p means that

a(ekD_{) − kDb(kD)}y = O((kD)p+1)y. We may abbreviatekD by a symbol κ. The above formula is equivalent to

a(eκ₎

b(eκ₎ = κ + O(κ

p+1_{) as κ → 0.} _(4.9)

We have the following theorem

Theorem 4.2. A multistep method with_{b(1) 6= 0 is of order p if and only if} a(z)

b(z) = log z + O((z − 1)

p+1_{) as z → 1.}

It is consistent if and only if

(18)

Proof. The first formula can be obtain from (4.9) by writingeκ _{= z. For the second formula, we}

expandlog z about 1. We can get

a(z) = b(z) (z − 1) − (z − 1) 2 2 + (z − 1)3 3 + · · · + O((z − 1)p+1). We also expanda(z) and b(z) about z = 1, we can get

a(1) + (z − 1)a′(1) = b(1)(z − 1) + O((z − 1)2). Thus, the scheme is consistent if and only ifa(1) = 0 and a′(1) = b(1).

Homeworks.

1. Consider the linear ODEy′ = λy, derive the finite difference equation using multistep method involvingyn+1_{, y}n_{, y}n−1_and_y′n_and_y′n−1_{for this linear ODE.}

2. Solve the linear finite difference equations derived from previous exercise.

1.5 Linear difference equation

Second-order linear difference equation. In the linear case y′ = λy, the above difference scheme results in a linear difference equation. Let us consider general second order linear difference equation with constant coefficients:

ayn+1+ byn+ cyn−1= 0, (5.10) where_{a 6= 0. To find its general solutions, we try the ansatz y}n= ρnfor some numberρ. Here, the n in ynis an index, whereas then in ρnis a power. Plug this ansatz into the equation, we get

aρn+1+ bρn+ cρn−1 = 0. This leads to

aρ2+ bρ + c = 0.

There are two solutionsρ1andρ2. In caseρ1 6= ρ2, these two solutions are independent. Since the

equation is linear, any linear combination of these two solutions is again a solution. Moreover, the general solution can only depend on two free parameters, namely, oncey0andy−1are known, then {yn_}

n∈Zis uniquely determined. Thus, the general solution is

yn= C1ρn1 + C2ρn2,

whereC1, C2are constants. In case ofρ1 = ρ2, then we can use the two solutions ρn2 andρn1 with

ρ2 6= ρ1, but very closed, to produce another nontrivial solution:

lim

ρ2→ρ1

ρn 2 − ρn1

ρ2− ρ1

This yields the second solution isnρn−1₁ . Thus, the general solution is C1ρn1 + C2nρn−11 .

(19)

1.5. LINEAR DIFFERENCE EQUATION 15

Linear finite difference equation of orderr . We consider general linear finite difference equa-tion of orderr:

aryn+r+ · · · + a0yn= 0, (5.11)

wherear6= 0. Since yn+r can be solved in terms ofyn+r−1, ..., ynfor alln, this equation together

with initial datay0, ..., y−r+1has a unique solution. The solution space isr dimensions.

To find fundamental solutions, we try the ansatz

yn= ρn

for some numberρ. Plug this ansatz into equation, we get

arρn+r+ · · · + a0ρn= 0,

for alln. This implies

a(ρ) := arρr+ · · · + a0 = 0. (5.12)

The polynomiala(ρ) is called the characteristic polynomial of (??) and its roots ρ1, ..., ρrare called

the characteristic roots.

• Simple roots (i.e. ρi 6= ρj, for alli 6= j): The fundamental solutions are ρni,i = 1, ..., r.

• Multiple roots: if ρi is a multiple root with multiplicitymi, then the corresponding

indepen-dent solutions

ρn_i, nρn−1_i , C₂nρn−2_i ..., C_mn_i₋₁ρn−mi+1

i

Here,C_kn_{:= n!/(k!(n − k)!). The solution C}₂nρ_in−2can be derived from differentiation _dρdC₁nρn−1 atρi.

In the case of simple roots, we can express general solution as

yn= C1ρn1 + · · · + Crρnr,

where the constantsC1, ..., Crare determined by

yi= C1ρi1+ · · · + Crρir, i = 0, ..., −r + 1.

System of linear difference equation. The aboverth order linear difference equation is equiva-lent to a first order linear difference system:

A0yn+1= Ayn (5.13) where yn=    y₁n .. . y_rn    =    yn−r+1 .. . yn   

(20)

A0 = I_{(r−1)×(r−1)} 0 0 ar , A =        0 1 0 _{· · ·} 0 0 0 1 _{· · ·} 0 .. . ... ... . .. ... 0 0 0 _{· · ·} 1

−a0 −a1 −a2 · · · −ar−1

       .

We may divide (5.13) by A0and get

yn+1 = Gyn.

We call G the fundamental matrix of (5.13). For this homogeneous equation, the solution is

yn= Gny0

Next, we compute Gnin terms of eigenvalues of G.

In the case that all eigenvaluesρi,i = 1, ..., r of G are distinct, then G can be expressed as

G= TDT−1, D = diag (ρ1, · · · , ρr),

and the column vectors of T are the corresponding eigenvectors.

When the eigenvalues of G have multiple roots, we can normalize it into Jordan blocks:

G= TJT−1, J = diag (J1, · · · , Js),

where the Jordan block Jicorresponds to eigenvalueρiwith multiplicitymi:

Ji =        ρi 1 0 · · · 0 0 ρi 1 · · · 0 .. . ... ... . .. ... 0 0 _{0 · · ·} 1 0 0 _{0 · · · ρ}i        mi×mi .

andPs_i=1mi = r. Indeed, this form also covers the case of distinct eigenvalues.

In the stability analysis below, we are concerned with whether Gnis bounnded. It is easy to see that Gn= TJnT−1, Jn= diag (Jn₁_{, · · · , J}n_s) Jn_i =        ρn_i nρn−1_i C₂nρn−2 _{· · · C}_mn_i₋₁ρn−mi+1 i 0 ρn i nρn−1i · · · Cmni−2ρ n−mi+2 i .. . ... ... . .. ... 0 0 0 _{· · ·} nρn−1_i 0 0 0 _{· · ·} ρn_i        mi×mi . whereC_kn:= _k!(n−k)!n! .

Definition 5.2. The fundamental matrix G is called stable if Gnremains bounded under certain norm_{k · k for all n.}

(21)

1.6. STABILITY ANALYSIS 17

Theorem 5.3. The fundamental matrix G is stable if and only if its eigenvalues satisfy the following

condition:

either_{|ρ| = 1 and ρ is a simple root,}

or_{|ρ| < 1} (5.14)

Proof. It is easy to see that thenth power of a Jordan form Jn

i is bounded if its eigenvalue|ρi| < 1

or ifρi| = 1 but simple. On the other hand, if |ρi| > 1 then Jinis unbounded; or ifρi| = 1 but not

simple, then the termnρn−1_i inJ_inwill be unbounded.

Nonhomogeneous linear finite difference system In general, we consider the nonhomogeneous linear difference system:

yn+1= Gyn+ fn (5.15)

with initial data y0. Its solution can be expressed as yn = Gyn−1+ fn−1 = G(Gyn−2+ fn−2) + fn−1 .. . = Gny0+ n−1 X m=0 Gn−1−mfm Homeworks.

1. Consider the linear ODE

y′ = λy

whereλ considered here can be complex. Study the linear difference equation derived for this ODE by forward Euler method, backward Euler, midpoint. Find its general solutions. 2. Consider linear finite difference equation with source term

ayn+1+ byn+ cyn−1= fn

Given initial datay¯0 andy¯1, find its solution.

3. Find the characteristic roots for the Adams-Bashforth and Adams-Moulton schemes with steps 1-3 for the linear equationy′ = λy.

1.6 Stability analysis

There are two kinds of stability concepts.

• Fix t = nk, the computed solution ynremains bounded as_{n → ∞ (or equivalently, k → 0).} • Fix k > 0, the computed solution ynremains bounded as_{n → ∞.}

(22)

1.6.1 Zero Stability

Our goal is to develop general convergence theory for multistep finite difference method for ODE: y′ = f (t, y) with initial condition y(0) = y0. Anr-step multistep finite difference scheme can be

expressed as Lyn= r X m=0 amyn+1−r+m− k r X m=0 bmf (tn+1−r+m, yn+1−r+m) = 0. (6.16)

Definition 6.3. The truncation errorτn_{(k) for the above finite difference scheme is defined by}

τn(k) := 1 k r X m=0 amy(tn+1−r+m) − k r X m=0 bmf (tn+1−r+m, y(tn+1−r+m)) ! ,

where_{y(·) is a true solution of the ODE.}

Definition 6.4. A difference scheme is called consistent if the corresponding truncation errorτn_{(k) →} 0 uniformly in n as the mesh size k → 0. The scheme is of order p if τn_{(k) = O(k}p_{) uniform in n.}

In the multistep method, the consistency is equivalent to_{τ (k) = O(k) because we assume y(·)} is smooth and the truncation error is a smooth function in_{k. The consistency is τ (k) → 0 as k → 0.} Thus the smoothness ofτ implies τ (k) = O(k).

Definition 6.5. A difference scheme is called zero stable if its solutions at time step n remain

bounded as the mesh size_{k → 0 (nk = T is fixed, according n → ∞).}

The main theorem is the follows. We will postpone its proof at the end of this section.

Theorem 6.4 (Dahlquist). For finite difference schemes for ODEy′ = f (t, y),

consistency _{+ zero-stability ⇔ convergence}

Stability criterion Let us investigate the condition on the coefficientsa and b of an explicit mul-tistep method for the stability

Lyn= 0

to be bounded. We assumear= 1 and br= 0. Let us write it in matrix form:

yn+1= Ayn+ kBfn where A=        0 1 0 1 . .. . .. 0 1

−a0 · · · −ar−2 −ar−1

       , yn=   y n−r · · · yn   ,

(23)

1.6. STABILITY ANALYSIS 19 B=        0 0 0 0 . .. ... 0 0 b0 · · · br−2 br−1        , fn=   f n−r · · · fn   ,

In order to have solution to be bounded for a multistep scheme_{Ly = 0 for arbitrary f , it has at least} to be valid when_{f ≡ 0. In this case, we need to invetigate the boundedness for the homogeneous} equation:

yn+1= Ayn We have seen in the last section that

Theorem 6.5. The necessary and sufficient condition for_kAn_{k to be bounded is that the}

charac-teristic rootsρiof the characteristic equationa(z) = 0 satisfies:

either|ρi| < 1

or_|ρi| = 1 but simple.

Convergence_{⇒ Stability}

Proof. We only need to find anf such that the corresponding multistep is not stable implies that it does not converge. We choose_{f ≡ 0. Since A}nis unbounded, which means there is an eigenvalue ρi with eigenvector yi such that|ρi| > 1 or |ρi| = 1 but not simple. We discuss the formal case.

The latter case can also be prove easily. In the former case, we choosey0 such that y0_{· y}i 6= 0.

Then the correspondingbf yn_{:= A}n_y0 _{will be unbounded. Hence it cannot converge to a constant,}

ask → 0. Here, we use that fact that y0 _{· y}

i 6= 0. We generate bf y0 = (y0r−1, · · · , y0)T by

some explicit scheme starting from y0. The point is that bf y0 depends on the mesh size k and y0_{(k) → (y}0, · · · , y0)T ask → 0. With this, given any yi, we can always construct y0(k) such

that y0_{(k) · y}i 6= 0 when k 6= 0.

Convergence_{⇒ Consistency}

Proof. We need to show thata(1) = 0 and a′(1) = b(1). To show the first, we consider y′= 0 with y(0) = 1. For the second, we consider y′ _{= 1 and y(0) = 0.}

• Show a(1) = 0: We choose y0 = (1, · · · , 1)T. From y1 = Ay0, we get yr_{= −a}0y0− · · · − ar−1yr−1 = −a0− · · · − ar−1.

Sinceyr _{is independent of}_{k, and we should have y}r _{→ 1 as k → 0 (by convergence), we}

(24)

• Show a′_{(1) = b(1). We choose f ≡ 1, y(0) = 0. The corresponding ODE solution is}

y(t) = t. The multistep method gives

a(Z)yn_{− kb(Z)1 = 0.} (6.17)

We write

a(Z) = a′_{(1)(Z − 1) + O((Z − 1)}2), b(Z)1 = b(1). Then the principal part of the above finite difference is

(Z − 1)y − k_ab(1)_′₍₁₎ = 0.

This is an arithmetic series. Its solution isyn= nk_ab(1)′₍₁₎. Indeed, this sequence also satisfies

(6.17) provided its initial data yn has the same form for _{0 ≤ n < r. Since nk = t, the} convergenceyn_{→ t as n → ∞ forces} b(1)

a′₍₁₎ = 1.

Stability + Consistency_{⇒ Convergence}

Proof. We recall that we can express the scheme as

yn+1 = Ayn+ kBfn.

LetY be an exact solution, then plug it into the above =scheme, we get Yn+1= AYn+ kBFn+ kτn. We substract these two and call en:= Yn_{− y}n_{. We get}

en+1= Aen+ kB (Fn_{− f}n) + kτn. The term Fn_{− f}ncan be repressed as

Fn_{− f}n = (f (Yn−r_{) − f (y}n−r_{), · · · , f (Y}n_{) − f (y}n))T = (L_−ren−r_{, · · · , L}0en)T := Lnen where L_−m := Z 1 0 f′(yn−m+ ten−m) dt. Thus, we get en+1 = (A + kBLn) en+ kτn = Gn(k)en+ kτn

withC independent of n and k. The reason is the follows. First, we assume that f is Lipschitz. Thus, the functionsL_−m above are uniformly bounded (independent of_{n). Hence the term kBLk} is uniformly bounded. Second we have a lemma

(25)

1.6. STABILITY ANALYSIS 21

Lemma 6.1. If_kAn_{k is bounded and kB}nk are uniformly bounded, then the product

k(A +_n1B1) · · · (A +

1 nBn)k

is also uniformly bounded.

We have en _{≤ G}_n−1en−1+ kτn−1 ≤ Gn−1Gn−2en−2+ k Gn−2τn−2+ τn−1 ≤ Gn−1Gn−2· · · G0e0 +k G_n−2_{· · · G}0τ0+ · · · + Gn−2τn−2+ τn−1

From the lemma, we get

(26)

(27)

Chapter 2

Finite Difference Methods for Linear

Parabolic Equations

2.1 Finite Difference Methods for the Heat Equation

2.1.1 Some discretization methods

Let us start from the simplest parabolic equation, the heat equation:

ut= uxx

Leth = ∆x, k = ∆t be the spatial and temporal mesh sizes. Define xj = jh, j ∈ Z and tn= nk,

n ≥ 0. Let us abbreviate u(xj, tn) by unj. We shall approximateujnbyUjn, whereUjnsatisfies some

finite difference equations.

Spatial discretization : The simplest one is that we use centered finite difference approximation foruxx:

uxx =

uj+1− 2uj + uj−1

h2 + O(h

2₎

This results in the following systems of ODEs

˙ Uj(t) = Uj+1(t) − 2Uj(t) + Uj−1(t) h2 or in vector form ˙ U = 1 h2AU whereU = (U0, U1, ...)t,A = diag (1, −2, 1). 23

(28)

Homeworks.

1. Derive the 4th order centered finite difference approximation foruxx:

uxx =

1

h2(−uj−2+ 16uj−1− 30uj+ 16uj+1− uj+2) + O(h 4_).

2. Derive a 2nd order centered finite difference approximation for(κ(x)ux)x.

Temporal discretization We can apply numerical ODE solvers • Forward Euler method:

Un+1= Un+ k h2AU

n _(1.1)

• Backward Euler method:

Un+1 = Un+ k h2AU n+1 _(1.2) • 2nd order Runge-Kutta (RK2): Un+1_{− U}n= k h2AU n+1/2_{, U}n+1/2_{= U}n₊ k 2h2AU n _(1.3) • Crank-Nicolson: Un+1_{− U}n= k 2h2(AU n+1_{+ AU}n_). _(1.4)

These linear finite difference equations can be solved formally as

Un+1= GUn where • Forward Euler: G = 1 +_hk2A, • Backward Euler: G = (1 − hk2A)−1, • RK2: G = 1 + _hk2A +1₂ _hk2 2 A2 • Crank-Nicolson: G = 1+ k 2h2A 1−_2h2k A

For the Forward Euler, We may abbreviate it as

U_jn+1= G(U_j−1n , U_jn, U_j+1n ), (1.5)

where

G(U_j−1, Uj, Uj+1) = Uj+ k

(29)

2.1. FINITE DIFFERENCE METHODS FOR THE HEAT EQUATION 25

2.1.2 Stability and Convergence for the Forward Euler method

Our goal is to show under what condition canU_jnconverges tou(xj, tn) as the mesh sizes h, k → 0.

To see this, we first see the local error a true solution can produce. Plug a true solutionu(x, t) into (1.1). We get un+1_j _{− u}n_j = k h2 u n j+1− 2unj + unj−1 + kτ_jn (1.6) where τ_jn= Dt,+unj − (ut)nj − (D+D−unj − (uxx)nj) = O(k) + O(h2).

Leten_j denote forun_j _{− U}_jn. Then substract (1.1) from (1.6), we get en+1_j _{− e}n_j = k h2 e n j+1− 2enj + enj−1 + kτ_jn. (1.7)

This can be expressed in operator form:

en+1 = Gen+ kτn. (1.8)

kenk ≤ kGen−1k + kkτn−1k

≤ kG2en−2_{k + k(kGτ}n−2_{k + kτ}n−1_k)

≤ kGne0_{k + k(kG}n−1τ0_{k + · · · + kGτ}n−2_{k + kτ}n−1_k) Suppose G satisfies the stability condition

kGnU k ≤ CkUk for someC independent of n. Then

kenk ≤ Cke0k + C max_m |τm|. If the local truncation error has the estimate

max

m kτ

m

k = O(h2) + O(k) and the initial errore0satisfies

ke0k = O(h2), then so does the global true error satisfies

kenk = O(h2) + O(k) for all n. The above analysis leads to the following definitions.

Definition 1.6. A finite difference method is called consistent if its local truncation errorτ satisfies

(30)

Definition 1.7. A finite difference schemeUn+1 _{= G}

h,k(Un) is called stable under the norm k · k

in a region_{(h, k) ∈ R if}

kGnh,kU k ≤ CkUk

for alln with nk fixed. Here, C is a constant independent of n.

Definition 1.8. A finite difference method is called convergence if the true error

keh,kk → 0 as h, k → 0.

In the above analysis, we have seen that for forward Euler method for the heat equation,

stability _{+ consistency ⇒ convergence.}

2.2 L

2

Stability – von Neumann Analysis

Since we only deal with smooth solutions in this section, theL2_{-norm or the Sobolev norm is a}

proper norm to our stability analysis. For constant coefficient and scalar case, the von Neumann analysis (via Fourier method) provides a necessary and sufficient condition for stability. For system with constant coefficients, the von Neumann analysis gives a necessary condition for statbility. For systems with variable coefficients, the Kreiss’ matrix theorem provides characterizations of stability condition.

Below, we giveL2_{stability analysis. We use two methods, one is the energy method, the other}

is the Fourier method, that is the von Neumann analysis. We describe the von Neumann analysis below.

Given_{Uj}j∈Z, we define

kUk2 =X

j

|Uj|2

and its Fourier transform

ˆ

U (ξ) = 1 2π

X

Uje−ijξ.

The advantages of Fourier method for analyzing finite difference scheme are

• the shift operator is transformed to a multiplier: d

T U (ξ) = eiξU (ξ),ˆ

where(T U )j := Uj+1;

• the Parseval equility

kUk2 = k ˆU k2 ≡

Z π −π| ˆU (ξ)|

(31)

2.2. L2STABILITY – VON NEUMANN ANALYSIS 27

If a finite difference scheme is expressed as

U_jn+1= (GUn)j = m X i=−l ai(TiUn)j, then [ Un+1_{= b}_{G(ξ) c}_Un_(ξ).

From the Parseval equality,

kUn+1k2 = k [Un+1_k2 = Z π −π| bG(ξ)| 2_{| c}_Un_(ξ)|2_dξ ≤ max ξ | bG(ξ)| 2Z π −π| c Un_(ξ)|2_dξ = | bG|2∞kUk2

Thus a sufficient condition for stability is

| bG|∞≤ 1. (2.9)

Conversely, suppose_{| b}G(ξ0)| > 1, from bG being a smooth function in ξ, we can find ǫ and δ such

that

| bG(ξ)| ≥ 1 + ǫ for all |ξ − ξ0| < δ.

Let us choose an initial dataU0inℓ2such that cU0(ξ) = 1 for |ξ − ξ0| ≤ δ. Then

k cUn_k2 ₌ Z | bG|2n(ξ)|cU0_|2 ≥ Z |ξ−ξ0|≤δ | bG|2n(ξ)|cU0_|2 ≥ (1 + ǫ)2nδ → ∞ as n → ∞

Thus, the scheme can not be stable. We conclude the above discussion by the following theorem. Theorem 2.6. A finite difference scheme

U_jn+1=

m

X

k=−l

akUj+kn

with constant coefficients is stable if and only if

b G(ξ) := m X k=−l ake−ikξ satisfies max −π≤ξ≤π| bG(ξ)| ≤ 1. (2.10)

(32)

Homeworks.

1. Compute the bG for the schemes: Forward Euler, Backward Euler, RK2 and Crank-Nicolson.

2.3 Energy method

We write the finite difference scheme as

U_jn+1= αU_j−1n + βU_jn+ γU_j+1n , (3.11)

where

α, β, γ ≥ 0 and α + β + γ = 1.

We multiply (3.11) byU_jn+1on both sides, apply Cauchy-Schwarz inequality, we get

(U_jn+1)2 = αU_j−1n U_jn+1+ βU_jnU_jn+1+ γU_j+1n U_jn+1

≤ α₂((U_j−1n )2+ (U_jn+1)2) +β 2((U n j)2+ (Ujn+1)2) + γ 2((U n j+1)2+ (Ujn+1)2)

Here, we have used_{α, β, γ ≥ 0. We multiply this inequality by h and sum it over j ∈ Z. Denote}

kUk2:=  X j |Uj|2h   1/2 . We get kUn+1k2 ≤ α 2(kU n_k2_{+ kU}n+1_k2_{) +} β 2(kU n_k2_{+ kU}n+1_k2_{) +}γ 2(kU n_k2_{+ kU}n+1_k2₎ = 1 2(kU n_k2_{+ kU}n+1_k2_).

Here,α + β + γ = 1 is applied. Thus, we get the energy estimate

kUn+1k2 ≤ kUnk2. (3.12)

Homeworks.

(33)

2.4. STABILITY ANALYSIS FOR MONTONE OPERATORS– ENTROPY ESTIMATES 29

2.4 Stability Analysis for Montone Operators– Entropy Estimates

Stbility in the maximum norm

We notice that the action ofG is a convex combinition of U_j−1, Uj, Uj+1, provided

0 < k h2 ≤ 1 2. (4.13) Thus, we get min{U_j−1n , U_jn, U_j+1n _{} ≤ U}_jn+1_{≤ max {U}_j−1n , U_jn, U_j+1n _}. This leads to minjUjn+1≥ minjUjn, maxjUjn+1 ≤ maxjUjn and maxj|Ujn+1| ≤ maxj|Ujn|

Such an operator G is called a monotone operator.

Definition 4.9. An operator G : ℓ∞ _{→ ℓ}∞_{satisfying U} _{≤ V ⇒ GUgV is called a monotone}

operator.

Entropy estimates

The property that Un+1 is a convex combination (average) of Un is very important. Given any convex functionη(u), called entropy function, by Jenson’s inequality, we get

η(U_jn+1_{) ≤ αη(U}_j−1n ) + βη(U_jn) + γη(U_j+1n ) (4.14) Summing over allj and using α + β + γ = 1, we get

X

j

η(U_jn+1_{) ≤}X

j

η(U_jn). (4.15)

This means that the “entropy” decreases in time. In particular, we choose • η(u) = |u|2, we recover theL2stability,

• η(u) = |u|p,_{1 ≤ p < ∞, we get} X j |Ujn+1|p ≤ X j |Ujn|p This leads to _ X j |Ujn+1|ph   1/p ≤  X j |Ujn|ph   1/p ,

(34)

• η(u) = |u − c| for any constant c, we obtain X j |Ujn+1− c| ≤ X j |Ujn− c|

This is called Kruzkov’s entropy estimate. We will see this in hyperbolic quations again.

Homeworks.

1. Show that the solution of the difference equation derived from the RK2 satisfies the entropy estimate. What is the condition required onh and k for such entropy estimate?

2.5 Entropy estimate for backward Euler method

In the backward Euler method, the amplification matrix is given by

G = (I − λA)−1 (5.16)

where

λ = k

h2, A = diag(1, −2, 1).

The matrix_{M := I − λA has the following property:} mii> 0, mij ≤ 0,

X

j6=i

|mij| ≤ mii (5.17)

Such a matrix is called an M-matrix.

Theorem 5.7. The inverse of an M-matrix is a nonnegative matrix, i.e. all its entries are

non-negative.

I shall not prove this general theorem. Instead, I will find the inverse ofM for the above specific M-matrix. Let us express

M = I − λ diag(1, −2, 1) = 1 + 2λ₂ diag_{(−a, 2, −a).} Here,

a = 2λ

1 + 2λ, and 0 < a < 1 if h, k > 0. The general solution of the difference equation

−auj−1+ 2uj − auj+1= 0 (5.18)

has the form:

(35)

2.5. ENTROPY ESTIMATE FOR BACKWARD EULER METHOD 31

whereρ1, ρ2are the characteristic roots, i.e. the roots of the polynomial equation

−aρ2+ 2ρ − a = 0. Thus, ρi= 1 ± √ 1 − a2 a .

From the assumption of the M-matrix,0 < a < 1, we have ρ1< 1 and ρ2 > 1.

Now, we define a fundamental solution:

gj =

(

ρj₁ for_{j ≥ 0} ρj₂ forj < 0 .

We can check thatgj → 0 as |j| → ∞. Moreover, gj satisfies the difference equation (5.18) for

|j| ≥ 1. For j = 0, we have

−ag−1+ 2g0− ag1 = −aρ2−1+ 2 − aρ1 = 2 − a(ρ1+ ρ−12 ) = d

We resetgj ← gj/d. Then we have

X

j

g_i−jmj,k= δi,k.

Thus,M−1 = (g_i−j) is a positive matrix (i.e. all its entries are positive). Further more, X

j

g_i−j = 1 for all i.

Such a matrix appears commonly in probability theory. It is called transition matrix of a Markov chain.

Let us go back to our backward Euler method for the heat equation. We get that Un+1_{= (1 − λA)}−1= GUn, where (GU )i = X j g_i−jUj.

Fromg_i−j > 0 andP_jg_i−j = 1. This is because that gj is the solutionUj1of the finite difference

equation

U_jn+1= U_jn_{+ λdiag(1, −2, 1)U}_jn+1

with initial dataU_j0 = δj. This finite difference equation is conservative, namely,PjUjnis

inde-pendent. This can easily be checked. Thus, we getP_jgj =Pjδj = 1. With this and the positivity

ofgj, we can thinkU_jn+1is a convex combination ofUjnwith weights gj. Thus,G is a monotone

operator.With this property, we can apply Jansen’s inequality to get the entropy estimates:

Theorem 5.8. Letη(u) be a convex function. Let U_jnbe a solution of the difference equation derived from the backward Euler method for the heat equation. Then we have

X

j

η(U_jn_{) ≤}X

j

(36)

Remark. It is important to note that there is no restriction on the mesh sizesh and k for stability for the Backward Euler method.

Homeworks.

1. Can the Crank-Nicolson method for the heat equation satisfy the entropy estimate? What is the condition onh and k?

2.6 Existence Theory

We can prove existence theorem of PDEs through finite difference approximation. In order to do so, let us define continuous and discrete Sobolev spaces and make a connection between them.

The continuous Sobolev space

Hm _{:= {u : R → R|u, u}′, ..., u(m) _{∈ L}2(R)} The discrete Sobolev space for functions defined on gridGh := {jh|j ∈ Z}.

H_hm_{:= {U : G}h → R|U, Dx,+U, ..., Dx,+m U ∈ ℓ2}.

Here,(Dx,+U )nj := (Uj+1n − Ujn)/h

For any discrete functionUj ∈ H_hmwe can construct a functionu in Hmdefined by

u(x) :=X

j

Ujφh(x − xj) (6.20)

whereφh(x) = sinc(x/h). We have

uh(xj) = Uj

It can be shown that

kDxmuhk ≡ kDmx,+U k. (6.21)

Similarly, the spaceL∞_k (H_hm) can be embeded into L∞(Hm) by defining

uh,k(x, t) = X n≥0 X j U_jnφk(t)φh(x)

The discrete norm and the continuous norm are equivalent.

2.6.1 Existence via forward Euler method

The forward Euler method for the heat equationut= uxxreads

U_jn+1= U_jn+ k h2 U n j−1− 2Ujn+ Uj+1n .

(37)

2.6. EXISTENCE THEORY 33

Here, We have seen that we can get the energy estimate:

kUnk ≤ kU0k.

We perform finite difference operation on the above equation, say the forward Euler equation, for instance, letVn

j = (Dx,+U )nj := (Uj+1n − Ujn)/h. Then Vjn satisfies the same finite difference

equation V_jn+1 = V_jn+ k h2 V n j−1− 2Vjn+ Vj+1n .

Thus, it also possesses the same energy estimate. Similar estimate forD2

x,+U . In general, we have

kDx,+m Unk ≤ kDmx,+U0k. (6.22)

If we assume the initial data_{f ∈ H}2, then we getUn_{∈ H}2

h for alln ≥ 0.

Theorem 6.9. If the initial datau0 ∈ Hm,m ≥ 2 and k/h2 ≤ 1/2, then the solution of forward

Euler equation has the estimate

kDx,+m Unk ≤ kDmx,+U0k, kDt,+Unk ≤ kDx,+2 U0k (6.23)

Further, the corresponding smoothing functionuh,k has the same estimate and has a subsequence

converges to a solutionu(x, t) of the original equation.

Proof. The functionsuh,kare unformly bounded in W1,∞(H2). Hence they have a subsequence

converges to a function _{u ∈ W}1,∞(H2) weakly in W1,∞(H2) and strongly in L∞(H1). The functionsuh,ksatisfy uh,k(xj, tn+1) − uh,k(xj, tn) = k h2(uh,k(xj−1, t n_{) − 2u} h,k(xj, tn) + uh,k(xj+1, tn))

Multiply a test smooth function φ, sum over j and n, take summation by part, we can get the subsubsequence converges to a solution ofut= uxxweakly.

2.6.2 A Sharper Energy Estimate for backward Euler method

In this subsection, we will get a sharper energy estimate for solutions obtained from the backward Euler method. Recall the backward Euler method for solving the heat equation is

U_jn+1_{− U}_jn= λ(U_j−1n+1_{− 2U}_jn+1+ U_j+1n+1) (6.24) whereλ = k/h2. An important technique is the summation by part:

X j (Uj− Uj−1)Vj = − X j Uj(Vj+1− Vj) (6.25)

(38)

We multiply both sides byU_jn+1, then sum overj. We get X j (U_jn+1)2_{− U}_jn+1U_jn_{= −λ}X j |Uj+1n+1− Ujn+1|2. The term U_jn+1U_jn_≤ 1 2((U n+1 j )2+ (Ujn)2)

by Cauchy-Schwartz. Hence, we get 1 2 X j (U_jn+1)2_{− (U}_jn)2)_{≤ −λ}X j |Uj+1n+1− Ujn+1|2 Or 1 2Dt,−kU n+1_k2_{≤ −}h2 k k h2kDx,+U n+1_{k = −kD} x,+Un+1k2. (6.26) where, D_t,−V_jn+1 := V n+1 j − Vjn k , Dx,+U n+1 j := U_j+1n+1_{− U}_jn+1 h ,

We sum inn from n = 1 to N , we get the following theorem.

Theorem 6.10. For the backward Euler method, we have the estimate

kUNk2+ C

N

X

n=1

kDx,+Unk2 ≤ kU0k2 (6.27)

This gives controls not only on_kUn_k2 but also on_kDx,+Unk.

Homeworks.

1. Show that the Crank-Nicolson method also has similar energy estimate.

2. Can forward Euler method have similar energy estimate?

2.7 Relaxation of errors

In this section, we want to study the evolution of an error. We consider

ut= uxx+ f (x) (7.28)

with initial dataφ. The error en

j := u(xj, tn) − Ujnsatisfies

(39)

2.7. RELAXATION OF ERRORS 35

We want to know how error is relaxed to zero from an initial errore0_{. We study the homogeneous}

finite difference quation first. That is

en+1_j = en_j + λ(e_j−1n _{− 2e}n_j + en_j+1). (7.30) oren+1 = G(un). The matrix is a tridiagonal matrix. It can be diagonalized by Fourier method. The eigenfunctions and eigenvalues are

vk,j = e2πijk/N, ρk = 1 − 2λ + 2λ cos(2πk/N) = 1 − 4λ sin2(πk/N ), k = 0, ..., N − 1.

When_{λ ≤ 1/2, all eigenvalues are negative except ρ}0:

1 = ρ0 > |ρ1| > |ρ2| > · · · .

The eigenfunction

v0 ≡ 1.

Hence, the projection of any discrete functionU onto this eigenfunction is the average:P_jUj.

Now, we decompose the error into

en=Xen_kvk Then en+1_k = ρkenk. Thus, en_k = ρn_ke0_k. We see thaten

kdecays exponentially fast excepten0, which is the average ofe0. Thus, the average of

initial error never decay unless we choose it zero. To guarantee the average ofe0 is zero, we may chooseU_jnto be the cell average ofu(x, tn) in the jth cell:

U_jn= 1 h

Z xj+1/2

x_j−1/2

u(x, tn) dx.

instead of the grid data. This implies the initial error has zero local averages, and thus so does the global average.

The contribution of the truncation to the true solution is:

en+1= ρkenk+ ∆tτkn Its solution is en_k = ρn_ke0_k+ ∆t n−1 X m=0 ρn−1−m_k τ_km

We see that the termen

0 does not tend to zero unlessτ0m = 0. This can be achieved if we choose Uj

(40)

Homeworks.

1. DefineUj := _h1

Rxj+1/2

x_j−1/2 u(x) dx. Show that if u(x) is a smooth periodic function on [0, 1],

then u′′(xj) = 1 h2(Uj−1− 2Uj+ Uj+1) + τ withτ = O(h2_).

2.8 Boundary Conditions

2.8.1 Dirichlet boundary condition Dirichlet boundary condition is

u(0) = a, u(1) = b. (8.31)

The finite difference approximation of uxx at x1 involves u at x0 = 0. We plug the boundary

contion: uxx(x1) = U0− 2U1+ U2 h2 + O(h 2_{) =} a − 2U1+ U2 h2 + O(h 2₎ Similarly, uxx(xN −1) = U_{N −2}_{− 2U}_{N −1}+ UN h2 + O(h 2_{) =} UN −2− 2UN −1+ b h2 + O(h 2₎

The unknowns areU₁n, ..., U_{N −1}n with_{N −1 finite difference equations at x}1, ..., xN −1. The discrete

Laplacian becomes A =      −2 1 _{0 · · · 0} 0 1 _{−2 1 · · · 0} 0 .. . ... ... . .. ... ... 0 0 _{0 · · · 1 −2}     . (8.32)

This discrete Laplacian is the same as a discrete Laplacian with zero Dirichlet boundary condition. We can have energy estimates, entropy estimates as the case of periodic boundary condition. Next, we exame how error is relaxed for the Euler method with zero Dirichlet boundary con-dition. From Fourier method, we observe that the eigenfunctions and eigenvalues for the forward Euler method are

vk,j= sin(2πjk/N ), ρk= 1 − 2λ + 2λ cos(2πk/N) = 1 − 4λ sin2(πk/N ), k = 1, ..., N − 1.

In the present case, all eigenvalues

ρi < 1, i = 1, ..., N − 1.

provided the stability condition

(41)

2.8. BOUNDARY CONDITIONS 37

Thus, the errorsen

i decays to zero exponentially for alli = 1, ..., N − 1. The slowest mode is ρ1

which is ρ1 = 1 − 4λ sin2(π/N ) ≈ 1 − 4 π N 2 and ρn₁ _≈ 1 − 4 π_N2 n ≈ e−4π2t where we have usedk/h2is fixed andnk = t.

2.8.2 Neumann boundary condition The Neumann boundary condition is

u′(0) = σ0, u′(1) = σ1. (8.33)

We may use the following disrete discretization methods:

• First order: U1− U0 h = σ0. • Second order-I: U1− U0 h = ux(x1/2) = ux(0) + h 2uxx(x0) = σ0+ h 2f (x0)

• Second order-II: we use extrapolation

3U0− 2U1+ U2

2h2 = σ0.

The knowns are U_jn with j = 0, ..., N . In the mean time, we add two more equations at the boundaries.

Homeworks.

1. Find the eigenfunctions and eigenvalues for the discrete Laplacian with the Neumann bound-ary condition (consider both first order and second order approximation at boundbound-ary). Notice that there is a zero eigenvalue.

Hint: You may use Matlab to find the eigenvalues and eigenvectors.

Here, I will provide another method. SupposeA is the discrete Laplacian with Neumann boundary condition. _{A is an (N + 1) × (N + 1) matrix. Suppose Av = λv. Then for j = 1, ..., N − 1, v} satisfies

(42)

Forv0, we have

−2v0+ 2v1= λv0.

ForvN, we have

−2vN + 2vN −1= λvN.

Then this matrix has the following eigenvectors:

v_jk= cos(πjk/N ), with eigenvalue λk_{= −2 + 2 cos(πk/N) = −4 sin}2 πk 2N . Homeworks.

1. Complete the calculation.

2. Consider

ut= uxx+ f (x)

on[0, 1] with Neumann boundary condition u′(0) = u′(1) = 0. IfR _{f (x) dx 6= 0. What wil} happen to_{u as t → ∞?}

2.9 The discrete Laplacian and its inversion

We consider the elliptic equation

uxx− αu = f (x), x ∈ (0, 1)

2.9.1 Dirichlet boundary condition Dirichlet boundary condition is

u(0) = a, u(1) = b (9.34)

The finite difference approximation of uxx at x1 involves u at x0 = 0. We plug the boundary

contion: uxx(x1) = U0− 2U1+ U2 h2 + O(h 2_{) =} a − 2U1+ U2 h2 + O(h 2₎ Similarly, uxx(xN −1) = U_{N −2}_{− 2U}_{N −1}+ UN h2 + O(h 2_{) =} UN −2− 2UN −1+ b h2 + O(h 2₎

(43)

2.9. THE DISCRETE LAPLACIAN AND ITS INVERSION 39

The unknowns areUn

1, ..., U_{N −1}n withN −1 finite difference at x1, ..., xN −1. The discrete Laplacian

becomes A =      −2 1 _{0 · · · 0} 0 1 _{−2 1 · · · 0} 0 .. . ... ... . .. ... ... 0 0 _{0 · · · 1 −2}     . (9.35)

This is the discrete Lalacian with Dirichlet boundary condition. In one dimension, we can solve A−1_{explicitly. Let us solve}_{(A − 2β)}−1_where_{β = αh}2_{/2. The difference equation}

U_j−1_{− (2 + 2β)U}j+ Uj+1 = 0

has two independent solutionsρ1andρ2, whereρiare roots of

ρ2_{− (2 + 2β)ρ + 1 = 0.} That is

ρ = 1 + β ±p(1 + β)2_{− 1}

Whenβ = 0, the two solutions are Uj = 1 and Uj = j. This gives the fundamental solution

Gi,j =

jCi j ≤ i

(N − j)Ci′ j ≥ i

FromG_i,i−1_−2Gi,i+ Gi,i+1 = 1 and iCi = (N −i)Ci′we getCi = −(N −i)/N and Ci′ = −i/N.

Whenβ > 0, the two roots are ρ1 < 1 and ρ2 > 1.

Homeworks.

1. Use matlab or maple to find the fundamental solutionGi,j := (A − 2β)−1withβ > 0.

2. Is it correct thatvi,jhas the following form?

Gi,j =

(

ρj−i₁ _{N − 1 > j ≥ i} ρj−i₂ 1 < j < i

Let us go back to the original equation:

uxx− αu = f (x)

The above study of the Green’s function of the discrete Laplacian helps us to quantify the the error produced from the source term. IfAu = f and A−1= G, then an error in f , say τ , will produce an error

(44)

If the off-diagonal part ofG decays exponentially (i.e. β > 0), then the error is “localized,” oth-erwise, it polutes everywhere. The error from the boundary also has the same behavior. Indeed, if β = 0, then The discrete solution is

u(xj) = aG0(j) + bG1(j) +

X

j

Gi,jfj

whereG(j) = jh, G1(j) = 1−jh and G = A−1, the Green’s function with zero Dirichlet boundary

condition. Here,G0 solves the equation

G0(i − 1) − 2G0(i) + G0(i + 1) = 0, i = 1, ..., N − 1,

for_{j = 1, ..., N − 1 with G}0(0) = 1 and G0(N ) = 0. And G1 solves the same equation with

G1(0) = 0 and G1(N ) = 1.

Ifβ > 0, we can see that both G0andG1are also localized.

Project 2. Solve the following equation

uxx− αu + f (x) = 0, x ∈ [0, 1]

numerically with periodic, Dirichlet and Neumann boundary condition. The equilibrium

1. A layer structure f (x) = −1 1/4 < x < 3/4 1 otherwise 2. An impluse f (x) = γ 1/2 − δ < x < 1/2 + δ 0 otherwise 3. A dipole f (x) =    γ _{1/2 − δ < x < 1/2} −γ 1/2 < x < 1/2 + δ 0 otherwise

You may chooseα = 0.1, 1, 10, observe how solutions change as you vary α. Project 3. Solve the following equation

−uxx+ f (u) = g(x), x ∈ [0, 1]

numerically with Neumann boundary condition. Here,f (u) = F′(u) and the potential is

F (u) = u4_{− γu}2.

Study the solution as a function ofγ. Choose simple g, say piecewise constant, a delta function, or a dipole.

(45)

Chapter 3

Finite Difference Methods for Linear

elliptic Equations

3.1 Discrete Laplacian in two dimensions

We will solve the Poisson equation

△u = f in a domain_{Ω ⊂ R}2with Dirichlet boundary condition

u = g on ∂Ω

Such a problem is a core problem in many applications. We may assumeg = 0 by substracting a suitable function fromu. Thus, we limit our discussion to the case of zero boundary condition. Let h be the spatial mesh size. For simplicity, let us assume Ω = [0, 1] × [0, 1]. But many discussion below can be extended to general smooth bounded domain.

3.1.1 Discretization methods

Centered finite difference The Laplacian is approximated by

A = 1

h2 (Ui−1,j + Ui+1,j+ Ui,j−1+ Ui,j+1− 4Ui,j) .

For the square domain, the indeces run from_{1 ≤ i, j ≤ N − 1 and} U0,j = UN,j = Ui,0= Ui,N = 0

from the boundary condition.

(46)

If we order the unknowns_{U by i + j ∗ (N − 1) with j being outer loop index and i the inner} loop index, then the matrix form of the discrete Laplacian is

A = 1 h2        T I I T I I T I . .. ... ... I T       

This is an_{(N − 1) × (N − 1) block tridiagonal matrix. The block T is an (N − 1) × (N − 1) matrix}

T =        −4 1 1 _{−4 −1} 1 ₋₄ 1 . .. ... ... 1 ₋₄       

Since this discrete Laplacian is derived by centered finite differencing over uniform grid, it is second order accurate, the truncation error

τi,j :=

1

h2 (u(xi−1, yj) + u(xi+1, yj) + u(xi, yj−1) + u(xi, yj+1) − 4u(xi, yj))

= O(h2).

3.1.2 The 9-point discrete Laplacian The Laplacian is approximated by

∇29 = 1 6h2   1 4 1 4 −20 4 1 4 1   One can show by Taylor expansion that

∇29u = ∇2u + 1 12h 2 ∇4u + O(h4). If_{u is a solution of ∇}2u = f , then ∇29u = f + 1 12h 2 ∇2f + O(h4). Thus, we get a 4th order method:

∇29Uij = fij+h 2

12∇

2_f ij

(47)

3.2. STABILITY OF THE DISCRETE LAPLACIAN 43

3.2 Stability of the discrete Laplacian

We have seen that the true solution of_{△u = f with Dirichlet boundary condition satisfies} Au = f + τ,

whereA is the discrete Laplacian and τ is the truncation error and satisfies τ = O(h2_{) in maximum}

norm. The numerical solutionU satisfies AU = f . Thus, the true error satisfies

Ae = τ,

where_{e = u − U. Thus, e satisfies the same equation with right-hand side τ and with the Dirichlet} boundary condition. To get the convergence result, we need an estimate ofe in terms of τ . This is the stability criterion of_{A. We say that A is stable if there exists some norm k · k and a constant C} such that

kek ≤ CkAek. 3.2.1 Fourier method

Since our domain_{Ω = [0, 1] × [0, 1] and the coefficients are constant, we can apply Fourier} trans-form. Let us see one dimensional case first. Consider the Laplaciand2/dx2on domain[0, 1] with Dirichlet boundary condition. The discrete Laplacian isA = _h12diag(1, −2, 1), where h = 1/N.

From

A sin(iπkh) = 1

h2 (sin((i + 1)πhk) + sin((i − 1)πhk) − 2 sin(iπhk))

= 2

h2 ((cos(πhk) − 1) sin(iπhk)) .

The above is valid fori = 1, ..., N − 1 and k = 1, ..., N − 1. We also see that sin(iπkh) sat-isfies Dirichlet boundary condition at i = 0 and i = N . we see that the eigenvectors of A are (sin(iπhk))N −1_i=1 for_{k = 1, ..., N − 1. The corresponding eigenvalues are} _h22(cos(πhk) − 1).

For two dimensional case, the eigenfunctions of the discrete Laplacian are

(Uk,ℓ)i,j = sin(iπkh) sin(jπℓh).

The corresponding eigenvalues are

λk,ℓ =

2

h2(cos(kπh) + cos(ℓπh) − 2)

= −_h4₂(sin2(kπh/2) + sin2(ℓπh/2)) The smallest eigenvalue (in magnitude) is

λ1,1 = −

8 h2 sin

2

(48)

To show the stability, we take Fourier transform ofU and A. We then have

h bA ˆU , ˆ_{U i}_{≥ 2π}2_{k ˆ}_{U k}2. The left-hand side has

h bA ˆU , ˆ_{U i}

≤ k bA ˆ_{U kk ˆ}_{U k.} Hence, theL2_{norm of b}_{A has the following estimate:}

k bA ˆ_{U k ≥ 2π}2_{k ˆ}_{U k.} Thus, we get

k ˆU k ≤ 1

2π2k bA ˆU k.

From Parseval equality, we have

kUk ≤ _2π1₂kAUk Applying this stability to the formula:Ae = τ , we get

kek ≤ _2π12kτk = O(h 2_).

Homeworks.

1. Compute th eigenvalues and eigenfunctions of the 9-point discrete Laplacian on the domain [0, 1] × [0, 1] with zero boundary condition.

3.2.2 Energy method

Below, we use energy method to prove the stability result for discrete Laplacian. We shall prive it for rectangular domain. However, it can be extended to more general domain. To perform energy estimate, we rewrite the discrete Laplacian as

AUi,j =

1

h2(Ui−1,j+ Ui+1,j+ Ui,j−1+ Ui,j+1− 4Ui,j) = ((Dx+Dx−+ Dy+Dy−)U )i,j

where

(Dx+U )i,j =

Ui+1,j− Ui,j

h

the forward differencing. We multiply the discrete Laplacian by Ui,j, then sum over all i, j. By

applying the ummation by part, we get

(AU, U ) = ((Dx+Dx−+ Dy+Dy−)U, U )

= −(Dx−U, Dx−U ) − (Dy−U, Dy−U )

(49)

3.2. STABILITY OF THE DISCRETE LAPLACIAN 45

Here, the discreteL2_{norm is defined by}

kUk2h=

X

i,j

|Ui,j|2h2.

The boundary term does not show up beause we consider the zero Dirichlet boundary problem. Thus, the discrete Poisson equation has the estimate

k∇hU k2h= |(f, U)| ≤ kf khkUkh. (2.1)

Next, for the zero Dirichlet boundary condition, we have Poincare inequality. Before stating the Poincare inequality, we need to clarify the meaning of zero boundary condition in the discrete sense. We define the Sobolev spaceH1

h,0to be the completion of the restriction of allC01functions

to the grid points under the discreteH1norm. Here,C₀1function is aC1function that is zero on the boundary; the discreteH1_{norm is}

kUkh,1:= kUkh+ k∇hU kh.

Lemma 2.2. LetΩ be a bounded domain in_R2_{, then there exist a constant}_d

Ω, which is the diameter

of the domain_{Ω, such that for any U ∈ H}_h,01 ,

kUkh≤ dΩk∇hU kh (2.2)

Proof. Let us take_{Ω = [0, X]×[0, Y ] as an example for the proof. We assume X = Mh, Y = Nh.} From zero boundary condition, we have

U_i,j2 = ( i X i′₌₁ D_x−Ui′_,jh)2 ≤ ( i X i′₌₁ 12_{) · (} i X i′₌₁ (D_x−Ui′_,j)2)h2(H ¨older’s inequa;ity) ≤ i( M X i′₌₁ (D_x−Ui′_,j)2)h2

multiply both sides byh2_{then sum over all}_{i, j, we get}

kUk2h = X i,j U_i,j2 h2 ≤ ( M X i=1 i)h2X i′_,j (D_x−Ui′_,j)2h2 ≤ M 2 2 h 2X i′_,j (D_x−Ui′_,j)2h2 = M 2 2 h 2 kDx−U k2h

(50)

Similarly, we have kUk2h ≤ N2 2 h 2_kD y−U k2h Thus, kUk2h ≤ h2 1 2max{M 2_{, N}2_}k∇ hU k2 ≤ d2Ωk∇hU k2h.

With the Poincare inequality, we can obtain two estimates forU .

Proposition 1. Consider the discrete Laplacian with zero boundary condition. We have

kUkh≤ d2Ωkf kh, (2.3)

k∇hU k ≤ dΩkf kh. (2.4)

Proof. From

k∇hU k2h ≤ kf kh· kUkh

We apply the Poincare inequality to the left-hand side, we obtain

kUk2h ≤ d2Ωk∇Uk2h≤ d2Ωkf khkUkh

This yields

kUkh ≤ d2Ωkf kh

If we apply the Poincare inequality to the right-hand side, we get

k∇hU k2h≤ kf kh· kUkh ≤ kf kh· dΩk∇hU kh

Thus, we obtain

k∇hU k ≤ dΩkf kh

When we apply this result toAe = τ , we get

kek ≤ d2Ωkτk = O(h2)

Finite Difference Methods for Solving Differential Equations