Baixe o app para aproveitar ainda mais
Prévia do material em texto
Numerical Methods for Elliptic PDE 1 Numerical Methods for Solving Equilibrium Problem (Elliptic PDE) Numerical Methods for Elliptic PDE 2 Consider 2-D Poisson equation of the form 2 2 2 2 f f S x y ∂ ∂ + = ∂ ∂ , 0 , 0x a y b≤ ≤ ≤ ≤ (1) xi-1/2 xi+1/2 yj+1/2 yj-1/2 Numerical Methods for Elliptic PDE 3 Using the 2nd-order central difference scheme, Eq. (1) can be discretized over a uniform mesh as , -1 , , 1 -1, , 1, ,2 2 ( 2 ) ( 2 )i j i j i j i j i j i j i j f f f f f f S y x + +− + − ++ = ∆ ∆ (2a) In compact form, 2 2 2 21, , 1 , 1, , 1 ,2(1 )i j i j i j i j i j i jf f f f f S xβ β β− − + ++ − + + + = ∆ (2b) i = 2, …, (I-1), j = 2, …., (J-1) where /x yβ = ∆ ∆ with / ( 1)x a I∆ = − and /( 1)y b J∆ = − . Numerical Methods for Elliptic PDE 4 Alternatively, we have , 1, , , 1 , , , 1, , , 1 , W S P E N i j i j i j i j i j i j i j i j i j i j i ja f a f a f a f a f b− − + ++ + + + = (2c) where 2 , ,i j i jb S x= ∆ , , 1 W i ja = , 2 , E i ja β= , , 1 S i ja = , 2 , N i ja β= , and , , , , ,( ) P W E S N i j i j i j i j i ja a a a a= − + + + . Numerical Methods for Elliptic PDE 5 Methods of Solving System of Linear Algebraic Equation Systems of linear algebraic equations can be expressed as 1 N ij j i j a f b = =∑ ( 1, , )i N= or [ ]{ } { }A f B= (1) where [ ]A is the coefficient matrix { }B is the non-homogeneous vector 11 12 1 21 22 2 1 2 [ ] N N N N NN a a a a a a A a a a = , 1 2{ } N f f f f = , and 1 2{ } N b b B b = Numerical Methods for Elliptic PDE 6 There are two basic approaches for solving system of linear algebraic equations: (a) Direct methods: Cramer’s rule (matrix inversion) Gauss elimination LU decomposition (b) Iterative methods: Jacobi iteration Gauss-Seidel iteration Relaxation Conjugate Gradient Methods user 註解 沒有用 Numerical Methods for Elliptic PDE 7 Direct Methods [LU Decomposition] The coefficient matrix [A] is decomposed or factorized into [ ] [ ][ ]A L U= (1) where 11 21 22 1 2 0 0 0 [ ] N N NN l l l L l l l = and 11 12 1 22 20[ ] 0 0 N N NN u u u u u U u = (2) (Lower triangular matrix) (Upper triangular matrix) Numerical Methods for Elliptic PDE 8 Now [ ]{ } { }[ ][ ]{ }A f L U f B= = and { }-1 -1 [ ] { } [ ] [ ][ ]{ } [ ] I B L L U f L B ′ = ⇒ [ ]{ } { }U f B′= ⇒ [ ]{ } { } {B } [ ][ ]{ }A f L U f B ′ = = ⇒ [ ]{ } { }L B B′ = It follows that [ ]{ } { } (Forward substitution) [ ]{ } { } (Backward substitution) L B B U f B ′ = ′= [LU Decomposition Algorithm] (a) Forward substitution: using [ ]{ } { }L B B′ = to solve for { }B′ . (b) Backward substitution: solve for { }f by [ ]{ } { }U f B′= . Numerical Methods for Elliptic PDE 9 Tri-Diagonal Matrix Algorithm (TDMA) (Thomas Algorithm) Consider a system of algebraic equations with tri-diagonal coefficient matrix 1 1 1 1 2 2 2 2 2 -1 -1 -1 -1 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 i i i i i N N N N N N N N N b c f s a b c f s a b c f s a b c f s a b f s = or Numerical Methods for Elliptic PDE 10 1 1 1 2 1 -1 1 -1 (1) , ( 2, 3, , 1) (2) (3) i i i i i i i N N N N N b f c f s a f b f c f s i N a f b f s + + = + + = = … − + = We postulate the existence of two vectors { }D and { }E such that 1i i i if d f e+= + (4) ⇒ -1 -1 -1i i i if d f e= + (5) Numerical Methods for Elliptic PDE 11 Substituting Eq. (5) into Eq. (2) and solving for fi gives -1 -1 1( )i i i i i i i i ia d f e b f c f s++ + + = ⇒ -1 1 -1( ) ( )i i i i i i i i ib a d f c f s a e++ = − + − ⇒ -11 -1 -1 ( ) ( ) ( ) i i i i i i i i i i i i i i c s a ef f b a d b a d d e + − − = + + + (6) From Eqs. (4) and (6) we have -1 -1 -1 ( ) ( ) ( ) i i i i i i i i i i i i cd b a d s a ee b a d − = + − = + 2 i N≤ ≤ (7) Numerical Methods for Elliptic PDE 12 Furthermore, to start the forward substitution of Eq. (7) the coefficients 1d and 1e can be evaluated from Eq. (1) as 1 1 1 2 1b f c f s+ = ⇒ 1 11 2 1 2 1 1 1 ( ) ( )c sf f d f e b b − = + = + ⇒ 1 1 1 1 1 1 ( ) ( ) cd b se b = − = (8) Similarly, from Eq. (3) -1N N N N Na f b f s+ = ⇒ -1 ( ) ( )N NN N N N b sf f a a = − + (9) Numerical Methods for Elliptic PDE 13 Further application of Eq. (5) for i = (N-1) and comparison with Eq. (9) gives -1 -1 -1 ( ) ( )N NN N N N N N N b sf d f e f a a = + = − + ⇒ -1 -1[ ( )] [( ) ]N NN N N N N b sd f e a a + = − ⇒ -1 -1 ( ) ( ) N N N N N N N s a ef b a d − = + (16) Next, the backward substitution can be implemented using Eqs. (5) and (16) as 1i i i if d f e+= + ( 1), ( 2), , 2, 1i N N= − − … (17) Numerical Methods for Elliptic PDE 14 Furthermore, let -1 -1 i i i i i i i i b b a d s s a e ′ = + ′ = − (18) Eq. (7) becomes then i i i i i i cd b se b = − ′ ′ = ′ 2 i N≤ ≤ (19) Numerical Methods for Elliptic PDE 15 Since 1 1 1 1 1 1 1 1 1 1 c cd b b s se b b = − = − ′ ′ = = ′ ⇒ 1 1 1 b b s s ′ = ′ = (20) Eq. (17) can be further expressed as 1( )i i ii i s c ff b +′ −= ′ ( 1), ( 2), , 2, 1i N N= − − … (22) where -1 -1 -1 /i i i i i i i i i i i r a b b b rc s s r s ′= ′ = − ′ ′= − (23) Numerical Methods for Elliptic PDE 16 Likewise, Eq. (16) becomes -1 -1 ( ) ( ) N N N N N N N N N s a e sf b a d b ′− = = ′+ (24) Numerical Methods for Elliptic PDE 17 [TDMA Algorithm] (a) Initialize the variables: 1 1b b′ = and 1s s′ = (b) Recurrently calculate the following equations in increasing order of 2,3, , i N= … : -1 -1 -1 / (Forward substitution) i i i i i i i i i i i r a b b b rc s s r s ′= ′ = − ′ ′= − (c) Calculate /N N Nf s b′ ′= (d) Calculate the following equations in decreasing order of i: 1( )i i ii i s c ff b +′ −= ′ ( 1), ( 2), , 2, 1i N N= − − … . (Backward substitution) Numerical Methods for Elliptic PDE 18 FORTRAN Subroutine of TDMA SUBROUTINE TRIDG(IL, IU, A, B, C, S) DIMENSION A(1), B(1), C(1), S(1) C….Subroutine to perform solution algorithm for the tri-diagonal C…..system of algebraicequations C….. C…..IL = subscript of first equation C…..IU = subscript of the last equation C…..A = coefficient to the left of the main diagonal C…..B = coefficient of the main diagonal C…..C = coefficient to right of the main diagonal C…..S = Element in the constant vector C…… C…..Forward substitution to establish upper tri-diagonal matrix C….. DO 10 I = LL+1, IU R = A(I)/B(I−1) B(I)= B(I)−R*C(I−1) S(I) = S(I)−R*S(I−1) 10 CONTINUE C….. C…..Backward substitution C…. S(IU) = S(IU)/B(IU) DO 20 I = IU−1, 1, −1 S(I) = (S(I)−C(I)*S(I+1))/B(I) 20 CONTINUE RETURN END Numerical Methods for Elliptic PDE 19 Iterative Methods [Basic Concepts] Consider the system of algebraic equations of [ ]{ } { }A f B= (1) Iteration Algorithm: (a) Assume an initial solution vector 0{ }f . (b) Generate an improved solution vector { }mf , m = 1, 2, …, based on some strategy for reducing the difference between the iterating solution vector { }mf and the true solution vector { }f , i.e., Numerical Methods for Elliptic PDE 20 After m iterations we have an appropriate solution vector, { }mf , there is a non-zero residual vector { }mR as [ ]{ } { } { } {0}m mA f B R− = ≠ (2) Substrating Eq. (2) from Eq. (1) gives [ ]({ } { }) { }m mA f f R− = (3) Let the convergence error vector { }mε defined as { } { } { }m mf fε = − (4) Eq. (3) becomes then [ ]{ } { }m mA Rε = (5) Numerical Methods for Elliptic PDE 21 (c) Repeat step (b) till convergence, i.e., { } {0}mR → ⇒ { } {0}mε → ⇒ { } { }mf f→ user 螢光標示 Numerical Methods for Elliptic PDE 22 Generalized Iterative Scheme The system of algebraic equations can be modified by a non-singular pre-conditioning matrix [ ]P as [ ]{ } { }[ ]( )P A f B= or { }* *[ ] { }A f B= (6) where *[ ] [ ][ ]A P A= and *{ } [ ]{ }B P B= . Decomposing the coefficient matrix to *[ ] [ ] [ ]A M N= − (7) { } *([ ] [ ]) { }M N f B− = (8) Numerical Methods for Elliptic PDE 23 An iterative scheme can then expressed as { } { }1 *[ ] [ ] { }m mM f N f B+ = + (9) Alternatively, Eq. (9) can be rearranged as { } { } { }1 * +1 [ ]( ) { } ([ ] [ ]) { } { } m m m m m M f f B M N f Rδ + − = − − +1[ ]{ } { }m mM Rδ = (10) where { } { }+1 1{ }m m mf fδ += − is the correction vector and is an approximation to the convergence vector { }mε . It follows that +1{ } {0}mδ → ⇒ { } {0}mR → ⇒ { } {0}mε → . user 註解 讓error逼近0 Numerical Methods for Elliptic PDE 24 [Jacobi Iteration Method] (a) Point Iteration Let [P] = [I] and [A] = [L]+[D]+[U] (11) where 1 0 0 0 1 0 [ ] 0 0 1 I = 21 1 2 0 0 0 0 0 [ ] 0N N a L a a = (Identity matrix) (Lower Triangular matrix) Numerical Methods for Elliptic PDE 25 11 22 0 0 0 0 [ ] 0 0 NN a a D a = 12 1 2 0 0 0 [ ] 0 0 0 N N a a a U = (Diagonal matrix) (Upper Triangular matrix) By setting [ ] [ ] [ ] ([ ] [ ]) M D N L U = = − + (12) Numerical Methods for Elliptic PDE 26 The iteration formula takes then the form of { } { }1 *[ ] { } ([ ] [ ])m mD f B L U f+ = − + (Point Jacobi iteration) (13) For the case of the Poisson equation, Eq. (2c), the Jacobi iteration formula can be expressed as 1, , , 1, , , 1 , 1, , , 1 , 1 [ ( )]m W m S m E m N mi j i j i j i j i j i j i j i j i j i jP i j f b a f a f a f a f a + − − + += − + + + (14) Alternatively, Eq. (14) can be rearranged as 1, , , , 1, , , 1 , , , 1, , , 1 , 1 [ ( )]m m W m S m P m E m N mi j i j i j i j i j i j i j i j i j i j i j i j i jP i j f f b a f a f a f a f a f a + − − + += + − + + + + Numerical Methods for Elliptic PDE 27 ,1, , , m i jm m i j i j P i j R f f a + = + (15) with the residual at the grid (i, j) as , , , 1, , , 1 , , , 1, , , 1[ ( )] m W m S m P m E m N m i j i j i j i j i j i j i j i j i j i j i j i jR b a f a f a f a f a f− − + += − + + + + (16) user 註解 每一個點都能同時進行 Numerical Methods for Elliptic PDE 28 Point Iteration Algorithm (1) Re-order the equations to provide diagonal-dominant coefficient matrix such that the magnitude of the diagonal elements are larger than those of other elements in the same row, i.e. ii ija a≥ ( )j i≠ (2) Re-arrange each equation as of the form ,1, , , m i jm m i j i j P i j R f f a + = + (17) with the residual at the grid (i, j) as , , , 1, , , 1 , , , 1, , , 1[ ( )] m W m S m P m E m N m i j i j i j i j i j i j i j i j i j i j i j i jR b a f a f a f a f a f− − + += − + + + + (18) where the superscript m (= 0, 1, …) denotes the level of iteration. Numerical Methods for Elliptic PDE 29 (3) Assume an initial guess for 0{ }f . (4) Calculate new values of 1{ }mf + from the known value { }mf , (m ≥ 0) from Eq. (17). (5) Repeat Step (4) until a prescribed convergence criterion is satisfied. Numerical Methods for Elliptic PDE 30 (b) Line Iteration To improve iteration convergence, the point Jacobi iteration formula, Eq. (14), may be rearranged as line Jacobi iteration formulas: 1 1 1, 1, , , , 1, , , , 1 , , 1( ) W m P m E m S m N m i j i j i j i j i j i j i j i j i j i j i ja f a f a f b a f a f + + + − + − ++ + = − + (i = 2,.., (I-1)) (19a) or +1 1 +1, , 1 , , , , 1 , , 1, , 1,( ) S m P m N m W m E m i j i j i j i j i j i j i j i j i j i j i ja f a f a f b a f a f + − + − ++ + = − + (j = 2,.., (J-1)) (19b) Then, a suitable mean, such as TDMA, is used to obtain the solution of 1{ }mf + along either ith column, Eq. (19a), or jth row, Eq. (19b). Numerical Methods for Elliptic PDE 31 [Gauss-Seidel Iteration Method] (a) Point Iteration By setting [ ] ([ ] [ ]) [ ] [ ] M D L N U = + = − (20) The iteration formula takes then the form of { } { }1 *([ ]+[ ]) { } ([ ])m mD L f B U f+ = − (21) (Point Gauss-Seidel iteration) user 註解 修正jacob Numerical Methods for Elliptic PDE 32 For Eq. (2c), the Gauss-Seidel iteration formula can be expressed as 1 1 1, , , 1, , , 1 , 1, , , 1 , 1 [ ( )]m W m S m E m N mi j i j i j i j i j i j i j i j i j i jP i j f b a f a f a f a f a + + + − − + += − + + + (22) Alternatively, 1 +1 +1, , , , 1, , , 1 , , , 1, , , 1 , 1 [ ( )]m m W m S m P m E m N mi j i j i j i j i j i j i j i j i j i j i j i j i jP i j f f b a f a f a f a f a f a + − − + += + − + + + + ,1, , , m i jm m i j i j P i j R f f a + = + (23) with the residual at the grid (i, j) as +1 +1, , , 1, , , 1 , , , 1, , , 1[ ( )] m W m S m P m E m N m i j i j i j i j i j i j i j i j i j i j i j i jR b a f a f a f a f a f− − + += − + + + + (24) user 註解 該點有新值就用新的算沒有就愈舊的 Numerical Methods for Elliptic PDE 33 (b) Line Iteration The line Gauss-Seidel iteration formulas take the forms of: 1 1 1 +1, 1, , , , 1, , , , 1 , , 1( ) W m P m E m S m N m i j i j i j i j i j i j i j i j i j i j i ja f a f a f b a f a f + + + − + − ++ + = − + (i = 2,.., (I-1)) (25a) (Row iteration) or +1 1 +1 +1, , 1 , , , , 1 , , 1, , 1,( ) S m P m Nm W m E m i j i j i j i j i j i j i j i j i j i j i ja f a f a f b a f a f + − + − ++ + = − + (j = 2,.., (J-1)) (25b) (Column iteration) Numerical Methods for Elliptic PDE 34 [Relaxation Techniques] The idea of relaxation technique is to vary the convergence rate of a (point or line) iteration scheme by modifying the propagation of the correction 1 1{ } { } { }m m mf f fδ + += − through the grid system as 1 1 * * 1 { } { } { } (1 ){ } { } ({ } { }) { } { } m m m m m m m f f f f f f f f f δ ω ω ω ω δ + + + = = + − = + − = + (26) where *{ }f denotes the newly calculated value from the iteration scheme Numerical Methods for Elliptic PDE 35 ω is the relaxation factor with 0 < ω < 2, 0 1, (Under-relaxation) 1, (No relaxation) 1 2, (Over-relaxation) ω ω ω < < = < < The relaxation technique may be further classified into Stationary relaxation in which ω = constant through the iteration. Non- Stationary relaxation in which ω ≠ constant through the iteration. Numerical Methods for Elliptic PDE 36 (a) Point Relaxation The point relaxation scheme may be expressed as { } { }* *[ ] [ ] { }mM f N f B= + and 1 *{ } { } (1 ){ }m mf f fω ω+ = + − ⇒ { } { } { }1 * *[ ] [ ]( (1 ){ })=[ ] { }m m mM f M f f N f Bω ω+ = + − + { } { } { } { } { } 1 * * [ ] * * [ ] ( [ ] (1 )[ ]) { } ([ ] ([ ] [ ])) { } ([ ] [ ]) { } ((1 )[ ] [ ]) { } m m m A m m M f N M f B M M N f B M A f B A N f B ω ω ω ω ω ω ω ω ω + = = + − + = − − + = − + = − + + As a result, the point relaxation iteration scheme can be generally expressed as { } { }1 *[ ] ((1 )[ ] [ ]) { }m mM f A N f Bω ω+ = − + + (27) Numerical Methods for Elliptic PDE 37 Gauss-Seidel Relaxation Scheme Incorporating the point relaxation technique, the Gauss-Seidel point iteration scheme takes the form of * 1 1, , , 1, , , 1 , 1, , , 1 , 1 [ ( )]W m S m E m N mi j i j i j i j i j i j i j i j i j i jP i j f b a f a f a f a f a + + − − + += − + + + (28a) and 1 *, , ,(1 ) m m i j i j i jf f fω ω + = + − (28b) For 1 < ω < 2, the Gauss-Seidel iteration scheme is known as a successive over-relaxation (SOR) scheme. user 螢光標示 user 螢光標示 Numerical Methods for Elliptic PDE 38 [Ex] Numerical Methods for Elliptic PDE 39 Numerical Methods for Elliptic PDE 40 Numerical Methods for Elliptic PDE 41 Numerical Methods for Elliptic PDE 42 (b) Line Relaxation [Uni-Directional Line Relaxation] The Gauss-Seidel line relaxation formula takes first either row or column iterations as: * * * +1, 1, , , , 1, , , , 1 , , 1( ) W P E S m N m i j i j i j i j i j i j i j i j i j i j i ja f a f a f b a f a f− + − ++ + = − + (i = 2,.., (I-1)) (29a) or * * * +1, , 1 , , , , 1 , , 1, , 1,( ) S P N W m E m i j i j i j i j i j i j i j i j i j i j i ja f a f a f b a f a f− + − ++ + = − + (j = 2,.., (J-1)) (29b) Then 1 *, , ,(1 ) m m i j i j i jf f fω ω + = + − (30) Numerical Methods for Elliptic PDE 43 [Alternating Direction Line Relaxation] This scheme consists of line relaxation alternately along row and column directions as (1) Row Line Relaxation: * * * +1, 1, , , , 1, , , , 1 , , 1( ) W P E S m N m i j i j i j i j i j i j i j i j i j i j i ja f a f a f b a f a f− + − ++ + = − + (i = 2,.., (I-1)) (31a) 1 *, , ,(1 ) m m i j i j i jf f fω ω + = + − (31b) (2) Column Line Relaxation: * * * +2 +1, , 1 , , , , 1 , , 1, , 1,( ) S P N W m E m i j i j i j i j i j i j i j i j i j i j i ja f a f a f b a f a f− + − ++ + = − + (j = 2,.., (J-1)) (32b) 2 * +1, , ,(1 ) m m i j i j i jf f fω ω + = + − (32b) Numerical Methods for Elliptic PDE 44 Convergence Criteria for Iterative Methods The convergence of an iterative process is achieved when the desired accuracy criterion (or criteria) is satisfied, i.e. +1{ } {0}mδ → ⇒ +1{ } {0}mε → where { } { }+1 1{ }m m mf fδ += − is the correction vector +1 +1{ } { } { }m mf fε = − is the convergence vector Numerical Methods for Elliptic PDE 45 Convergence criteria can be specified in terms of Absolute error: 1 1, , , ,( ) m m m i j a i j i j i jf f f fδ δ + += = − Relative error: 1 1 , , , , 1 1 , , ( ) m m m i j i j i j i j r m m i j i j f f f f f f δ δ + + + + − = = Numerical Methods for Elliptic PDE 46 [Absolute Convergence Criteria] (a) Infinity norm: 1, , max{ } m m a i j i jf f fδ ε + ∞ = − ≤ (b) The 1st norm: 1, ,1{ } m m a i j i j j i f f fδ ε+= − ≤∑∑ (c) The 2nd norm: 1 2 1/ 2, ,2{ } [ ( ) ] m m a i j i j j i f f fδ ε+= − ≤∑∑ Numerical Methods for Elliptic PDE 47 [Relative Convergence Criteria] (a) Infinity norm: 1 , , max 1 , max { } m m i j i j r m i j f f f f δ ε + ∞ + − = ≤ (b) The 1st norm: 1 , , 11 , { } m m i j i j r m i jj i f f f f δ ε + + − = ≤∑∑ (c) The 2nd norm: 1 , , 2 1/ 2 12 , { } [ ( ) ] m m i j i j r m i jj i f f f f δ ε + + − = ≤∑∑ Numerical Methods for Elliptic PDE 48 Comparison of the Methods for Solving Systems of Linear Algebraic Equations Method Approx. Max. No. of Equations Stability Precision Breadth of Applications Programming Effort Comments Cramer’s rule 3 - Affected by round-off error Limited - Excessive computational effort required for more than three equations LU Decomposition 100 - Affected by round-off error General Moderate Preferred elimination method Gauss-Seidel 1000 May not converge if not diagonally dominant Excellent Appropriate only for diagonally dominant systems Easy Numerical Methods for Elliptic PDE 49 END
Compartilhar