Some discretization methods

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

u_t= u_xx

Leth = ∆x, k = ∆t be the spatial and temporal mesh sizes. Define x_j = jh, j ∈ Z and tⁿ= nk, n ≥ 0. Let us abbreviate u(xj, tⁿ) by uⁿ_j. We shall approximateuⁿ_j byU_jⁿ, whereU_jⁿsatisfies some finite difference equations.

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

u_xx = u_j+1− 2uj + u_j−1

h² + O(h²) This results in the following systems of ODEs

U˙_j(t) = U_j+1(t) − 2Uj(t) + U_j−1(t) h²

or in vector form

U =˙ 1 h²AU whereU = (U₀, U₁, ...)^t,A = diag (1, −2, 1).

23

Homeworks.

1. Derive the 4th order centered finite difference approximation foru_xx: u_xx = 1

h²(−uj−2+ 16u_j−1− 30uj+ 16u_j+1− uj+2) + O(h⁴).

2. Derive a 2nd order centered finite difference approximation for(κ(x)u_x)_x. Temporal discretization We can apply numerical ODE solvers

• Forward Euler method:

Uⁿ⁺¹= Uⁿ+ k

h²AUⁿ (1.1)

• Backward Euler method:

Uⁿ⁺¹ = Uⁿ+ k

h²AUⁿ⁺¹ (1.2)

• 2nd order Runge-Kutta (RK2):

Uⁿ⁺¹− Uⁿ= k

h²AU^n+1/2, U^n+1/2= Uⁿ+ k

2h²AUⁿ (1.3)

• Crank-Nicolson:

Uⁿ⁺¹− Uⁿ= k

2h²(AUⁿ⁺¹+ AUⁿ). (1.4)

These linear finite difference equations can be solved formally as Uⁿ⁺¹= GUⁿ

where

• Forward Euler: G = 1 +_h^k2A,

• Backward Euler: G = (1 − _h^k2A)⁻¹,

• RK2: G = 1 + _h^k2A +¹_{2 h}^k2

2

A²

• Crank-Nicolson: G = ¹⁺

k 2h2A 1−_2h2^k A

For the Forward Euler, We may abbreviate it as

U_jⁿ⁺¹= G(U_j−1ⁿ , U_jⁿ, U_j+1ⁿ ), (1.5) where

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

h²(U_j−1− 2U^j+ Uj+1)

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_jⁿconverges tou(xj, tⁿ) 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

uⁿ⁺¹_j − uⁿj = k

h² uⁿ_j+1− 2uⁿj + uⁿ_j−1

+ kτ_jⁿ (1.6)

where

τ_jⁿ= D_t,+uⁿ_j − (ut)ⁿ_j − (D+D₋uⁿ_j − (uxx)ⁿ_j) = O(k) + O(h²).

Leteⁿ_j denote foruⁿ_j − Ujⁿ. Then substract (1.1) from (1.6), we get eⁿ⁺¹_j − eⁿj = k

h² eⁿ_j+1− 2eⁿj + eⁿ_j−1

+ kτ_jⁿ. (1.7)

This can be expressed in operator form:

eⁿ⁺¹ = Geⁿ+ kτⁿ. (1.8)

keⁿk ≤ kGeⁿ⁻¹k + kkτⁿ⁻¹k

≤ kG²eⁿ⁻²k + k(kGτⁿ⁻²k + kτⁿ⁻¹k)

≤ kGⁿe⁰k + k(kGⁿ⁻¹τ⁰k + · · · + kGτⁿ⁻²k + kτⁿ⁻¹k) Suppose G satisfies the stability condition

kGⁿU k ≤ CkUk for someC independent of n. Then

keⁿk ≤ Cke⁰k + C max_m |τ^m|.

If the local truncation error has the estimate

maxm kτ^mk = O(h²) + O(k) and the initial errore⁰satisfies

ke⁰k = O(h²), then so does the global true error satisfies

keⁿk = O(h²) + 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 kτh,kk → 0 as h, k → 0.

Definition 1.7. A finite difference schemeUⁿ⁺¹ = G_h,k(Uⁿ) is called stable under the norm k · k in a region(h, k) ∈ R if

kGⁿh,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

²

Stability – von Neumann Analysis

Since we only deal with smooth solutions in this section, theL²-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 giveL²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

kUk² =X

j

|Uj|² and its Fourier transform

U (ξ) =ˆ 1 2π

XU_je^−ijξ.

The advantages of Fourier method for analyzing finite difference scheme are

• the shift operator is transformed to a multiplier:

T U (ξ) = ed ^iξU (ξ),ˆ where(T U )_j := U_j+1;

• the Parseval equility

kUk² = k ˆU k²

≡ Z _π

−π| ˆU (ξ)|²dξ.

2.2. L²STABILITY – VON NEUMANN ANALYSIS 27 If a finite difference scheme is expressed as

U_jⁿ⁺¹= (GUⁿ)j = Xm i=−l

ai(TⁱUⁿ)j,

then U[ⁿ⁺¹= bG(ξ) cUⁿ(ξ).

From the Parseval equality,

kUⁿ⁺¹k² = k [Uⁿ⁺¹k²

= Z π

−π| bG(ξ)|²| cUⁿ(ξ)|²dξ

≤ max

ξ | bG(ξ)|² Z π

−π| cUⁿ(ξ)|²dξ

= | bG|²_∞kUk² Thus a sufficient condition for stability is

| bG|∞≤ 1. (2.9)

Conversely, suppose| bG(ξ₀)| > 1, from bG being a smooth function in ξ, we can find ǫ and δ such that

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

Let us choose an initial dataU₀inℓ²such that cU⁰(ξ) = 1 for |ξ − ξ0| ≤ δ. Then k cUⁿk² =

Z

| bG|²ⁿ(ξ)|cU⁰|²

≥ Z

|ξ−ξ0|≤δ| bG|²ⁿ(ξ)|cU⁰|²

≥ (1 + ǫ)²ⁿδ → ∞ 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_jⁿ⁺¹=

Xm k=−l

a_kU_j+kⁿ

with constant coefficients is stable if and only if G(ξ) :=b

Xm k=−l

a_ke^−ikξ satisfies

−π≤ξ≤πmax | bG(ξ)| ≤ 1. (2.10)

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_jⁿ⁺¹= αU_j−1ⁿ + βU_jⁿ+ γU_j+1ⁿ , (3.11) where

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

We multiply (3.11) byU_jⁿ⁺¹on both sides, apply Cauchy-Schwarz inequality, we get

(U_jⁿ⁺¹)² = αU_j−1ⁿ U_jⁿ⁺¹+ βU_jⁿU_jⁿ⁺¹+ γU_j+1ⁿ U_jⁿ⁺¹

≤ α

2((U_j−1ⁿ )²+ (U_jⁿ⁺¹)²) +β

2((U_jⁿ)²+ (U_jⁿ⁺¹)²) + γ

2((U_j+1ⁿ )²+ (U_jⁿ⁺¹)²) Here, we have usedα, β, γ ≥ 0. We multiply this inequality by h and sum it over j ∈ Z. Denote

kUk²:=



X

j

|U^j|²h





1/2

.

We get

kUⁿ⁺¹k² ≤ α

2(kUⁿk²+ kUⁿ⁺¹k²) + β

2(kUⁿk²+ kUⁿ⁺¹k²) +γ

2(kUⁿk²+ kUⁿ⁺¹k²)

= 1

2(kUⁿk²+ kUⁿ⁺¹k²).

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

kUⁿ⁺¹k² ≤ kUⁿk². (3.12)

Homeworks.

1. Can the RK-2 method possess an energy estimate?

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, U_j, U_j+1, provided 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 Uⁿ⁺¹ is a convex combination (average) of Uⁿ is very important. Given any convex functionη(u), called entropy function, by Jenson’s inequality, we get

η(U_jⁿ⁺¹) ≤ αη(U_j−1ⁿ ) + βη(U_jⁿ) + γη(U_j+1ⁿ ) (4.14)

This means that the “entropy” decreases in time. In particular, we choose

• η(u) = |u|², we recover theL²stability, the generalL^pstability. Takingp → ∞, we recover L^∞stability.

• η(u) = |u − c| for any constant c, we obtain X

j

|Ujⁿ⁺¹− c| ≤X

j

|Ujⁿ− 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)⁻¹ (5.16)

where

λ = k

h², A = diag(1, −2, 1).

The matrixM := I − λA has the following property:

m_ii> 0, m_ij ≤ 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λ

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+ 2u_j − auj+1= 0 (5.18)

has the form:

u_j = C₁ρ^j₁+ C₂ρ^j₂

2.5. ENTROPY ESTIMATE FOR BACKWARD EULER METHOD 31 whereρ₁, ρ₂are the characteristic roots, i.e. the roots of the polynomial equation

−aρ²+ 2ρ − a = 0.

Now, we define a fundamental solution:

gj = (

ρ^j₁ forj ≥ 0 ρ^j₂ forj < 0 .

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

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

Thus,M⁻¹ = (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 Uⁿ⁺¹= (1 − λA)⁻¹= GUⁿ,

with initial dataU_j⁰ = δ_j. This finite difference equation is conservative, namely,P

jU_jⁿis inde-pendent. This can easily be checked. Thus, we getP

jg_j =P

jδ_j = 1. With this and the positivity ofg_j, we can thinkU_jⁿ⁺¹is a convex combination ofU_jⁿwith weights g_j. 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_jⁿbe a solution of the difference equation derived from the backward Euler method for the heat equation. Then we have

X

j

η(U_jⁿ) ≤X

j

η(U_j⁰). (5.19)

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

H^m := {u : R → R|u, u^′, ..., u^(m) ∈ L²(R)}

The discrete Sobolev space for functions defined on gridG_h := {jh|j ∈ Z}.

H_h^m:= {U : Gh → R|U, Dx,+U, ..., D_x,+^m U ∈ ℓ²}.

Here,(D_x,+U )ⁿ_j := (U_j+1ⁿ − Ujⁿ)/h

For any discrete functionUj ∈ Hh^mwe can construct a functionu in H^mdefined by u(x) :=X

j

U_jφ_h(x − xj) (6.20)

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

u_h(x_j) = U_j It can be shown that

kDx^mu_hk ≡ kD^mx,+U k. (6.21)

Similarly, the spaceL^∞_k (H_h^m) can be embeded into L^∞(H^m) by defining u_h,k(x, t) =X

n≥0

X

j

U_jⁿφ_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 equationu_t= u_xxreads U_jⁿ⁺¹= U_jⁿ+ k

h² U_j−1ⁿ − 2Ujⁿ+ U_j+1ⁿ
.

2.6. EXISTENCE THEORY 33 Here, We have seen that we can get the energy estimate:

kUⁿk ≤ kU⁰k.

We perform finite difference operation on the above equation, say the forward Euler equation, for instance, letV_jⁿ = (D_x,+U )ⁿ_j := (U_j+1ⁿ − Ujⁿ)/h. Then V_jⁿ satisfies the same finite difference equation

V_jⁿ⁺¹ = V_jⁿ+ k

h² V_j−1ⁿ − 2Vjⁿ+ V_j+1ⁿ
.

Thus, it also possesses the same energy estimate. Similar estimate forD_x,+² U . In general, we have kDx,+^m Uⁿk ≤ kD^mx,+U⁰k. (6.22) If we assume the initial dataf ∈ H², then we getUⁿ∈ Hh²for alln ≥ 0.

Theorem 6.9. If the initial datau₀ ∈ H^m,m ≥ 2 and k/h² ≤ 1/2, then the solution of forward Euler equation has the estimate

kDx,+^m Uⁿk ≤ kD^mx,+U⁰k, kDt,+Uⁿk ≤ kDx,+² U⁰k (6.23) Further, the corresponding smoothing functionu_h,k has the same estimate and has a subsequence converges to a solutionu(x, t) of the original equation.

Proof. The functionsu_h,kare unformly bounded in W^1,∞(H²). Hence they have a subsequence converges to a function u ∈ W^1,∞(H²) weakly in W^1,∞(H²) and strongly in L^∞(H¹). The functionsu_h,ksatisfy

u_h,k(x_j, tⁿ⁺¹) − uh,k(x_j, tⁿ) = k

h²(u_h,k(x_j−1, tⁿ) − 2uh,k(x_j, tⁿ) + u_h,k(x_j+1, tⁿ)) Multiply a test smooth function φ, sum over j and n, take summation by part, we can get the subsubsequence converges to a solution ofu_t= u_xxweakly.

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_jⁿ⁺¹− Ujⁿ= λ(U_j−1ⁿ⁺¹− 2Ujⁿ⁺¹+ U_j+1ⁿ⁺¹) (6.24) whereλ = k/h². An important technique is the summation by part:

X

j

(U_j− Uj−1)V_j = −X

j

U_j(V_j+1− Vj) (6.25)

There is no boundary term because we consider periodic condition in the present case.

We multiply both sides byU_jⁿ⁺¹, then sum overj. We get X

j

(U_jⁿ⁺¹)²− Ujⁿ⁺¹U_jⁿ= −λX

j

|Uj+1ⁿ⁺¹− Ujⁿ⁺¹|².

The term

U_jⁿ⁺¹U_jⁿ≤ 1

2((U_jⁿ⁺¹)²+ (U_jⁿ)²) by Cauchy-Schwartz. Hence, we get

1 2

X

j

(U_jⁿ⁺¹)²− (Ujⁿ)²)

≤ −λX

j

|Uj+1ⁿ⁺¹− Ujⁿ⁺¹|²

Or

1

2D_t,−kUⁿ⁺¹k²≤ −h² k

k

h²kDx,+Uⁿ⁺¹k = −kDx,+Uⁿ⁺¹k². (6.26) where,

D_t,−V_jⁿ⁺¹ := V_jⁿ⁺¹− Vjⁿ

k , Dx,+U_jⁿ⁺¹:= U_j+1ⁿ⁺¹− Ujⁿ⁺¹

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

kU^Nk²+ C XN n=1

kD^x,+Uⁿk² ≤ kU⁰k² (6.27)

This gives controls not only onkUⁿk² but also onkD^x,+Uⁿk.

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

u_t= u_xx+ f (x) (7.28)

with initial dataφ. The error eⁿ_j := u(x_j, tⁿ) − Ujⁿsatisfies

eⁿ⁺¹_j = eⁿ_j + λ(eⁿ_j−1− 2eⁿj + eⁿ_j+1) + kτ_jⁿ (7.29)

2.7. RELAXATION OF ERRORS 35 We want to know how error is relaxed to zero from an initial errore⁰. We study the homogeneous finite difference quation first. That is

eⁿ⁺¹_j = eⁿ_j + λ(eⁿ_j−1− 2eⁿj + eⁿ_j+1). (7.30) oreⁿ⁺¹ = G(uⁿ). The matrix is a tridiagonal matrix. It can be diagonalized by Fourier method.

The eigenfunctions and eigenvalues are

v_k,j = e^2πijk/N, ρ_k = 1 − 2λ + 2λ cos(2πk/N) = 1 − 4λ sin²(πk/N ), k = 0, ..., N − 1.

Whenλ ≤ 1/2, all eigenvalues are negative except ρ⁰:

1 = ρ0 > |ρ¹| > |ρ²| > · · · . The eigenfunction

v₀ ≡ 1.

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

jU_j. Now, we decompose the error into

eⁿ=X eⁿ_kv_k Then

eⁿ⁺¹_k = ρ_keⁿ_k. Thus,

eⁿ_k = ρⁿ_ke⁰_k.

We see thateⁿ_kdecays exponentially fast excepteⁿ₀, which is the average ofe⁰. Thus, the average of initial error never decay unless we choose it zero. To guarantee the average ofe⁰ is zero, we may chooseU_jⁿto be the cell average ofu(x, tⁿ) in the jth cell:

U_jⁿ= 1 h

Z _xj+1/2

x_j−1/2

u(x, tⁿ) 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:

eⁿ⁺¹= ρ_keⁿ_k+ ∆tτ_kⁿ Its solution is

eⁿ_k = ρⁿ_ke⁰_k+ ∆t

n−1X

m=0

ρ^n−1−m_k τ_k^m

We see that the termeⁿ₀ does not tend to zero unlessτ₀^m = 0. This can be achieved if we choose U_j as well asf_jto be the cell averages instead the grid data.

Homeworks.

1. DefineU_j := _h¹R_xj+1/2

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

u^′′(x_j) = 1

h²(U_j−1− 2Uj+ U_j+1) + τ withτ = O(h²).

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:

u_xx(x₁) = U₀− 2U1+ U₂

h² + O(h²) = a − 2U1+ U₂

h² + O(h²) Similarly,

u_xx(x_{N −1}) = U_{N −2}− 2UN −1+ U_N

h² + O(h²) = U_{N −2}− 2UN −1+ b

h² + O(h²)

The unknowns areU₁ⁿ, ..., U_{N −1}ⁿ withN −1 finite difference equations at x1, ..., x_{N −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

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

In the present case, all eigenvalues

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

provided the stability condition

λ ≤ 1/2.

2.8. BOUNDARY CONDITIONS 37 Thus, the errorseⁿ_i decays to zero exponentially for alli = 1, ..., N − 1. The slowest mode is ρ1

which is

ρ₁ = 1 − 4λ sin²(π/N ) ≈ 1 − 4 π N

2

and

ρⁿ₁ ≈



1 − 4 π N

2n

≈ e^−4π2t where we have usedk/h²is fixed andnk = t.

2.8.2 Neumann boundary condition

The Neumann boundary condition is

u^′(0) = σ₀, u^′(1) = σ₁. (8.33) We may use the following disrete discretization methods:

• First order:

U₁− U0

h = σ₀.

• Second order-I:

U₁− U0

h = u_x(x_1/2) = u_x(0) + h

2u_xx(x₀) = σ₀+h 2f (x₀)

• Second order-II: we use extrapolation

3U₀− 2U1+ U₂ 2h² = σ0.

The knowns are U_jⁿ 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

v_j−1− 2vj+ v_j+1= λv_j, j = 1, ..., N − 1.

Forv₀, we have

−2v⁰+ 2v1= λv0. Forv_N, we have

−2vN + 2v_{N −1}= λv_N. Then this matrix has the following eigenvectors:

v_j^k= cos(πjk/N ), with eigenvalue

λ^k= −2 + 2 cos(πk/N) = −4 sin²

π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 tou as t → ∞?

2.9 The discrete Laplacian and its inversion

We consider the elliptic equation

u_xx− α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 u_xx at x₁ involves u at x₀ = 0. We plug the boundary contion:

u_xx(x₁) = U₀− 2U1+ U₂

h² + O(h²) = a − 2U1+ U₂

h² + O(h²) Similarly,

u_xx(x_{N −1}) = U_{N −2}− 2UN −1+ U_N

h² + O(h²) = U_{N −2}− 2UN −1+ b

h² + O(h²)

2.9. THE DISCRETE LAPLACIAN AND ITS INVERSION 39 The unknowns areU₁ⁿ, ..., U_{N −1}ⁿ withN −1 finite difference at x1, ..., x_{N −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⁻¹explicitly. Let us solve(A − 2β)⁻¹whereβ = αh²/2. The difference equation

U_j−1− (2 + 2β)Uj+ U_j+1 = 0 has two independent solutionsρ1andρ2, whereρiare roots of

ρ²− (2 + 2β)ρ + 1 = 0.

That is

ρ = 1 + β ±p

(1 + β)²− 1

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

 jCi j ≤ i

(N − j)Ci^′ j ≥ i

FromG_i,i−1−2Gi,i+ G_i,i+1 = 1 and iC_i = (N −i)Ci^′we getC_i = −(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 solutionG_i,j := (A − 2β)⁻¹withβ > 0.

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

G_i,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⁻¹= G, then an error in f , say τ , will produce an error

e = Gτ.

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(x_j) = aG₀(j) + bG₁(j) +X

j

G_i,jf_j

whereG(j) = jh, G₁(j) = 1−jh and G = A⁻¹, the Green’s function with zero Dirichlet boundary condition. Here,G₀ solves the equation

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

forj = 1, ..., N − 1 with G0(0) = 1 and G₀(N ) = 0. And G₁ solves the same equation with G1(0) = 0 and G1(N ) = 1.

Ifβ > 0, we can see that both G₀andG₁are also localized.

Project 2. Solve the following equation

u_xx− α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) = u⁴− γu².

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

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²with 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

h² (U_i−1,j + U_i+1,j+ U_i,j−1+ U_i,j+1− 4Ui,j) . For the square domain, the indeces run from1 ≤ i, j ≤ N − 1 and

U_0,j = U_N,j = U_i,0= U_i,N = 0 from the boundary condition.

41

If we order the unknownsU 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

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

τ_i,j := 1

h² (u(x_i−1, y_j) + u(x_i+1, y_j) + u(x_i, y_j−1) + u(x_i, y_j+1) − 4u(xi, y_j))

= O(h²).

3.1.2 The 9-point discrete Laplacian

The Laplacian is approximated by

∇²9 = 1

One can show by Taylor expansion that

∇²9u = ∇²u + 1

12h²∇⁴u + O(h⁴).

Ifu is a solution of ∇²u = f , then

∇²9u = f + 1

12h²∇²f + O(h⁴).

Thus, we get a 4th order method:

∇²9Uij = fij+h² 12∇²fij

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(h²) in maximum norm. The numerical solutionU satisfies AU = f . Thus, the true error satisfies

Ae = τ,

wheree = 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 ofA. 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 Laplaciand²/dx²on domain[0, 1] with Dirichlet boundary condition. The discrete Laplacian isA = _h¹2diag(1, −2, 1), where h = 1/N.

From

A sin(iπkh) = 1

h² (sin((i + 1)πhk) + sin((i − 1)πhk) − 2 sin(iπhk))

= 2

h² ((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 fork = 1, ..., N − 1. The corresponding eigenvalues are _h²²(cos(πhk) − 1).

For two dimensional case, the eigenfunctions of the discrete Laplacian are (U^k,ℓ)_i,j = sin(iπkh) sin(jπℓh).

The corresponding eigenvalues are λ_k,ℓ = 2

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

= − 4

h²(sin²(kπh/2) + sin²(ℓπh/2)) The smallest eigenvalue (in magnitude) is

λ_1,1 = − 8

h² sin²(πh/2) ≈ −2π².

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

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



 ≤ k bA ˆU kk ˆU k.

Hence, theL²norm of bA has the following estimate:

k bA ˆU k ≥ 2π²k ˆU k.

Thus, we get

k ˆU k ≤ 1

2π²k bA ˆU k.

From Parseval equality, we have

kUk ≤ 1

2π²kAUk Applying this stability to the formula:Ae = τ , we get

kek ≤ 1

2π²kτk = O(h²).

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

AU_i,j = 1

h²(U_i−1,j+ U_i+1,j+ U_i,j−1+ U_i,j+1− 4Ui,j) = ((D_x+D_x−+ D_y+D_y−)U )_i,j where

(D_x+U )_i,j = U_i+1,j− Ui,j

h

the forward differencing. We multiply the discrete Laplacian by U_i,j, then sum over all i, j. By applying the ummation by part, we get

(AU, U ) = ((D_x+D_x−+ D_y+D_y−)U, U )

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

= −k∇hU k²h

3.2. STABILITY OF THE DISCRETE LAPLACIAN 45 Here, the discreteL²norm is defined by

kUk²h=X

i,j

|Ui,j|²h².

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 k²h= |(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 spaceH_h,0¹ to be the completion of the restriction of allC₀¹functions to the grid points under the discreteH¹norm. Here,C₀¹function is aC¹function that is zero on the boundary; the discreteH¹norm is

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

Lemma 2.2. LetΩ be a bounded domain inR², then there exist a constantd_Ω, which is the diameter of the domainΩ, such that for any U ∈ Hh,0¹ ,

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,j² = (

multiply both sides byh²then sum over alli, j, we get kUk²h = X

Similarly, we have

kUk²h ≤ N²

2 h²kDy−U k²h

Thus,

kUk²h ≤ h²1

2max{M², N²}k∇hU k²

≤ d²Ωk∇hU k²h.

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

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

kUkh≤ d²Ωkf kh, (2.3)

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

Proof. From

k∇hU k²h ≤ kf kh· kUkh

We apply the Poincare inequality to the left-hand side, we obtain kUk²h ≤ d²Ωk∇Uk²h≤ d²Ωkf khkUkh

This yields

kUkh ≤ d²Ωkf kh

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

k∇hU k²h≤ 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 ≤ d²Ωkτk = O(h²) k∇hek ≤ dΩkτk = O(h²).

Chapter 4

Finite Difference Theory For Linear Hyperbolic Equations

4.1 A review of smooth theory of linear hyperbolic equations

Hyperbolic equations appear commonly in physical world. The propagation of acoustic wave, electric-magnetic waves, etc. obey hyperbolic equations. Physical characterization of hyperbol-icity is that the signal propagates at finite speed. Mathematically, it means that compact-supported initial data yield compact-supported solutions at all time. This hyperbolicity property has been char-acterized in terms of coefficients of the corresponding linear partial differential equations through Fourier method.

They are two techaniques for hyperbolic equations, one is based on Fourier method (Garding et

They are two techaniques for hyperbolic equations, one is based on Fourier method (Garding et

在文檔中 Finite Difference Methods for Solving Differential Equations (頁 27-0)