5 Numerical Solution of Elliptic Boundary Value Problems

5 Numerical Solution of Elliptic Boundary Value Problems

5.1 Some Theory and Examples The Poisson equation

− ∆u = f in Ω ⊆ R

with f ∈ C

(Ω) (P) is widely used and solved in applications. We first consider the Dirichlet boundary condition

u = g auf Γ := ∂Ω mit g ∈ C

(Γ) (D) and make some assumptions on Ω:

Ω is bounded,

Γ = ∂Ω is piecewise smooth

Ω lies locally on one side of Γ.

Under these conditions there exists an unique solution u ∈ C

(Ω) ∩ C

(Ω), such that (P) and (D) are satisfied.

The solution depends continuously on the data: u

, u

are solutions of

−∆u = f with the boundary values g

and g

, then sup

x∈Ω

|u

(x) − u

(x)| = sup

x∈Γ

|g

(x) − g

(x)|

(continuous dependence on boundary data).

Further, there is a constant c, which only depends on Ω, such that sup

x∈Ω

|u(x)| ≤ sup

x∈Γ

|u(x)| + c sup

x∈Ω

|∆u(x)|

for all u ∈ C

(Ω) ∩ C

(Ω).

Why the latter inequality implies the continuous dependence of the

solution on the right hand side f ?

For (P) with Neumann boundary conditions

∂

u(x) = h(x) on Γ = ∂Ω with h ∈ C

(Γ) (N) one must require

Z

Γ

h(x)dS + Z

Ω

f(x)dx = 0

in order to ensure the existence of a solution.

The reason is the Gaussian

formula Z

Ω

∆u(x) dx = Z

Γ

∂

u(x) dS.

If u is a solution, so is u + c (c ∈ R).

1

Carl Friedrich Gauß (1777–1855)

In applied analysis, there are many methods for solving the Poisson equation in special cases. We restrict ourselves to a few examples which are useful for tests.

Example 1. Let Ω be the open unit disc in R

. If the Dirichlet boundary condition is given in a form of a Fourier

series

g(cos(ϕ), sin(ϕ)) = α

+

∞

X

j=1

[α

cos(jϕ) + β

sin(jϕ)],

then the solution of ∆u = 0 in Ω can be expressed as g(r cos(ϕ), r sin(ϕ)) = α

+

∞

X

j=1

r

[α

cos(jϕ) + β

sin(jϕ)].

2

Jean Baptiste Joseph Fourier (1768–1813)

The solution of ∆u = 0 with Neumann boundary condition (N) has the same form as

h(cos(ϕ), sin(ϕ)) =

∞

X

j=1

[α

cos(jϕ) + β

sin(jϕ)]

(Note R

−π

h(ϕ)dϕ = 0).

The constant α

∈ R in the expression for g can be chosen arbitrarily.

−3 −2 −1 0 1 2 3

−8

−6

−4

−2 0 2 4 6 8

r=1 r=.75 r=.5 r=.25

−1

−0.5 0

0.5 1

−1

−0.5 0 0.5 1

−8

−6

−4

−2 0 2 4 6 8

Example 2. Let g(x) = P

j=1

α

sin(jπx) be the Fourier expansion of a function g ∈ C

[0, 1].

The Laplace equation in the unit square Ω = (0, 1)

with boundary values

u(x, 0) = u(0, y) = u(1, y) = 0 und u(x, 1) = g(x) has the solution

u(x, y) =

∞

X

j=1

α

sinh(jπ) sin(jπx) sinh(jπy),

in Ω \ {(x, y) : y = 1} from C

.

Further examples with inhomogeneous Dirichlet boundary data also on

the other boundary segments can be constructed by superposition.

Example 3. We consider the Poisson equation

−∆u = f auf Ω = (−1, 1)

with the source term f (x, y) ≡ 1 and a homogeneous Dirichlet boundary condition applied at ∂Ω.

By separation of variables we construct the series solution u(x, y) = 1 − x

2 − 16 π

X

k∈N, kodd

sin(kπ(1 + x)/2) k

sinh(kπ)

·

sinh kπ(1 + y)

2 + sinh kπ(1 − y) 2

.

Perform the calculation via separation of variables for Example 2.

Transform Example 3 into a BVP for the Laplace equation with

inhomogeneous Dirichlet boundary conditions.

−1

−0.5 0

0.5 1

−1

−0.5 0 0.5 1

−0.05 0 0.05 0.1 0.15 0.2 0.25 0.3

−1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1

−1

−0.8

−0.6

−0.4

−0.2 0 0.2 0.4 0.6 0.8 1

Surface and contour plot of the solution from Example 3.

5.2 Laplace’s Equation on a Square 5.2.1 The Dirichlet Problem

We consider the following boundary value problem (BVP) for Laplace’s equation:

∆u = u

+ u

= 0 on Ω := (0, 1)

, (1a)

u = g along Γ := ∂Ω. (1b)

with a given boundary function g = g(x, y) defined on Γ.

Just as in our discretization of ordinary differential equations, we shall

apply finite difference methods to seek approximations to the solution

u = u(x, y) of (1) not at all points of Ω, but rather at a finite number

of grid points.

Consider each coordinate separately: introducing the partitions 0 = x

< x

< x

< · · · < x

< x

= 1, 0 = y

< y

< y

< · · · < y

< y

= 1 of [0, 1], we define the grid (tensor product grid) or mesh

Ω

:=

(x

, y

) : 0 ≤ i ≤ n

+ 1, 0 ≤ j ≤ n

+ 1 consisting of (n

+ 2)(n

+ 2) points.

A grid which is equidistant (in both directions) is one in which the grid spacing in each direction is constant:

x

= i∆x, i = 0, 1, . . . , n

+ 1, ∆x = 1/(n

+ 1), y

= j∆y, j = 0, 1, . . . , n

+ 1, ∆y = 1/(n

+ 1).

For equal spacing in both directions (uniform grid) we set

h := ∆x = ∆y.

0 0.2 0.4 0.6 0.8 1 0

0.2 0.4 0.6 0.8 1

tensor product grid

0 0.2 0.4 0.6 0.8 1

0 0.2 0.4 0.6 0.8 1

equidistant grid

In the following we consider the uniform case ∆x = ∆y = h (uniform grid):

0 0.2 0.4 0.6 0.8 1

0 0.2 0.4 0.6 0.8 1

Notation:

u

:= u(x

, y

) function value of exact solution at point (x

, y

) U

≈ u

approximation of u(x

, y

).

Taylor-expansion at (x

, y

):

u

= u

+

hu

+ h

2 u

+ h

6 u

+ h

24 u

i,j

+ O(h

),

u

= u

+

−hu

+ h

2 u

− h

6 u

+ h

24 u

i,j

+ O(h

).

Summing yields (the terms with h

also cancel) u

+ u

= 2u

+

h

u

+ h

12 u

i,j

+ O(h

).

Analogously, in the y-direction, u

+ u

= 2u

+

h

u

+ h

12 u

i,j

+ O(h

)

and, summing all 4 expansions and dividing by h

,

ui+1,j+ui−1,j+ui,j+1+ui,j−1−4ui,j

h²

| {z }

=:(∆_hu)_i,j

=

∆u+h²

12(uxxxx+uyyyy)

i,j

+O(h⁴).

The difference formula ∆

is called the 5-point-stencil or 5-point discrete Laplacian.

Summary: In each grid point (x

, y

), ∆

approximates the Laplace operator with a local discretisation error of order O(h

) as h → 0, i.e.,

(∆

u)

= (∆u)

+ O(h

).

Note: This requires sufficient smoothness of the solution, in this case

u ∈ C

(Ω).

In analogy with the differential equation (1a) we now require that the approximate solutions U

satisfy the discrete Laplace equation

∆

U

:= (∆

U )

= 0, 1 ≤ i ≤ n

, 1 ≤ j ≤ n

. (2a) (2a) represents n

n

equations in (n

+ 2)(n

+ 2) unknowns, since the 5-point stencil is not applicable in the boundary points.

But there, the boundary conditions supply the missing equations:

U

= g(x

, y

) if (x

, y

) ∈ ∂Ω. (2b) Written out in detail, (2b) reads

U

= g(0, y

), U

= g(1, y

), 0 ≤ j ≤ n

+ 1,

U

= g(x

, 0), U

= g(x

, 1), 0 ≤ i ≤ n

+ 1.

Determine the linear system for the numerical solutions U

of the elliptic BVP

−∆u(x, y) = 0,

u(x, 0) = u(0, y) = 0, u(x, 1) = x

, u(1, y) = y,

on a uniform grid with the following enumeration of the unknowns:

0 1

1 x

y

U1,1 U1,2 U1,3

U2,1 U2,2 U2,3

U3,1 U3,2 U3,3

As in our example, equations (2) together represent a linear system of equations

Au = b (3)

to determine the approximate values U

at the interior grid points.

The entries of the coefficient matrix A and the right hand side b depend on our enumeration of the unknowns. A common one is the lexicographic order:

u = [U

, U

, . . . , U

, U

, . . . , U

, . . . , U

, . . . , U

]

or

u =









 u

u

.. . u











, where u

=









 U

U

.. . U











, j = 1, 2, . . . , n

.

The Matrix A possesses the block structure

A = 1 h















 T I

I T I

I . .. ...

. .. I

I T

















∈ R

,

in which I denotes the n

× n

identity matrix and T is the tridiagonal matrix

T =

















−4 1 1 −4 1

1 . .. ...

. .. 1

1 −4

















∈ R

.

For unknowns U

near the boundary, i.e.,

i ∈ {1, n

} or j ∈ {1, n

},

the 5-point stencil (2a) contains neighboring approximation values already specified by the Dirichlet boundary condition (2b).

Moving these to the r.h.s. of each corresponding equation, we obtain for the (nonzero) entries of the vector b in (3) for the example shown

f

= − 1

h

g(x

, y

) + g(x

, y

)

or f

= − 1

h

g(x

, y

),

Representation of A via Kronecker products

Due to the

simple geometry of Ω,

the simple structure of the Laplacian and

the same type of boundary conditions on all boundaries,

the block tridiagonal matrix A possesses additional structure: it can be built up from discretization matrices arising in the simpler

one-dimensional BVP.

We therefore consider the discretization of the (ordinary) one-dimensional BVP

u

(x) = 0, u(0) = g

, u(1) = g

using central differences with uniform mesh width h = 1/(n + 1).

Its discrete approximation leads to the linear system of equations A

u = b

, A ∈ R

, b

∈ R

,

where

A

= 1

h

tridiag(1, −2, 1) := 1 h













−2 1 1 . .. ...

. .. 1

1 −2











 ,

b

= − 1

h

[g

, 0, . . . , 0, g

]

and vector of unknowns

u = [U

, U

, . . . , U

]

, U

≈ u(x

), i = 1, 2, . . . , n.

The discretization of the two-dimensional Laplace operator can now be expressed as

A = T

⊗ I + I ⊗ T

(4)

with T

=

tridiag(−1, 2, −1) and the identity matrix I

∈ R n × n.

In equation (4), the Kronecker product (or tensor product) M ⊗ N of two matrices M ∈ R

and N ∈ R

is defined by

M ⊗ N =







m

N . . . m

N

.. . .. .

m

N . . . m

N





 ∈ R

5.2.2 The Neumann Problem

We now consider a BVP for the Laplace operator in which the Dirichlet BC are replaced by the Neumann BC

∂u

∂n = h(x), x ∈ Γ. (5)

Since in our model problem the four segments of the boundary Γ lie parallel to the coordinate axes, the discretization of (5) is easy.

In the boundary points (x

, y

), for example, one may use either U

− U

h = h(x

, y

), (backward difference), (6a) U

− U

2h = h(x

, y

) (central difference) (6b)

with analogous formulas at the remaining boundaries.

Note:

The central difference formula (6b) introduced so-called ghost points U

, which, strictly speaking, lie outside the domain Ω.

These can, however, immediately be eliminated from the equations by solving (6b) for U

and inserting this expression into the 5-point stencil centered at (x

, y

).

In contrast to the discretized Dirichlet problem with the same mesh width h, for the Neumann problem the unknowns U

on the domain boundary Γ are not fixed by the boundary conditions alone.

They but must also be determined, along with the interior unknowns, by solving the coupled linear system of equations.

In this case the system has (n

+ 2) · (n

+ 2) unknowns.

The discretization matrix (same ordering) is now obtained as the Kronecker product

A = ˜ T

⊗ I + I ⊗ T ˜

with

T ˜

= 1 h



















−1 1 1 −2 . ..

. ..

. .. −2 1 1 −1

















 .

For illustration discretize the one-dimensional BVP

u

(x) = x (0 < x < 1), u(0) = 0, u

(1) = 1

with h =

once using backward differences and once using central

differences for the Neumann conditions.

5.2.3 Eigenvalues and Eigenvectors

One reason why (1) is often chosen as a model problem is that the spectral decomposition of the discretization matrix is available in closed form.

Beginning with the n × n matrix T

=

tridiag(1, −2, 1), we have T

v

= λ

v

, j = 1, 2, . . . , n,

with eigenvalues λ

= 2

h

[cos(jπh) − 1] = − 4

h

sin

jπh

2 , j = 1, 2, . . . , n, (7) and (orthonormal) eigenvectors

[v

]

= r 2

n + 1 sin(jkπh), k = 1, 2, . . . , n; j = 1, 2, . . . , n.

Using properties of the Kronecker product of matrices, the eigenvalues and eigenvectors of the 2D problem can also be inferred from those of the 1D problem.

More precisely: the eigenvalues λ

of A in (4) are given by

λ

= λ

+ λ

, 1 ≤ i, j ≤ n, (8) with associated eigenvectors

v

= v

⊗ v

, 1 ≤ i, j ≤ n.

Determine the smallest and the largest eigenvalue of the discretization

matrix A for the 2D Dirichlet case and discuss the conditioning of A in

dependence on the mesh width h.

5.2.4 Stability and Convergence

As we have already shown in the derivation of the 5-point approximation

∆

of the Laplace operator, there holds

∆

u

= ∆u

+ h

12 (u

+ u

)

+ O(h

).

Setting

[d

]

:= ∆

u

− ∆u

(local discretization error) [e

]

:= u

− U

(global discretization error) we conclude that, due to ∆u = 0 for the exact solution u of the BVP, e

solves the linear system of equations

A

e

= d

= O(h

)

with A

= A the matrix in (4).

We therefore have e

= A

d

and, because of ke

k = kA

d

k ≤ kA

kkd

k

the O(h

) behavior as h → 0 is transferred from the local to the global discretization error if kA

k remains uniformly bounded for all

sufficiently small values of h > 0.

Fixing the norm to be the Euclidean norm k · k = k · k

, we have kA

k = 1

|λ

(A

)| .

The uniform boundedness property can now be inferred from (7) and (8), since

|λ

(A

)| = 4

h

· 2 sin

πh

2 = 2π

+ O(h

), (h → 0).

Summary

The (global) discretization error of the approximate solution of the model problem (1) obtained by finite difference approximation using central differences with uniform mesh size h satisfies

ke

k = O(h

) as h → 0,

assuming the solution u is four times continuously differentiable.

The uniform boundedness of kA

k as h → 0 is a stability property of the difference scheme.

Check the O(h

) behaviour of the local discretization error with

Matlab for the test problem on page 180 with g(x) = 2 sin(3πx).

5.3 The Poisson Equation

Replacing the Laplace equation (1a) in the BVP (1) by the inhomogeneous (Poisson) equation

∆u = f,

with a given function f = f (x, y) defined on Ω, applying the central finite difference discretization leads to a discrete problem in which (2a) is modified to

∆

U

= f

, with f

= f (x

, y

) (1 ≤ i ≤ n

, 1 ≤ j ≤ n

).

For the Dirichlet problem, the right hand side f of the resulting linear

system of equations (3) consists of the sum of the Dirichlet boundary

values and the corresponding function values f

.

Excursus: The 9-Point Stencil

The formulas for difference schemes such as ∆

u can become cumbersome, particularly in 2 or more space dimensions.

Therefore, the following notation, in which the weights defining the scheme are given in a table corresponding to their spatial arrangement, is sometimes practical:

∆

, 1 h





0 1 0

1 −4 1

0 1 0



 .

Using this notation, the approximation

∆

, 1 6h





1 4 1

4 −20 4

1 4 1





of the Laplace operator is known as the 9-point stencil.

Using the same analysis based on Taylor series expansion reveals that

∆

u = ∆u + h

12 (u

+ 2u

+ u

) + O(h

).

As an approximation of ∆, the 9-point stencil ∆

therefore has the

same order of consistency as the 5-pont stencil ∆

.

However, it can still be used to construct a higher order approximation as follows: The leading error term is observed to contain the biharmonic operator applied to u:

u

+ 2u

+ u

= ∆(u

+ u

) = ∆(∆u) = ∆

u.

Since the exact solution u of the Poisson equation satisfies ∆u = f , we conclude that the local discretization error of the approximate equation

∆

U

= ˜ f

, 1 ≤ i ≤ n

, 1 ≤ j ≤ n

,

with the modified right-hand side

f ˜

:= f (x

, y

) + h

12 ∆f (x, y)

has order O(h

) as h → 0.

Example: We consider the Poisson equation ∆u = −1 on Ω = (−1, 1)

with homogeneous Dirichlet BCs (cf. pg. 181).

We discretize the problem using both ∆

and ∆

and compare the absolute value of the error in the point (0, 0) for different values of the mesh size h.

n × n ∆

∆

7 × 7 3.6e − 03 1.8e − 05 15 × 15 9.0e − 04 1.1e − 06 31 × 31 2.3e − 04 7.0e − 08 63 × 63 5.7e − 05 4.3e − 09 127 × 127 1.4e − 05 2.7e − 10 255 × 255 3.5e − 06 1.7e − 11

511 × 511 8.9e − 07 1.4e − 12

10^ï8 10^ï6 10^ï4 10^ï2

Number of grid points in each direction

Error in (0,0)

5ïpoint stencil 9ïpoint stencil

Note:

The last problem requires solving a linear system of dimension 261121 and

takes 5 s on a really old laptop.

Goals of Chapter 5

You should know how test examples for the Laplace and Poisson equation on simple domains can be constructed and be able to do this in some simple cases.

You should be able to discretize Dirichlet and Neumann BVPs for the Laplace and Poisson equation on simple domains.

You should be able to write Matlab programs for solving such problems.

You should know the 1D and 2D discretization matrices for the Laplacian with Dirichlet and Neumann BCs and their connection via Kronecker products.

You should know the order of convergence for the 5-point stencil.