Electrostatic PDEs via Finite Elements* Industrial Strength = Finite - PowerPoint PPT Presentation

Electrostatic PDEs via Finite Elements* Industrial Strength � = Finite Differences Rubin H Landau Sally Haerer, Producer-Director Based on A Survey of Computational Physics by Landau, Páez, & Bordeianu with Support from the National Science Foundation Course: Computational Physics II 1 / 1

Finite-Element Method (FEM) Divide Space into Elements; Approximate Solution on Each x Match element edges Active development, � i, j-1 y ↑ convergence, robust, i+1, j i-1, j i, j precision i, j+1 Flexible, variable domains FEM: ↓ computer time Finite Difference: how FEM: ↑ analysis approximate derivatives Not just “on” grid FEM: solution on elements 2 / 1

Sample 1-D Basis & 2-D Elements (Wikipedia) 1-D elements, 1-D Elemental Solutions x N element b �� � �� � x x x 0 node 2-D 3 / 1

Electric Field V ( x ) from Charge Density = ? 2 Conducting Plates x N element b �� � �� � x x x 0 node d 2 U ( x ) Solve Poisson Equation = − 4 πρ ( x ) = − 1 dx 2 BC: U ( x = a = 0 ) = 0, U ( x = b = 1 ) = 1 √ Analytic solution: U ( x ) = x 2 ( 3 − x ) 4 / 1

Finite Element Method Analytic + Numeric x N element b �� � �� � x x x 0 node Split domain into subdomains = elements Lines = subdomains, Dots = nodes x i Hypothesize trial solution in subdomain Global fit trial solution to exact (matching subdomains) � Via PDE weak form ≡ min (approximate - exact) Implement BC, Solve linear equations 5 / 1

Approximate Subdomain Solutions φ i ( x ) Solution on Element = Basis � � � . . . 0 1 N x 0 x 1 . . . x N-1 Left: Piecewise-linear Overlapping elemental solutions φ i ( x ) Right: -quadratic Essentially basis functions Match φ i s together Elements cover space 1-D elements: lines φ i ( x ) within 1 element 2-D: rectangles, triangles 6 / 1

Weak Form of PDE (Global Best Fit) − d 2 U ( x ) U ( x ) = Solution: = 4 πρ ( x ) dx 2 Best possible global agreement with exact Assume basis Φ( a ) = Φ( b ) = 0, adjust for BC later Convert to integral/weak form, ODE × Φ , integrate by parts − d 2 U ( x ) Φ( x ) = 4 πρ ( x )Φ( x ) (1) dx 2 � b � b − dU ( x ) dx dU ( x ) Φ( x ) | b a + Φ ′ ( x ) = dx 4 πρ ( x ) Φ( x ) (2) dx dx a a � b � b dx dU ( x ) Φ ′ ( x ) = ⇒ dx 4 πρ ( x ) Φ( x ) (3) dx a a 7 / 1

Spectral Decomposition of Solution U ( x ) (Galerkin) Hat Basis: Linear Interpolation of U ( nodes ) � � � . . . N − 1 0 1 N � U ( x ) ≃ α j φ j ( x ) (1) j = 0 x 0 x 1 . . . x N-1 φ i : elemental solution =?? Simple: Piecewise-cont α j = ? φ i ( A , B ) = 0, otherwise x − x i − 1  for x i − 1 ≤ x ≤ x i , h i − 1 ,  φ i ( x ) = ( h i = x i + 1 − x i ) (2) x i + 1 − x for x i ≤ x ≤ x i + 1 , ,  h i φ i ( x i ) = 1 ⇒ α i = U ( x i ) (nodes) N − 1 � U ( x ) ≃ U ( x j ) φ j ( x ) (3) j = 0 8 / 1

Determine α j via Linear Equations � α j φ j ( x ) = � α j U ( x node ) U ( x ) ≃ Sub U ( x ) , Φ( x ) = φ i into integral Poisson Eq: � b N − 1 � b α j d φ j ( x ) d φ i � = dx 4 πρ ( x ) φ i ( x ) , i = 0 , N − 1 dx (1) dx dx a a j = 0 ⇒ N linear equations for α j ’s � b � b � b φ ′ i φ ′ 0 dx + α 1 φ ′ i φ ′ 1 dx + · · · + α N − 1 φ ′ i φ ′ N − 1 dx α 0 a a a � b = 4 πρφ i dx , i = 0 , N − 1 (2) a � φ i = known hat functions ⇒ φ ′ i φ ′ N − 1 dx = h i = known 9 / 1

As Matrix Equations for α j Stn Form: Stiffness Matrix A , Unknown Load Vector y A y = b (1) � x 1 x 0 dx 4 πρ ( x ) φ 0 ( x )    α 0  � x 2 x 1 dx 4 πρ ( x ) φ 1 ( x ) α 1       y = b = (2)  ,  ...  ...         � x N  α N − 1 x N − 1 dx 4 πρ ( x ) φ N − 1 ( x ) 1 1 − 1 − 1  h 0 + 0 . . .  h 1 h 1 h 0 − 1 1 1 − 1 h 1 + 0 . . .   A = h 1 h 2 h 2   ... ...   1 1 1 1 − − h N − 2 + h N − 1 h N − 2 h N − 1 A : known, no ∆ with ρ , b : known analytic or quadrature 10 / 1

Impose Dirichlet (Function) Boundary Conditions φ ends = 0 ⇒ U ends = 0 φ N ( x ) satisfies homogeneous eq (Laplace) ∇ 2 φ N = 0 Add particular solution to satisfy BC @ A (B later) N − 1 � U ( x ) = α j φ j ( x ) + U a φ N ( x ) j = 0 Sub U ( x ) − U a φ N ( x ) into weak form (homo BC): A 0 , 0 · · · A 0 , N − 1 0 b 0 − A 0 , 0 U a     ... ...   b ′ =   A =  ,         A N − 1 , 0 · · · A N − 1 , N − 1 0 b N − 1 − A N − 1 , 0 U a    0 0 · · · 1 U a b ′ b ′ i = b i − A i , 0 U a , i = 1 , . . . , N − 1 , N = U a 11 / 1

Implementation LaplaceFEM.py • 3 elements (dots) poor Solve Ay = b ′ 1-D: ∼ 100–1000 eqs 1 N = 3 (scaled) 3-D: ∼ millions N = 11 U Number of calcs ∼ N 2 0 ⇒ keep N small 0 1 x • N = 11 excellent 12 / 1

FEM Exercise Examine solution for a = 0 , 1 b = 1 , U a = 0 , U b = 1 Compare piecewise-quadratic basis 2 Solve 3 ≤ N ≤ 1000 elements 3 Check stiffness matrix A triangular 4 Verify accuracy of integrations in load vector b 5 Verify solution of Ay = b 6 Plot U ( x ) ; N = 10, 100, 1000; compare analytic 7 � b � � 1 U FEM ( x ) − U exact ( x ) � � # Rel error = log 10 dx � � b − a U exact ( x ) � � a Plot error versus x , N = 10, 100, 1000 8 13 / 1