麻省理工学院：偏微分方程式数字方法（英文版）：lec5

SMA-HPC ?2002 NUS 16.920J/SMA 5212 Numerical Methods for PDEs Thanks to Franklin Tan Finite Differences: Parabolic Problems B. C. Khoo Lecture 5 SMA-HPC ?2002 NUS Outline ? Governing Equation ? Stability Analysis ? 3 Examples ? Relationship between σ and λh ? Implicit Time-Marching Scheme ? Summary 2 SMA-HPC ?2002 NUS [ 2 2 0, u x t υ π ? = ? ? Governing Equation 3 Consider the Parabolic PDE in 1-D 0 subject to at 0, atuu x u u x π π= = = ? If υ≡ viscosity → Diffusion Equation ? If υ≡ thermal conductivity → Heat Conduction Equation 0x = x π= 0 u π ( ) , ux t = ] u x ? ∈ = u ? SMA-HPC ?2002 NUS Stability Analysis Discretization 4 Keeping time continuous, we carry out a spatial discretization of the RHS of There is a total of 1 grid points such that , 0,1, 2,...., j N j x j + =? = 0x = x π= 0 x 1 x 3 x 1N x ? N x 2 2 u t υ ? ? = ? ? x N u x SMA-HPC ?2002 NUS Stability Analysis Discretization 5 2 1 2 2 2 ( j j j u u u Ox x + ? ? ? = ? ? ? ? 2 2 Use the Central Difference Scheme for u x ? ? which is second-order accurate. ? Schemes of other orders of accuracy may be constructed. 1 2 ) j u x ? + ? + ? ? ? SMA-HPC ?2002 NUS 6 Stability Analysis Discretization We obtain at 0 0 Note that we need not evaluate at and since and are given as boundary conditions. N N u x x x u = = 1 1 2 2 : 2 ) o du x u u dt x υ = + ? 2 2 2 3 2 : 2 ) du x u u dt x υ = + ? 1 2 : 2 ) j j j j du x u u dt x υ ? + = + ? 1 1 1 2 : 2 ) N N N N du x u u dt x υ ? ? ? = + ? x u 1 ( u ? 1 ( u ? 1 ( j u ? 2 ( N u ? ? SMA-HPC ?2002 NUS 7 Stability Analysis Matrix Formulation Assembling the system of equations, we obtain 1 1 2 2 2 2 1 1 2 21 0 1 1 0 1 1 1 0 1 o jj N N N du u u dt x du u dt udu x dt u du u x dt υ υ υ ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? = ? ? ? ? ? ?? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?? ? ? ? ? ?? ? ? ? ? ? ? ? 0 0 A 2 2 2 ? ??? ??? ??? ??? ??? ??? ??? +? ?? ??? ??? ??? ??? ??? ??? ??? ? ? ? ? SMA-HPC ?2002 NUS 8 Stability Analysis PDE to Coupled ODEs Or in compact form du Aub dt = null null null We have reduced the 1-D PDE to a set of Coupled ODEs! [ ] 1 1 where T N u u u ? = null 2 0 0 T o u b x υ ? ? = ? ? ? ? ? null + 2 u 2 0 N u x υ ? SMA-HPC ?2002 NUS 9 Stability Analysis Eigenvalue and Eigenvector of Matrix A If A is a nonsingular matrix, as in this case, it is then possible to find a set of eigenvalues { 1 1 , ,...., ,...., j λ λ λ λ ? = For each eigenvalue , we can evaluate the eigenvector consisting of a set of mesh point values , i.e. j j j i V v λ 1 1 T j j j N V v v ? ? = ? ( from det 0.A λ? } 2 N λ 2 j v ? ? ) I = SMA-HPC ?2002 NUS 10 Stability Analysis Eigenvalue and Eigenvector of Matrix A The ( 1) ( 1) matrix formed by the ( 1) columns diagonalizes the matrix by j N E N V ?× ? ? 1 E AE ? = Λ 1 2 1 where N λ λ λ ? ? ? ? ? ? ? ? ? Λ= ? ? ? ? ? ? ? ? 0 0 N A SMA-HPC ?2002 NUS 11 Stability Analysis Coupled ODEs to Uncoupled ODEs Starting from du Au b dt = + null null null 1 Premultiplication by yieldsE ? 1 1 du E EAu E b dt ? ? = null null null ( 1 1 1 du E EA EE u E b dt ? ? ? = null null null I ( 1 1 1 du E EAE E u E b dt ? ? ? = null null null Λ 1 ? + ) 1 ? + ) 1 ? + SMA-HPC ?2002 NUS 12 Stability Analysis Coupled ODEs to Uncoupled ODEs 1 Let and , we haveUE u F E b ? = null null null d U F dt = Λ+ nullnullnullnull null which is a set of Uncoupled ODEs! 1 1 du E u E b dt ? ? =Λ + null null null Continuing from 1? = null U 1 E ? SMA-HPC ?2002 NUS 13 Stability Analysis Coupled ODEs to Uncoupled ODEs Expanding yields Since the equations are independent of one another, they can be solved separately. The idea then is to solve for and determineU EU= null null null 1 11 1 dU U dt λ= + 2 2 2 dU U dt λ= + j j j dU U dt λ= + 1 1 1 N N N dU U dt λ ? ? ? = u F 2 F j F 1 N F ? + SMA-HPC ?2002 NUS 14 Considering the case of independent of time, for the general equation, th b j null Stability Analysis Coupled ODEs to Uncoupled ODEs 1 jt j j j U e F λ λ = is the solution for j = 1,2,….,N–1. Evaluating, ( ) 1 t uEU E ce E E b λ ? = ? Λ nullnullnullnullnull null null Complementary (transient) solution Particular (steady-state) solution ( 11 1 1 where j N T t tt t j ce c e c e c e c e λ λλ λ ? ? ? ? = ? ? nullnullnullnullnull j c ? 1? = null ) 2 2 t N λ SMA-HPC ?2002 NUS We can think of the solution to the semi-discretized problem 15 Stability Analysis Stability Criterion ( ) 1 t uE ce E E b λ ? = Λ nullnullnullnullnull null null This is the criterion for stability of the space discretization (of a parabolic PDE) keeping time continuous. Since the transient solution must decay with time, for all j ( ) Real 0 j λ ≤ Each mode contributes a (transient) time behaviour of the form to the time-dependent part of the solution. j t j e λ as a superposition of eigenmodes of the matrix operator A. 1? ? SMA-HPC ?2002 NUS 16 Stability Analysis Use of Modal (Scalar) Equation It may be noted that since the solution is expressed as a contribution from all the modes of the initial solution, which have propagated or (and) diffused with the eigenvalue , and a contribution fr j u λ null om the source term , all the properties of the time integration (and their stability properties) can be analysed separately for each mode with the scalar equation j b j dU UF dt λ ? = ? ? ? + ? ? SMA-HPC ?2002 NUS 17 Stability Analysis Use of Modal (Scalar) Equation The spatial operator A is replaced by an eigenvalue λ, and the above modal equation will serve as the basic equation for analysis of the stability of a time-integration scheme (yet to be introduced) as a function of the eigenvalues λ of the space-discretization operators. This analysis provides a general technique for the determination of time integration methods which lead to stable algorithms for a given space discretization. SMA-HPC ?2002 NUS 18 1 11 1 12 2 2 21 1 22 2 du au a u dt du au a u dt = = 1 1 12 2 1 22 Let , u a u u a ?? ? ? = ? ?? ? ? ?? ? ? null du Au dt = null null Consider a set of coupled ODEs (2 equations only): Example 1 Continuous Time Operator + + 1 2 a A a = SMA-HPC ?2002 NUS 19 Proceeding as before, or otherwise (solving the ODEs directly), we can obtain the solution 1 1 1 11 2 12 2 21 2 22 t t u e c e u e c e λ λ λ ξ ξ = = 11 21 1 21 22 1 where and are eigenvalues of and and are eigenvectors pertaining to and respectively. A ξ λ ξ λ ?? ? ? ?? ? ? ?? ? ? () j As the transient solution must decay with time, it is imperative that Real 0 for 1, 2.jλ ≤ Example 1 Continuous Time Operator 2 2 1 1 t t c c λ ξ ξ + + 2 2 ξ λ ξ λ = SMA-HPC ?2002 NUS 20 Suppose we have somehow discretized the time operator on the LHS to obtain 1 1 1 1 12 2 1 2 1 1 22 2 n n n n u u a u u u a u ? ? ? ? = = where the superscript n stands for the n th time level, then 1 11 12 1 21 22 where and T n n n a u u u u u A a ? ? ? ? = = ? ? ? ? ? null null Since A is independent of time, 1 0 .... n n n u u AAu A u ? = = = null null null Example 1 Discrete Time Operator 1 1 1 2 n n a a + + 2 n n a A a ? = ? null 2 n A ? = null SMA-HPC ?2002 NUS 21 Example 1 Discrete Time Operator 10 1 2 0 where = 0 n n n n u E u λ λ ? ? ? =Λ Λ ? ? ? ? nullnullnull null ' 1 1 11 1 2 12 2 ' 2 1 21 1 2 22 2 n n n n u c u c λξ λ ξ λξ λ ξ = = 1 10 2 ' where are constants. ' c Eu c ? ?? = ?? ?? nullnullnull 1 1 1 0 , .... n AE E u E E E E E u ? ? ? =Λ =Λ ? Λ ? ? Λ ? nullnullnull null As A A A n E ' ' n n c c + + 1 E ? SMA-HPC ?2002 NUS 22 Example 1 Comparison Comparing the solution of the semi-discretized problem where time is kept continuous [ 1 2 1 1 12 1 2 1 22 t t u e c u e λ λ ξ ξ ? ? ?? ? ? = ? ? ?? ? ? ?? ? ? ? ? to the solution where time is discretized [ 1 1 12 1 1 2 1 22 2 ' n n n u c u ξ ξ ξ λ ? ? ?? ? ? = ? ? ?? ? ? ?? ? ? ? ? ? ? coTh nte di inuofference equation where time is hus exponential solution as . t e λ The difference equation where time is hasdiscretized power solution . n λ ] 1 2 2 c ξ ξ ] 1 2 2 ' c λ ξ SMA-HPC ?2002 NUS 23 Example 1 Comparison In equivalence, the transient solution of the difference equation must decay with time, i.e. for this particular form of time discretization. 1 n λ < SMA-HPC ?2002 NUS 24 Consider a typical modal equation of the form t j du uae dt μ λ ? = ? ? where is the eigenvalue of the associated matrix . j Aλ (For simplicity, we shall henceforth drop the subscript j). We shall apply the “leapfrog” time discretization scheme given as Substituting into the modal equation yields ( 1 2 n t tnh u uae h μ λ + = ? = n n u e μ λ= 1 where 2 n du u u h dt h + ? = ? Example 2 Leapfrog Time Discretization ? + ? ? ) 1 n u ? + h a+ 1 n t ? = SMA-HPC ?2002 NUS Solution of u consists of the complementary solution c n , and the particular solution p n , i.e. u n =c n + p n There are several ways of solving for the complementary and particular solutions. shift operator S and characteristic polynomial. The time shift operator S operates on c n such that Sc n = c n+1 S 2 c n = S(Sc n ) = Sc n+1 = c n+2 25 ( 1 1 2 2 2 n n n n n n hn u u e u h u u ha e h μ λ + + ? = ? ? ? = Example 2 Leapfrog Time Discretization Time Shift Operator One way is through use of the ) 1 1 n h u a μ λ ? ? + SMA-HPC ?2002 NUS 26 The complementary solution c n satisfies the homogenous equation 1 2 2 2 2 1 ( ) 0 ( 1) 0 n n n n n n n c c c c Sc h c S Sc h Sc c S c S S S λ λ λ λ + ? = ? = ? = ? = 2 () ( 2 1) 0pS S h Sλ= ? = characteristic polynomial Example 2 Leapfrog Time Discretization Time Shift Operator 1 0 0 2 2 n n n h h ? ? ? ? ? ? SMA-HPC ?2002 NUS 27 The complementary solution to the modal equation would then be 11 2 2 n n c βσ σ= The particular solution to the modal equation is 2 2 2 hn h n h ahe e p e e μ μ λ = ? ? Combining the two components of the solution together, () ( ) n n u p= ( ) ( ) 22 22 1 2 2 1 2 hn h n h ahe e h h h e e μ μ βλ λ β λ λ λ ? ? = + + ? + + ? ? ? ? ? ? The solution to the characteristic polynomial is 22 () 1h h hσλ λ λ== ± + σ 1 and σ 2 are the two roots Example 2 Leapfrog Time Discretization Time Shift Operator n β+ 1 h h μ μ n c + 2 1 1 n h h h μ μ ? ? + ?? ? ? S SMA-HPC ?2002 NUS 28 For the solution to be stable, the transient (complementary) solution must not be allowed to grow indefinitely with time, thus implying that is the stability criterion for the leapfrog time discretization scheme used above. ( ) ( ) 22 1 22 2 1 1 h h σ λ σ λ = + < = ? < Example 2 Leapfrog Time Discretization Stability Criterion 1 1 h h λ λ + + SMA-HPC ?2002 NUS 29 Im(σ ) Re(σ ) -1 1 Region of Stability The stability diagram for the leapfrog (or any general) time discretization scheme in the σ-plane is Example 2 Leapfrog Time Discretization Stability Diagram SMA-HPC ?2002 NUS 30 In particular, by applying to the 1-D Parabolic PDE 2 2 u t υ ? ? = ? ? the central difference scheme for spatial discretization, we obtain which is the tridiagonal matrix. Example 2 Leapfrog Time Discretization 2 21 1 1 1 1 A x υ ? ? ? ? ? ? ? ? ? ? = ? ? ? ? ? ? ? ? ? ? ? ? ? ? 0 0 u x 2 2 SMA-HPC ?2002 NUS 31 Example 2 Leapfrog Time Discretization According to analysis of a general triadiagonal matrix B(a,b,c), the eigenvalues of the B are 2 2 cos , 1,..., 1 22 cos j j j b c j N N j N π λ π λ ? = + ? ? ? ? ? =? + ? ? ? ? ? The most “dangerous” mode is that associated with the eigenvalue of largest magnitude max 2 4 x υ λ =? ? ( ( 2 max 1 ax max 2 max 2 ax max 1 1 h h h h σλ λ λ σλ λ λ = + = + i.e. which can be plotted in the absolute stability diagram. a x υ ? = ? ? ? ? ? ? ? ? ) ) 2 m 2 m h h + ? SMA-HPC ?2002 NUS 32 Example 2 Leapfrog Time Discretization Absolute Stability Diagram for σ As applied to the 1-D Parabolic PDE, the absolute stability diagram for σ is Region of stability Unit circle Region of instability σ 2 at h = ?t = 0 σ 2 with h increasing σ 1 with h increasing σ 1 at h = ?t = 0 Re(σ) Im(σ) SMA-HPC ?2002 NUS 33 Stability Analysis Some Important Characteristics Deduced A few features worth considering: 1. Stability analysis of time discretization scheme can be carried out for all the different modes . 2. If the stability criterion for the time discretization scheme is j λ valid for all modes, then the overall solution is stable (since it is a linear combination of all the modes). 3. When there is more than one root , then one of them is the principal root which represents σ ( 0 an approximation to the physical behaviour. The principal root is recognized by the fact that it tends towards one as 0, i.e. lim 1. (The other roots are spurious, which affect the stability h h λ λ λ → → but not the accuracy of the scheme.) ) h σ = SMA-HPC ?2002 NUS 34 Stability Analysis Some Important Characteristics Deduced 1 4. By comparing the power series solution of the principal root to , one can determine the order of accuracy of the time discretization scheme. In this example of leapfrog time discretization, 1 h e h λ σ = ( ( 1 22 22 44 2 22 1 22 1 . 1 2 1 2 ! 1 ... 2 and compared to 1 .. 2! is identical up to the second order of . Hence, the above scheme is said to be second-order accurate. h h h h h h h e h λ λ λ λ λ σ λ λ λ ? + + + + =+ + + =+ + + λ + ) ) 1 2 . 2 . h h λ λ = SMA-HPC ?2002 NUS 35 Analyze the stability of the explicit Euler-forward time discretization 1n du u u dt t + ? = ? as applied to the modal equation du u dt λ= 1 1 Substituting where into the modal equation, we obtain (1 ) 0 n n du u h h t dt u uλ + + = = ? ? + Euler-Forward Time Discretization Stability Analysis Example 3 n n n u h + = SMA-HPC ?2002 NUS 36 Making use of the shift operator S 1 (1 ) (1 ) [ (1 )] 0 n n n n c c Sc h c S h cλ λ + ?+ = ? + = ?+ = Therefore ( ) 1 and n h c σ λ βσ = + = characteristic polynomial The Euler-forward time discretization scheme is stable if Euler-Forward Time Discretization Stability Analysis Example 3 1 hσ ≡ + or bounded by 1 s.t. 1 in the -plane.h λ σ λ=? < n h λ n h λ 1λ < hσ SMA-HPC ?2002 NUS 37 Euler-Forward Time Discretization Stability Diagram Example 3 Im(λh) 0-1-2 Unit Circle Region of Stability Re(λh) The stability diagram for the Euler-forward time discretization in the λh-plane is SMA-HPC ?2002 NUS 38 Euler-Forward Time Discretization Absolute Stability Diagram Example 3 max 2 4 As applied to the 1-D Parabolic PDE, x υ λλ= ? ? max 2 The stability limit for largest h λ ? ≡? = 1-1 σ leaves the unit circle at λh = ?2 σ at h (=?t) = 0 Re(σ) Im(σ) σ with h increasing = t SMA-HPC ?2002 NUS 39 Relationship between σ and λh σ = σ(λh) Thus far, we have obtained the stability criterion of the time discretization scheme using a typical modal equation. generalize the relationship between σ and λh as follows: ? Starting from the set of coupled ODEs du Aub dt = + null null null ? Apply a specific time discretization scheme like the “leapfrog” time discretization as in Example 2 1 2 n du u u dt h + ? ? = We can 1 n SMA-HPC ?2002 NUS 40 Relationship between σ and λh σ = σ(λh) ? The above set of ODEs becomes 1 2 n n n u Au h + ? = + null null null ? Introducing the time shift operator S 1 2 2 n n n n n u Su hAu hb S SS A Iu b h ? = + ? ? ? ? ? ? null null null null null 1 1 Premultiplying on the LHS and RHS and introducing operating on n E IEE u ? ? = i null 1 1 1 1 2 n SS E AE E E E u E b h ? ? ? ? ? ? ? ? ? ? null null Λ 1 n u b ? null 2 n + ? = ? ? null 1 ? ? = ? ? SMA-HPC ?2002 NUS 41 Relationship between σ and λh σ = σ(λh) ? Putting 1 , n n n U u F E b ? = null null null 1 2 j j SS U h λ ? ? ? ? =? ? ? we obtain 1 1 2 n SS E U F h ? ? ? ? Λ ? ? ? ? null null 1 2 SS h ? ? i.e. 1 2 n SS U h ? ? ? Λ? = ? ? ? null null which is a set of uncoupled equations. Hence, for each j, j = 1,2,….,N-1, 1 n E ? = null j F ? ? ? n E ? = ? ? n F ? ? ? SMA-HPC ?2002 NUS 42 Relationship between σ and λh σ = σ(λh) Note that the analysis performed above is identical to the analysis carried out using the modal equation j dU UF dt λ ? = ? ? All the analysis carried out earlier for a single modal equation is applicable to the matrix after the appropriate manipulation to obtain an uncoupled set of ODEs. Each equation can be solved independently for and the 's can then be coupled through . th n n n j j U u EU= null null ? + ? ? n j U SMA-HPC ?2002 NUS 43 Relationship between σ and λh σ = σ(λh) 1. Uncoupling the set, 2. Integrating each equation in the uncoupled set, 3. Re-coupling the results to form the final solution. These 3 steps are commonly referred to as the ISOLATION THEOREM Hence, applying any “consistent” numerical technique to each equation in the set of coupled linear ODEs is mathematically equivalent to SMA-HPC ?2002 NUS 44 Implicit Time- Marching Scheme Thus far, we have presented examples of explicit time-marching methods and these may be used to integrate weakly stiff equations. Implicit methods are usually employed to integrate very stiff ODEs efficiently. solution of a set of simultaneous algebraic equations at each time-step (i.e. matrix inversion), whilst updating the variables at the same time. Implicit schemes applied to ODEs that are inherently stable will be unconditionally stable or A-stable. However, use of implicit schemes requires SMA-HPC ?2002 NUS 45 Implicit Time- Marching Scheme Euler-Backward Consider the Euler-backward scheme for time discretization 1 1 n n du u u dt h + + ? ? = ? ? t du uae dt μ λ= ( ( ( 1 1 1 1 1 1 n n n n n u u e h hu u ahe μ μ λ λ + + + + + ? ? ? = ? ? ? = Applying the above to the modal equation for Parabolic PDE yields n ? ? ? + ) ) ) n h h n u a+ ? SMA-HPC ?2002 NUS 46 Implicit Time- Marching Scheme Euler-Backward Applying the S operator, ( ) ( 1 1 n n hS u ahe μ λ + ? ? = ? the characteristic polynomial becomes ( ) ( ) ( ) 1 0S Sσ ? Ρ Ρ = ? ? = ? The principal root is therefore 22 1 which, upon comparison with 1 .... , is only 2 first-order accurate. h e h λ λ =+ + + 22 1 1 .... 1 h h σ λ λ = =+ + + ? The solution is ( ( 1 1 1 1 n u n h ahe U h e μ μ β λ + ? = ? ? ? ? ) 1 h ?? ? 1 h λ ?= ? h λ hλ ) ) 1 h h λ ? + ? ? ? SMA-HPC ?2002 NUS 47 Implicit Time- Marching Scheme Euler-Backward For the Parabolic PDE, λ is always real and < 0. Therefore, the transient component will always tend towards zero for large n irregardless of h (≡?t). The time-marching scheme is always numerically stable. In this way, the implicit Euler/Euler-backward time discretization scheme will allow us to resolve different time-scaled events with the use of different time-step sizes. scaled events, and then a large time-step size used for the longer time-scaled events. h max . A small time-step size is used for the short time- There is no constraint on SMA-HPC ?2002 NUS 48 Implicit Time- Marching Scheme Euler-Backward However, numerical solution of u requires the solution of a set of simultaneous algebraic equations or matrix inversion, which is computationally much more intensive/expensive compared to the multiplication/ addition operations of explicit schemes. SMA-HPC ?2002 NUS Summary 49 ? Stability Analysis of Parabolic PDE null Uncoupling the set. null Integrating each equation in the uncoupled set → modal equation. null Re-coupling the results to form final solution. ? Use of modal equation to analyze the stability |σ(λh)| < 1. ? Explicit time discretization versus Implicit time discretization.