The Model Problem Scientific Computing I 2D Poisson Equation on - - PDF document

SLIDE 1

Scientific Computing I

Module 8: Discretisation of PDEs Michael Bader

Lehrstuhl Informatik V

Winter 2007/2008

The Model Problem

2D Poisson Equation on unit square: ∂ 2 ∂x2u(x,y)+ ∂ 2 ∂y2u(x,y) = f(x,y) in Ω = (0,1)2 Dirichlet boundary conditions: u(x,y) = g(x,y)

• n ∂Ω

Part I Finite Differences Grid Generation

generate a grid on the given domain

xi,j xi−1,j xi+1,j xi,j+1 xi,j−1 hx hy hx hz hy

Compute values of unknown function u at each grid point: uij ≈ u(xij) uijk ≈ u(xijk)

Finite Difference Discretisation

Replace derivatives (at each grid point) by difference quotients: ∂ 2u ∂x2 (xi,j) ≈ u(xi+1,j)−2u(xi,j)+u(xi−1,j) h2

x

∂ 2u ∂y2 (xi,j) ≈ u(xi,j+1)−2u(xi,j)+u(xi,j−1) h2

y

leads to linear system of equations (h := hx = hy):

1 h2

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

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

• =

f(xi,j) xi,j ∈ (0,1)2 u(xi,j) = g(xi,j) xi,j ∈ ∂Ω

System of Linear Equations

• bjective: write linear system in matrix-vector-form:

Ahxh = fh xh contains the unknowns uij ⇒ requires sequentialisation of the unknowns for example, with simple row-wise numbering of the grid points: xh = (u1,1,...,u1,n,u2,1,...,u2,n,...,un,1,...,un,n)

SLIDE 2

System of Linear Equations (2)

Ah is a sparse matrix (five-diagonal) Ah is a block-tridiagonal matrix: Ah = 1 h2       Bh I I ... ... ... ... I I Bh       Bh = tridiag(1,−4,1) I the identity matrix boundary values to right-hand side

Stencil Notation

illustrate matrix structure as a discretisation stencil represents one line of the matrix matrix elements placed according to their corresponding geometrical position stencils for the Poisson equation (h2 factors ignored): 1D :

• 1

−2 1

• 2D :

  1 1 −4 1 1  

Accuracy

for 5-point Poisson stencil;

• rder of accuracy: uh −u = O(h2) = O(N−2)

curse of dimension: for that, we need O(Nd) points Possibilities of an improvement: use higher-order discretisation

via higher order terms of Taylor series leads to larger stencils (involving more neighbouring grid points)

use locally refined (adaptive) grids

Part II Finite Element Methods Finite Elements – Main Ingredients

1

compute a function as numerical solution; search in a function space Wh: uh = ∑

j

ujϕj(x), span{ϕ1,...,ϕJ} = Wh

2

solve weak form of PDE to reduce regularity properties u′′ = f − →

• v′u′ dx =
• vf dx

3

choose basis functions with local support, e.g.: ϕj(xi) = δij

Choose Test and Ansatz Space

search for solution functions uh of the form uh = ∑

j

ujϕj(x) the basis (ansatz) functions ϕj(x) build a vector space (or function space) Wh span{ϕ1,...,ϕJ} = Wh the “best” solution uh in this function space is wanted

SLIDE 3

Example: Nodal Basis

ϕi(x) :=     

1 h(x −xi−1)

xi−1 < x < xi

1 h(xi+1 −x)

xi < x < xi+1

• therwise

0,6 0,4 0,8 0,2 1 1 0,4 0,2 x 0,8 0,6

Weak Forms and Weak Solutions

consider a PDE Lu = f (e.g. Lu = ∆u) transformation to the weak form: Lu,v =

• vLudx =
• vf dx = f,v

∀v ∈ V V a certain class of functions “real solution” u also solves the weak form Lu,v a bilinear form; often written as: a(u,v) = f,v ∀v ∈ V

Weak Form of the Poisson Equation

Poisson equation with Dirichlet conditions: −∆u = f in Ω,u = 0

• n δΩ

weak form: −

• Ω v∆udΩ =
• Ω vf dΩ

∀v apply Green’s formula: −

• Ω v∆udΩ =
• Ω ∇v ·∇udΩ−
• ∂Ω v ·∇uds

choose functions v such that v = 0 on ∂Ω:

• Ω ∇v ·∇udΩ =
• Ω vf dΩ

∀v

Weak Form of the Poisson Equation (2)

Poisson equation with Dirichlet conditions: −∆u = f in Ω,u = 0

• n δΩ

transformed into weak form:

• Ω ∇v ·∇udΩ =
• Ω vf dΩ

∀v weaker requirements for a solution u: twice differentiabale → first derivative integrable remember use of nodal basis: availability of first vs. second derivative!

Choose Test and Ansatz Space

search for solutions uh in a function space Wh: uh = ∑

j

ujϕj(x) where span{ϕj} = Wh (“ansatz space”) insert into weak solution

• vL
• ∑

j

ujϕj(x)

• dx =
• vf dx

∀v ∈ V

Choose Test and Ansatz Space (2)

choose a basis {ψi} of the test space V then: if all basis functions ψi satisfy

• ψiL
• ∑

j

ujϕj(x)

• dx =
• ψif dx

∀ψi then all v ∈ V satisfy the equation leads to system of equations for unknowns uj (one equation per test basis function ψi) V is often chosen to be identical to Wh (Ritz-Galerkin method)

SLIDE 4

Discretisation – Finite Elements

L linear ⇒ system of linear equations

∑

j

• ψiLϕj(x)dx
• =:Aij
• uj =
• ψif dx

∀ψi aim: make matrix A sparse → most Aij = 0 approach: local basis functions on a discretisation grid ψj,ϕj zero everywhere except in grid cells adjacent to grid point xj Aij = 0, if ψi and ϕj don’t overlap

Example Problem: Poisson 1D

in 1D: u′′(x) = f(x) on Ω = (0,1),

• hom. Dirichlet boundary cond.: u(0) = u(1) = 0

weak form:

1

0 v′(x)·u′(x)dx =

1

0 v(x)f(x)dx

∀v compuational grid: xi = ih, (for i = 1,...,n−1); mesh size h = 1/n V = W: piecewise linear functions (on intervals [xi,xi+1])

Nodal Basis

ϕi(x) :=     

1 h(x −xi−1)

xi−1 < x < xi

1 h(xi+1 −x)

xi < x < xi+1

• therwise

0,6 0,4 0,8 0,2 1 1 0,4 0,2 x 0,8 0,6

Nodal Basis – System of Equations

stiffness matrix: 1 h       2 −1 −1 2 ... ... ... −1 −1 2       right hand sides (assume f(x) = α ∈ R):

1

0 ϕi(x)f(x)dx =

1

0 ϕi(x)α dx = αh

system of equations very similar to finite differences

Element Stiffness Matrices

domain Ω splitted into finite elements Ω(k): Ω = Ω(1) ∪Ω(2) ∪···∪Ω(n)

• bservation: basis functions are defined

element-wise use:

b

a f(x)dx =

c

a f(x)dx+

b

c f(x)dx

element-wise evaluation of the integrals:

• Ω ∇v ·∇udx

= ∑

k

• Ω(k) ∇v ·∇udx
• Ω vf dx

= ∑

i

• Ω(i) vf dx

Element Stiffness Matrices (2)

leads to local stiffness matrices for each element:

• Ω(k) ∇φi ·∇φj dx
• =:A(k)

ij

and respective element systems: A(k)x = b(k) accumulate to obtain global system:

∑

k

A(k)

=:A

x = ∑

k

b(k)

SLIDE 5

Element Stiffness Matrices (3)

Some comments on notation: assume: 1D problem, n elements (i.e. intervals) in each element only two basis functions are non-zero! hence, almost all A(k)

ij

are zero: A(k)

ij

=

• Ω(k) ∇φi ·∇φj dx
• nly 2×2 elements of A(k) are non-zero

therefore convention to omit zero columns/rows ⇒ leaves only unknowns that are in Ω(k)

Example: 1D Poisson

Ω = [0,1] splitted into Ω(k) = [xk−1,xk] nodal basis; leads to element stiffness matrix: A(k) =

• 1

−1 −1 1

• consider only two elements:

A(1) +A(2) =   1 −1 −1 1  +   1 −1 −1 1   in stencil notation: [−1 1∗ ] + [1∗ −1] → [−1 2 −1]

Example: 2D Poisson

−∆u = f on domain Ω = [0,1]2 splitted into Ω(i,j) = [xi−1,xi]×[xj−1,xj] bilinear basis functions ϕij(x,y) = ϕi(x)ϕj(y) “pagoda” functions

Example: 2D Poisson (2)

leads to element stiffness matrix: A(k) =       2 −1

2

−1

2

−1 −1

2

2 −1 −1

2

−1

2

−1 2 −1

2

−1 −1

2

−1

2

2       accumulation leads to 9-point stencil    −1 −1 −1 −1 8 −1 −1 −1 −1   

Typical workflow

1

choose elements:

quadratic or cubic cells triangles (structured, unstructured) tetrahedra, etc.

2

set up basis functions for each element Ω(k); for example, at all nodes xi ∈ Ω(k) ϕi(xi) = 1 ϕi(xj) = for all j = i

3

for element stiffness matrix, compute all A(k)

ij

=

• Ω(k) ϕiLϕj dΩ

4

accumulate global stiffness matrix

Time-Dependent Problems

Example: 1D Heat Equation ut = uxx +f on domain Ω = [0,1] for t ∈ [0,tend] spatial discretisation: weak form

• vut dx

=

• vuxx dx+
• vf dx

∂ ∂t

• vudx
• =
• vuxx dx+
• vf dx

spatial discretisation – finite elements:

∂ ∂t (Mhuh) = Ahuh +fh

Mh: mass matrix, Ah: stiffness matrix, uh = uh(t)

SLIDE 6

Time-Dependent Problems (2)

Solve a system of ordinary differential equations: after spatial discretisation (Mh constant):

∂ ∂tMh (uh) = Ahuh +fh

uh a vector of time-dependent functions: uh = (u1(t),...,ui(t),...,un(t))T usually: approximate Mh by a simpler matrix (diagonal matrix, e.g.) → “mass lumping”

A Road to Theory

weak formulation is equivalent to variational approach: solution u minimises an energy functional best approximation property:

• u−uFE

h

• a ≤ inf

uh∈V u−uha

in terms of the norm induced by the bilinear form a (energy norm) thus: error bounded by interpolation error (in energy norm)