Scientific Computing I Module 8: Discretisation of PDEs Michael - - PowerPoint PPT Presentation

Scientific Computing I

Module 8: Discretisation of PDEs Michael Bader

Lehrstuhl Informatik V

Winter 2005/2006

Part I: Finite Differences

1

Grid Generation

2

Discretisation

3

System of Linear Equations

4

Discretisation Stencil

5

Accuracy

Part II: Finite Element Methods

6

FEM Main Ingredients

7

Weak Forms and Weak Solutions

8

Test and Ansatz Space

9

Discretisation

10 A Road to Theory 11 Choosing Basis Functions

Nodal Basis Hierarchical Basis

12 Element Stiffness Matrices

Example: 1D Poisson Example: 2D Poisson Workflow

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 ∂x2 (xi,j) ≈ u(xi+1,j)−2u(xi,j)+u(xi−1,j) h2

x

leads to linear system of equations:

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

write linear system as Ahxh = fh Ah is a sparse matrix (band structure): Ah = 1 h2       Bh −I −I ... ... ... ... −I −I Bh       Bh = tridiag(−1,4,−1) Ah block-tridiagonal (5-diagonal) matrix boundary values to right-hand side

Stencil Notation

illustrate matrix structure as a discretisation stencil represents one line of the matrix elements placed according to the geometrical position stencils for the Poisson equation: 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 Vh: uh = ∑

j

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

2

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

• v′u′ dx =
• vf dx

3

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

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(Lu,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 (and boundary conditions):

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

∀v weaker requirements for a solution u: twice differentiabale → first derivative integrable

Choose Test and Ansatz Space

search for weak solutions u in a certain function space W

• vL
• ∑

j

ujϕj(x)

• dx =
• vf dx

∀v ∈ V where span{ϕj} = W (“ansatz space”) choose a basis {ψi} of the test space V; then:

• ψiL
• ∑

j

ujϕj(x)

• dx =
• ψif dx

∀ψi leads to system of equations for unknowns uj usually V and W chosen identically (Ritz-Galerkin method)

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

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 interpoplation error (in energy norm)

Example Problem: Poisson 1D

1D Poisson’s equation on Ω = [0,1], homogeneous Dirichlet boundary conditions weak form:

1

0 ∇v ·∇udx =

1

0 vf 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

Hierarchical Basis

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

leads to diagonal stiffness matrix! (for 1D Poisson) solution function identical to that with nodal basis (same function space)

Element Stiffness Matrices

domain Ω splitted into finite elements Ω(i) element-wise evaluation of the integrals:

• Ω ∇v ·∇udx

= ∑

i

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

= ∑

i

• Ω(i) vf dx

leads to system of equations for each element: A(i)x = b(i) accumulate to obtain global system:

∑

i

A(i)

=:A

x = ∑

i

b(i)

Example: 1D Poisson

Notation: notation: omit zero columns/rows (leaves only unknowns that are in Ω(i)) 1D Poisson: Ω = [0,1] splitted into Ω(i) = [xi−1,xi] nodal basis; leads to element stiffness matrix: A(i) =

• 1

−1 −1 1

• 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(i) =       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 Ωh; at all nodes xi ∈ Ωh ϕi(xi) = 1 ϕi(xj) = for all j = i

3

for element stiffness matrix, compute all Aij =

• Ωh

ϕiLϕj dΩ

4

accumulate global stiffness matrix