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 Grid Generation 1 Discretisation 2 System of Linear Equations 3 Discretisation Stencil 4 Accuracy 5

Part II: Finite Element Methods FEM Main Ingredients 6 Weak Forms and Weak Solutions 7 Test and Ansatz Space 8 Discretisation 9 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: ∂ x 2 u ( x , y )+ ∂ 2 ∂ 2 in Ω = ( 0 , 1 ) 2 ∂ y 2 u ( x , y ) = f ( x , y ) Dirichlet boundary conditions: u ( x , y ) = g ( x , y ) on ∂ Ω

Part I Finite Differences

Grid Generation generate a grid on the given domain x i,j+1 h y x i−1,j x i,j x i+1,j x i,j−1 h z h y h x h x Compute values of unknown function u at each grid point: u ij ≈ u ( x ij ) u ijk ≈ u ( x ijk )

Finite Difference Discretisation Replace derivatives (at each grid point) by difference quotients: ∂ 2 u u ( x i + 1 , j ) − 2 u ( x i , j )+ u ( x i − 1 , j ) ∂ x 2 ( x i , j ) ≈ h 2 x ∂ 2 u u ( x i + 1 , j ) − 2 u ( x i , j )+ u ( x i − 1 , j ) ∂ x 2 ( x i , j ) ≈ h 2 x leads to linear system of equations: � 1 u i + 1 , j + u i , j + 1 − 4 u i , j h 2 � x i , j ∈ ( 0 , 1 ) 2 + u i , j − 1 + u i − 1 , j = f ( x i , j ) u ( x i , j ) = g ( x i , j ) x i , j ∈ ∂ Ω

System of Linear Equations write linear system as A h x h = f h A h is a sparse matrix (band structure):   B h − I ... ...   A h = 1 − I     ... ... h 2   − I   − I B h B h = tridiag ( − 1 , 4 , − 1 ) A h 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:   1 � � 1 D : 1 − 2 1 2 D : 1 − 4 1   1

Accuracy for 5-point Poisson stencil; order of accuracy: � u h − u � = O ( h 2 ) = O ( N − 2 ) curse of dimension: for that, we need O ( N d ) 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 compute a function as numerical solution; 1 search in a function space V h : u h = ∑ u j ϕ j ( x ) , span { ϕ 1 ,... ϕ J } = V h j solve weak form of PDE to reduce regularity 2 properties � � u ′′ = f v ′ u ′ dx = − → vf dx choose basis function with a local support , e.g.: 3 ϕ j ( x i ) = δ ij

Weak Forms and Weak Solutions consider a PDE Lu = f (e.g. Lu = ∆ u ) transformation to the weak form : � � � Lu , v � = vLu dx = 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 on δ Ω weak form: � � − Ω v ∆ u d Ω = Ω vf d Ω ∀ v apply Green’s formula (and boundary conditions): � � Ω ∇ v · ∇ u d Ω = Ω 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 u j ϕ j ( x ) dx = vf dx ∀ v ∈ V j where span { ϕ j } = W (“ansatz space”) choose a basis { ψ i } of the test space V ; then: � � � � ∑ ψ i L u j ϕ j ( x ) dx = ψ i f dx ∀ ψ i j leads to system of equations for unknowns u j usually V and W chosen identically (Ritz-Galerkin method)

Discretisation – Finite Elements L linear ⇒ system of linear equations � � � � ∑ ψ i L ϕ j ( x ) dx u j = ψ i f dx ∀ ψ i j � �� � =: A ij aim: make matrix A sparse → most A ij = 0 approach: local basis functions on a discretisation grid ψ j , ϕ j zero everywhere except in grid cells adjacent to grid point x j A ij = 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 − u FE � � a ≤ inf u h ∈ V � u − u h � a � h 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 � 1 0 ∇ v · ∇ u dx = 0 vf dx ∀ v compuational grid: x i = ih , (for i = 1 ,..., n − 1); mesh size h = 1 / n V = W : piecewise linear functions (on intervals [ x i , x i + 1 ] )

Nodal Basis  1 h ( x − x i − 1 ) x i − 1 < x < x i   1 ϕ i ( x ) := h ( x i + 1 − x ) x i < x < x i + 1   0 otherwise 1 0,8 0,6 0,4 0,2 0 0 0,2 0,4 0,6 0,8 1 x

Nodal Basis – System of Equations stiffness matrix:   2 − 1 ...   1 − 1 2     ... ... h   − 1   − 1 2 right hand sides (assume f ( x ) = α ∈ R ): � 1 � 1 0 ϕ i ( x ) f ( x ) dx = 0 ϕ i ( x ) α dx = α h system of equations very similar to finite differences

Hierarchical Basis 1 0,8 0,6 0,4 0,2 0 0 0,2 0,4 0,6 0,8 1 x 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 · ∇ u dx Ω ( i ) ∇ v · ∇ u dx i � � = ∑ Ω vf dx Ω ( i ) vf dx i leads to system of equations for each element: A ( i ) x = b ( i ) accumulate to obtain global system: x = ∑ ∑ A ( i ) b ( i ) i i � �� � =: A

Example: 1D Poisson Notation: notation: omit zero columns/rows (leaves only unknowns that are in Ω ( i ) ) 1D Poisson: Ω = [ 0 , 1 ] splitted into Ω ( i ) = [ x i − 1 , x i ] nodal basis; leads to element stiffness matrix: � � 1 − 1 A ( i ) = − 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 ) = [ x i − 1 , x i ] × [ x j − 1 , x j ] bilinear basis functions ϕ ij ( x , y ) = ϕ i ( x ) ϕ j ( y ) “pagoda” functions

Example: 2D Poisson (2) leads to element stiffness matrix:   − 1 − 1 2 − 1 2 2  − 1 − 1  2 − 1 A ( i ) =   2 2   − 1 − 1   − 1 2 2 2   − 1 − 1 − 1 2 2 2 accumulation leads to 9-point stencil   − 1 − 1 − 1   − 1 8 − 1   − 1 − 1 − 1

Typical workflow choose elements: 1 quadratic or cubic cells triangles (structured, unstructured) tetrahedra, etc. set up basis functions for each element Ω h ; 2 at all nodes x i ∈ Ω h ϕ i ( x i ) = 1 ϕ i ( x j ) = 0 for all j � = i for element stiffness matrix, compute all 3 � A ij = ϕ i L ϕ j d Ω Ω h accumulate global stiffness matrix 4