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

Scientific Computing I Michael Bader Outlines

Part I: Finite Differences Part II: Finite Element Methods

Scientific Computing I

Module 8: Discretisation of PDEs Michael Bader

Lehrstuhl Informatik V

Winter 2006/2007

Scientific Computing I Michael Bader Outlines

Part I: Finite Differences Part II: Finite Element Methods

Part I: Finite Differences

1

Grid Generation

2

Discretisation

3

System of Linear Equations

4

Discretisation Stencil

5

Accuracy

Scientific Computing I Michael Bader Outlines

Part I: Finite Differences Part II: Finite Element Methods

Part II: Finite Element Methods

6

FEM Main Ingredients

7

Weak Forms and Weak Solutions

8

Ansatz Functions

9

Weak Solutions

10 Test and Ansatz Space 11 Discretisation 12 A Road to Theory 13 Choosing Basis Functions

Nodal Basis

14 Element Stiffness Matrices

Example: 1D Poisson Example: 2D Poisson Workflow

Scientific Computing I Michael Bader Outlines

Part I: Finite Differences Part II: Finite Element Methods

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 ∂Ω

Scientific Computing I Michael Bader Grid Generation Discretisation System of Linear Equations Discretisation Stencil Accuracy

Part I Finite Differences

Scientific Computing I Michael Bader Grid Generation Discretisation System of Linear Equations Discretisation Stencil Accuracy

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)

Scientific Computing I Michael Bader Grid Generation Discretisation System of Linear Equations Discretisation Stencil Accuracy

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 ∈ ∂Ω

Scientific Computing I Michael Bader Grid Generation Discretisation System of Linear Equations Discretisation Stencil Accuracy

System of Linear Equations

write linear system as Ahxh = fh xh = (u1,1,...,u1,n,u2,1,...,u2,n,...,un,1,...,un,n) 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

Scientific Computing I Michael Bader Grid Generation Discretisation System of Linear Equations Discretisation Stencil Accuracy

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  

Scientific Computing I Michael Bader Grid Generation Discretisation System of Linear Equations Discretisation Stencil Accuracy

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

Scientific Computing I Michael Bader FEM Main Ingredients Weak Forms and Weak Solutions Ansatz Functions Weak Solutions Test and Ansatz Space Discretisation A Road to Theory Choosing Basis Functions

Nodal Basis

Element Stiffness Matrices

Example: 1D Poisson Example: 2D Poisson Workflow

Part II Finite Element Methods

Scientific Computing I Michael Bader FEM Main Ingredients Weak Forms and Weak Solutions Ansatz Functions Weak Solutions Test and Ansatz Space Discretisation A Road to Theory Choosing Basis Functions

Nodal Basis

Element Stiffness Matrices

Example: 1D Poisson Example: 2D Poisson Workflow

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} = Wh

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

Scientific Computing I Michael Bader FEM Main Ingredients Weak Forms and Weak Solutions Ansatz Functions Weak Solutions Test and Ansatz Space Discretisation A Road to Theory Choosing Basis Functions

Nodal Basis

Element Stiffness Matrices

Example: 1D Poisson Example: 2D Poisson Workflow

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} = Vh the “best” solution uh in this function space is wanted

Scientific Computing I Michael Bader FEM Main Ingredients Weak Forms and Weak Solutions Ansatz Functions Weak Solutions Test and Ansatz Space Discretisation A Road to Theory Choosing Basis Functions

Nodal Basis

Element Stiffness Matrices

Example: 1D Poisson Example: 2D Poisson Workflow

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

Scientific Computing I Michael Bader FEM Main Ingredients Weak Forms and Weak Solutions Ansatz Functions Weak Solutions Test and Ansatz Space Discretisation A Road to Theory Choosing Basis Functions

Nodal Basis

Element Stiffness Matrices

Example: 1D Poisson Example: 2D Poisson Workflow

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

Scientific Computing I Michael Bader FEM Main Ingredients Weak Forms and Weak Solutions Ansatz Functions Weak Solutions Test and Ansatz Space Discretisation A Road to Theory Choosing Basis Functions

Nodal Basis

Element Stiffness Matrices

Example: 1D Poisson Example: 2D Poisson Workflow

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

Scientific Computing I Michael Bader FEM Main Ingredients Weak Forms and Weak Solutions Ansatz Functions Weak Solutions Test and Ansatz Space Discretisation A Road to Theory Choosing Basis Functions

Nodal Basis

Element Stiffness Matrices

Example: 1D Poisson Example: 2D Poisson Workflow

Choose Test and Ansatz Space

search for weak solutions uh in a certain function space Wh

• vL
• ∑

j

ujϕj(x)

• dx =
• vf dx

∀v ∈ V where span{ϕj} = Wh (“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 V is often chosen to be identical to Wh (Ritz-Galerkin method)

Scientific Computing I Michael Bader FEM Main Ingredients Weak Forms and Weak Solutions Ansatz Functions Weak Solutions Test and Ansatz Space Discretisation A Road to Theory Choosing Basis Functions

Nodal Basis

Element Stiffness Matrices

Example: 1D Poisson Example: 2D Poisson Workflow

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

Scientific Computing I Michael Bader FEM Main Ingredients Weak Forms and Weak Solutions Ansatz Functions Weak Solutions Test and Ansatz Space Discretisation A Road to Theory Choosing Basis Functions

Nodal Basis

Element Stiffness Matrices

Example: 1D Poisson Example: 2D Poisson Workflow

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)

Scientific Computing I Michael Bader FEM Main Ingredients Weak Forms and Weak Solutions Ansatz Functions Weak Solutions Test and Ansatz Space Discretisation A Road to Theory Choosing Basis Functions

Nodal Basis

Element Stiffness Matrices

Example: 1D Poisson Example: 2D Poisson Workflow

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])

Scientific Computing I Michael Bader FEM Main Ingredients Weak Forms and Weak Solutions Ansatz Functions Weak Solutions Test and Ansatz Space Discretisation A Road to Theory Choosing Basis Functions

Nodal Basis

Element Stiffness Matrices

Example: 1D Poisson Example: 2D Poisson Workflow

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

Scientific Computing I Michael Bader FEM Main Ingredients Weak Forms and Weak Solutions Ansatz Functions Weak Solutions Test and Ansatz Space Discretisation A Road to Theory Choosing Basis Functions

Nodal Basis

Element Stiffness Matrices

Example: 1D Poisson Example: 2D Poisson Workflow

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

Scientific Computing I Michael Bader FEM Main Ingredients Weak Forms and Weak Solutions Ansatz Functions Weak Solutions Test and Ansatz Space Discretisation A Road to Theory Choosing Basis Functions

Nodal Basis

Element Stiffness Matrices

Example: 1D Poisson Example: 2D Poisson Workflow

Element Stiffness Matrices

domain Ω splitted into finite elements Ω(k): Ω = Ω(1) ∪Ω(2) ∪···∪Ω(n) 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

Scientific Computing I Michael Bader FEM Main Ingredients Weak Forms and Weak Solutions Ansatz Functions Weak Solutions Test and Ansatz Space Discretisation A Road to Theory Choosing Basis Functions

Nodal Basis

Element Stiffness Matrices

Example: 1D Poisson Example: 2D Poisson Workflow

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)

Scientific Computing I Michael Bader FEM Main Ingredients Weak Forms and Weak Solutions Ansatz Functions Weak Solutions Test and Ansatz Space Discretisation A Road to Theory Choosing Basis Functions

Nodal Basis

Element Stiffness Matrices

Example: 1D Poisson Example: 2D Poisson Workflow

Example: 1D Poisson

Notation:

• mit zero columns/rows

(leaves only unknowns that are in Ω(k)) 1D Poisson: Ω = [0,1] splitted into Ω(k) = [xk−1,xk] nodal basis; leads to element stiffness matrix: A(k) =

• 1

−1 −1 1

• stencil notation: [1∗ −1], [−1

1∗ ] [−1 1∗ ] + [1∗ −1] → [−1 2 −1]

Scientific Computing I Michael Bader FEM Main Ingredients Weak Forms and Weak Solutions Ansatz Functions Weak Solutions Test and Ansatz Space Discretisation A Road to Theory Choosing Basis Functions

Nodal Basis

Element Stiffness Matrices

Example: 1D Poisson Example: 2D Poisson Workflow

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

Scientific Computing I Michael Bader FEM Main Ingredients Weak Forms and Weak Solutions Ansatz Functions Weak Solutions Test and Ansatz Space Discretisation A Road to Theory Choosing Basis Functions

Nodal Basis

Element Stiffness Matrices

Example: 1D Poisson Example: 2D Poisson Workflow

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   

Scientific Computing I Michael Bader FEM Main Ingredients Weak Forms and Weak Solutions Ansatz Functions Weak Solutions Test and Ansatz Space Discretisation A Road to Theory Choosing Basis Functions

Nodal Basis

Element Stiffness Matrices

Example: 1D Poisson Example: 2D Poisson Workflow

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