16.920J/SMA 5212 Numerical Methods for Partial Differential Equations Lecture 5 Finite Differences: Parabolic Problems B. C. Khoo Thanks to Franklin Tan 19 February 2003 16.920J/SMA 5212 Numerical Methods for PDEs 2 OUTLINE ? Governing Equation ? Stability Analysis ? 3 Examples ? Relationship between σ and λh ? Implicit Time-Marching Scheme ? Summary Slide 2 GOVERNING EQUATION Consider the Parabolic PDE in 1-D a0 If υ ≡ viscosity → Diffusion Equation a0 If υ ≡ thermal conductivity → Heat Conduction Equation Slide 3 STABILITY ANALYSIS Discretization Keeping time continuous, we carry out a spatial discretization of the RHS of [ ]2 2 0,u u xt xυ pi? ?= ∈? ? 0subject to at 0, at u u x u u xpi pi= = = = 0x = x pi= 0u upi ( ), ?u x t = 2 2 u u t xυ ? ?= ? ? 0x = x pi= 0x 1x 2x 1Nx ? Nx 16.920J/SMA 5212 Numerical Methods for PDEs 3 Slide 4 STABILITY ANALYSIS Discretization which is second-order accurate. ? Schemes of other orders of accuracy may be constructed. Slide 5 Construction of Spatial Difference Scheme of Any Order p The idea of constructing a spatial difference operator is to represent the spatial differential operator at a location by the neighboring nodal points, each with its own weightage. The order of accuracy, p of a spatial difference scheme is represented as ( )pO x? . Generally, to represent the spatial operator to a higher order of accuracy, more nodal points must be used. Consider the following procedure of determining the spatial operator j du dx a0 a1 a2 a3 a4 a5 up to the order of accuracy ( )2O x? : There is a total of 1 grid points such that , 0,1,2,...., jN x j x j N + = ? = 2 2Use the Central Difference Scheme for u x ? ? 2 1 1 2 2 2 2 ( )j j j j u u uu O x x x + ?? + a6 a7 ? = + ?a8 a9 ? ?a10 a11 j?2 ? j?1 ? j ? j+1 ? j+2 ? j du dx a12 a13 a14 a15 a16 a17 16.920J/SMA 5212 Numerical Methods for PDEs 4 1. Let j du dx a0 a1 a2 a3 a4 a5 be represented by u at the nodes j?1, j, and j+1 with 1α? , 0α and 1α being the coefficients to be determined, i.e. ( )1 1 0 1 1 pj j j j du u u u O x dx α α α? ? + a6 a7 + + + = ?a8 a9a10 a11 2. Seek Taylor Expansions for 1ju ? , ju and 1ju + about ju and present them in a table as shown below. (Note that p is not known a priori but is determined at the end of the analysis when the α’s are made known.) uj uj′ uj′′ uj′′′ ju ′ 0 1 0 0 1 1juα? ? 1α? 1x α??? ? 2 1 1 2 x α?? ? 3 1 1 6 x α?? ? ? 0 juα 0α 0 0 0 1 1juα + 1α 1x α? ? 2 1 1 2 x α? ? 3 1 1 6 x α? ? 1 1 k j k j k k u uα = + =? ′ + a12 1S 2S 3S 4S ( 1 ) This column consists of all the terms on the LHS of (1). Each cell in this row comprises the sum of its corresponding column. 16.920J/SMA 5212 Numerical Methods for PDEs 5 where ∴ 1 1 2 3 4 1 .... k j k j k k u u S S S Sα = + =? ′ + = + + + +a0 3. Make as many iS ’s as possible vanish by choosing appropriate kα ’s. In this instance, since we have three unknowns 1α? , 0α and 1α , we can therefore set: 1 2 3 0 0 0 S S S = = = (Note that in the Taylor Series expansion, one starts off with the lower-order terms and progressively obtain the higher-order terms. We have deliberately set the iS pertaining to the lower-order terms to zero, thereafter followed by increasingly higher-order terms.) Hence, 1 0 1 01 1 1 11 0 1 1 0 1 0 x α α α ? a1 a2 a1 a2a3a1 a2 a4 a5 a4 a5 a4 a5 a4 a5 ? = ?a4 a5 a4 a5 ?a4 a5 a4 a5 a4 a5 a6 a7a3a6 a7 a4 a5 a6 a7 Solving the system of equations, we obtain 1 0 1 1 2 0 1 2 x x α α α ? = ? = = ? ? ( ) ( ) 1 1 0 1 2 1 1 2 2 3 1 1 3 3 4 1 1 1 1 1 2 2 1 1 6 6 j j j j S u S x x u S x x u S x x u α α α α α α α α α ? ? ? ? = + + ′= ? ? ? + ? ? a8 a9 ′′= ? ? + ? ?a10 a11 a12 a13 a8 a9 ′′′= ? ? ? + ? ?a10 a11 a12 a13 16.920J/SMA 5212 Numerical Methods for PDEs 6 4. Substituting the kα ’s into 1 1 2 3 4 1 .... k j k j k k u u S S S Sα = + =? ′ + = + + + +a0 yields ( ) 21 11 12 6j j j ju u u x ux + ?′ ′′′? ? = ? ? +? higher-order terms In other words, ( )1 1 2 ....2j jj j u uduu O x dx x + ?? a1 a2 ′ = = + ? +a3 a4 ?a5 a6 i.e. the above representation is accurate up to ( )2O x? . Some useful points to note: 1. These 4 steps are the general procedure used to obtain the representation of the spatial operator up to the order of accuracy ( )pO x? . 2. For other spatial operators, say 2 2 j d u dx a7 a8 a9 a10 a11 a12 , we simply replace j du dx a13 a14 a15 a16 a17 a18 in (1) with the said spatial operator. 3. For one-sided representations, one can choose nodal points , 0j ku k+ ≥ . This may be important especially for representations on a boundary. For example ( )0 1 1 2 2 .... pj j j j du u u u O x dx α α α+ + a19 a20 + + + + = ?a21 a22a23 a24 One possibility is ( )1 2 23 42j j j j u u udu O x dx x + +? + a25 a26 + = ?a27 a28 ?a29 a30 which is also second-order accurate. (We can also use a similar procedure to construct the finite difference scheme of Hermitian type for a spatial operator. This is not covered here). ( ), 2pO x p? = 16.920J/SMA 5212 Numerical Methods for PDEs 7 STABILITY ANALYSIS Discretization We obtain at 11 1 22: ( 2 )odux u u udt xυ= ? +? 22 1 2 32: ( 2 )dux u u udt xυ= ? +? 1 12: ( 2 )jj j j jdux u u udt xυ ? += ? +? 11 2 12: ( 2 )NN N N Ndux u u udt xυ?? ? ?= ? +? 0 0 Note that we need not evaluate at and since and are given as boundary conditions. N N u x x x x u u = = Slide 6 STABILITY ANALYSIS Matrix Formulation Assembling the system of equations, we obtain Slide 7 1 1 2 2 2 2 1 1 2 2 1 01 2 1 0 1 2 1 1 0 1 2 o jj NN N du uu dt x du udt udu x dt udu u xdt υ υ υ ? ? a0 a1 a0 a1 ? a0 a1a0 a1a2 a3 a2 a3? a2 a3a2 a3a2 a3 a2 a3 a2 a3a2 a3a2 a3 a2 a3 a2 a3a2 a3a2 a3 ? a2 a3a2 a3a2 a3a2 a3 a2 a3 a2 a3a2 a3a2 a3 a2 a3 a2 a3a2 a3a2 a3 a2 a3 = +a2 a3a2 a3a2 a3 a2 a3?? a2 a3a2 a3a2 a3 a2 a3 a2 a3a2 a3a2 a3 a2 a3 a2 a3a2 a3a2 a3 a2 a3 a2 a3a2 a3a2 a3 a2 a3 a2 a3a2 a3a2 a3 a2 a2 a3a2 a3a2 a3 a2 a2 a3a2 a3 ?a2 a3 a4 a5a6a4 a5 a2a2 ?a4 a5a2 a3a4 a5 a3 a3 a3 a3 0 0 A 16.920J/SMA 5212 Numerical Methods for PDEs 8 STABILITY ANALYSIS PDE to Coupled ODEs Or in compact form We have reduced the 1-D PDE to a set of Coupled ODEs! Slide 8 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 2 1, ,...., ,....,j Nλ λ λ λ λ ?= ( )from det 0.A Iλ? = For each eigenvalue , we can evaluate the eigenvector consisting of a set of mesh point values , i.e. j j j i V v λ Slide 9 STABILITY ANALYSIS Eigenvalue and Eigenvector of Matrix A The ( 1) ( 1) matrix formed by the ( 1) columns diagonalizes the matrix byj N N E N V A ? × ? ? 1E AE? = Λ [ ]1 2 1where TNu u u u ?=a0 2 20 0 0 T o Nu ub x x υ υa1 a2= a3 a4? ?a5 a6 a7 du Au b dt = + a8 a8a9a8 1 2 1 Tj j j j NV v v v ? a10 a11 =a12 a13 16.920J/SMA 5212 Numerical Methods for PDEs 9 Slide 10 STABILITY ANALYSIS Coupled ODEs to Uncoupled ODEs Starting from du Au bdt = + a8 a8a9a8 1Premultiplication by yieldsE? 1 1 1duE E Au E bdt? ? ?= + a0 a0 a0 Slide 11 STABILITY ANALYSIS Coupled ODEs to Uncoupled ODEs Continuing from 1 1 1duE E u E bdt? ? ?= Λ + a0 a0 a0 1 1Let and , we haveU E u F E b? ?= = a1a1 a1a1 1 2 1 where N λ λ λ ? a2 a3 a4 a5 a4 a5 a4 a5 Λ = a4 a5 a4 a5 a4 a5 a6 a7 0 0 ( )1 1 1 1duE E A EE u E bdt? ? ? ?= + a8 a8 a8 I Λ ( )1 1 1 1duE E AE E u E bdt? ? ? ?= + a9 a9 a9 16.920J/SMA 5212 Numerical Methods for PDEs 10 d U U Fdt = Λ + a0 a1 a0 a1a1 which is a set of Uncoupled ODEs! Slide 12 STABILITY ANALYSIS Coupled ODEs to Uncoupled ODEs Expanding yields 1 1 1 1dU U Fdt λ= + 2 2 2 2dU U Fdt λ= + j j j jdU U Fdt λ= + 1 1 1 1N N N NdU U Fdt λ? ? ? ?= + Since the equations are independent of one another, they can be solved separately. The idea then is to solve for and determine U u EU= a2 a2 a2 Slide 13 STABILITY ANALYSIS Coupled ODEs to Uncoupled ODEs Considering the case of independent of time, for the general equation,th b j a3 1jtj j j j U c e Fλ λ= ? is the solution for j = 1,2,….,N?1. 16.920J/SMA 5212 Numerical Methods for PDEs 11 Evaluating, ( ) 1 1tu EU E ce E E bλ ? ?= = ? Λ a0a1a0a1a0 a0a2 a2 a2a2 ( ) 11 21 2 1where j N Tt tt tt j Nce c e c e c e c eλ λλ λλ ??a3 a4= a5 a6a7a1a7a1a7a8a7a9 The stability analysis of the space discretization, keeping time continuous, is based on the eigenvalue structure of A. The exact solution of the system of equations is determined by the eigenvalues and eigenvectors of A. Slide 14 STABILITY ANALYSIS Coupled ODEs to Uncoupled ODEs We can think of the solution to the semi-discretized problem as a superposition of eigenmodes of the matrix operator A. Each mode contributes a (transient) time behaviour of the form to the time-dependent part of the solution.jt j eλ Since the transient solution must decay with time, ( )Real 0jλ ≤ for all j This is the criterion for stability of the space discretization (of a parabolic PDE) keeping time continuous. Slide 15 Complementary (transient) solution Particular (steady-state) solution ( ) 1 1tu E ce E E bλ ? ?= ? Λa10a1a10a1a10a8a10a11a11 a11 16.920J/SMA 5212 Numerical Methods for PDEs 12 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 frj u λ a0 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 jb Slide 16 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. Slide 17 EXAMPLE 1 Continuous Time Operator Consider a set of coupled ODEs (2 equations only): 1 11 1 12 2 2 21 1 22 2 du a u a u dt du a u a u dt = + = + 1 11 12 2 21 22 Let , u a a duu A Auu a a dt a1a3a2 a1 a2 = = a4 =a5 a6 a5 a6 a7a3a8 a7 a8 a9 a9 a9 Slide 18 j dU U F dt λ a10 a11 = +a12 a13a14 a15 16.920J/SMA 5212 Numerical Methods for PDEs 13 EXAMPLE 1 Continuous Time Operator Proceeding as before, or otherwise (solving the ODEs directly), we can obtain the solution 1 2 1 2 1 1 11 2 12 2 1 21 2 22 t t t t u c e c e u c e c e λ λ λ λ ξ ξ ξ ξ = + = + 11 21 1 2 21 22 1 2 where and are eigenvalues of and and are eigenvectors pertaining to and respectively. A ξ ξλ λ ξ ξ λ λ a0a2a1 a0 a1 a3 a4 a3 a4 a5a2a6 a5 a6 ( )j As the transient solution must decay with time, it is imperative that Real 0 for 1, 2.jλ ≤ = Slide 19 EXAMPLE 1 Discrete Time Operator Suppose we have somehow discretized the time operator on the LHS to obtain 1 1 1 11 1 12 2 1 1 2 21 1 22 2 n n n n n n u a u a u u a u a u ? ? ? ? = + = + where the subscript n stands for the nth time level, then 1 11 121 2 21 22 where and Tn n n n n a au Au u u u A a a? a7 a8a7 a8= = =a9 a10a11 a12 a11 a12 a13 a13 a13 Since A is independent of time, 1 2 0....n n n nu Au AAu A u? ?= = = = a14 a14 a14 a14 In later examples, we shall apply specific time discretization schemes such as the “leapfrog” and Euler-forward time discretization schemes. Slide 20 16.920J/SMA 5212 Numerical Methods for PDEs 14 EXAMPLE 1 Discrete Time Operator As 1 0 1 2 0 where = 0 nn n n nu E E u λ λ ? a0 a1 = Λ Λ a2 a3 a4 a5 a6a7a6a8 a8 ' ' 1 1 11 1 2 12 2 ' ' 2 1 21 1 2 22 2 n n n n n n u c c u c c λ ξ λ ξ λ ξ λ ξ = + = + 1 1 0 2 'where are constants. ' c E u c ? a9a11a10 =a12 a13a14a11a15 a16a7a16a17 Slide 21 Alternative View Alternatively, one can view the solution as: 0 1 1 1 2 0 2 2 n n n n U U U Uλ λ a18 a19 a18 a19 a18 a19 =a20 a21 a20 a21a22 a23 a22 a23 a22 a23 0 1where n nU U U E u?= Λ = a24 a24 a24 a24 EXAMPLE 1 Comparison Comparing the solution of the semi-discretized problem where time is kept continuous [ ] 1 2 1 11 12 1 2 2 21 22 t t u ec c u e λ λ ξ ξ ξ ξ a25 a26a25a27a26 a25 a26 = a28 a29a28 a29 a28 a29a30a27a31 a30 a31 a30 a31 to the solution where time is discretized [ ]1 11 12 11 2 2 21 22 2 ' ' n n n u c c u ξ ξ λ ξ ξ λ a32 a33 a32a27a33 a32 a33 = a34 a35 a34 a35 a34 a35 a36a27a37 a36 a37 a34 a35 a36 a37 1 1 1 1 0 , ....n A E E u E E E E E E u ? ? ? ? = Λ = Λ ? Λ ? ? Λ ?a38a39a38a40 a40 A A A 16.920J/SMA 5212 Numerical Methods for PDEs 15 difference equation where time is continuous has exponential solution The .teλ The difference equation where time is discretized has power solution .nλ Slide 22 EXAMPLE 1 Comparison In equivalence, the transient solution of the difference equation must decay with time, i.e. 1nλ < for this particular form of time discretization. Slide 23 EXAMPLE 2 Leapfrog Time Discretization Consider a typical modal equation of the form t j du u ae dt μλ a0 a1 = +a2 a3a4 a5 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 1 1 where 2 n ndu u u h tdt h + ?? = = ? Substituting into the modal equation yields 1 1 2 n nu u h + ?? ( )t t nh u aeμλ = = + n hnu aeμλ= + Slide 24 16.920J/SMA 5212 Numerical Methods for PDEs 16 Reminder Recall that we are considering a typical modal equation which had been obtained from the original equation du Au b dt = + a0 a0a1a0 EXAMPLE 2 Leapfrog Time Discretization: Time Shift Operator ( )1 1 1 1 2 22n n n hn n n n hnu u u ae u h u u ha eh μ μλ λ+ ? + ?? = + a2 ? ? = Solution of u consists of the complementary solution nc , and the particular solution np , i.e. n n nu c p= + There are several ways of solving for the complementary and particular solutions. One way is through use of the shift operator S and characteristic polynomial. The time shift operator S operates on nc such that 1n nSc c += ( )2 1 2n n n nS c S Sc Sc c+ += = = Slide 25 EXAMPLE 2 Leapfrog Time Discretization: Time Shift Operator The complementary solution nc satisfies the homogenous equation 1 12 0 2 0 n n n n n n c h c c cSc h c S λ λ + ?? ? = ? ? = 16.920J/SMA 5212 Numerical Methods for PDEs 17 2 2 1( 2 ) 0 ( 2 1) 0 n n n n S c h Sc c S cS h S S λ λ ? ? = ? ? = Slide 26 EXAMPLE 2 Leapfrog Time Discretization: Time Shift Operator The solution to the characteristic polynomial is 2 2( ) 1h S h hσ λ λ λ= = ± + The complementary solution to the modal equation would then be 1 1 2 2n nnc β σ β σ= + The particular solution to the modal equation is 2 2 2 1 hn h n h h ahe ep e h e μ μ μ μλ= ? ? . Combining the two components of the solution together, nu ( ) ( )n nc p= + ( ) ( )2 2 2 21 2 2 21 1 2 1 hn hn n h h ahe eh h h h e h e μ μ μ μβ λ λ β λ λ λ a0 a1 a0 a1 = + + + ? + + a2 a3a2 a3 ? ?a4 a5 a4 a5 Slide 27 EXAMPLE 2 Leapfrog Time Discretization: Stability Criterion For the solution to be stable, the transient (complementary) solution must not be allowed to grow indefinitely with time, thus implying that ( ) ( ) 2 2 1 2 2 2 1 1 1 1 h h h h σ λ λ σ λ λ = + + < = ? + < characteristic polynomial 2( ) ( 2 1) 0p S S h Sλ= ? ? = σ1 and σ2 are the two roots 16.920J/SMA 5212 Numerical Methods for PDEs 18 is the stability criterion for the leapfrog time discretization scheme used above. Slide 28 EXAMPLE 2 Leapfrog Time Discretization: Stability Diagram The stability diagram for the leapfrog (or any general) time discretization scheme in the σ-plane is Slide 29 Stability Diagram in the λh-plane Alternatively, we can express the stability criterion for the leapfrog time discretization scheme as 1 1 s.t. 12hλ σ σσ a0 a1 = ? <a2 a3a4 a5 Since 1 and exp( )iσ σ θ< = , sinh iλ θ= for stability. The stability diagram for the leapfrog time discretization scheme in the λh-plane would therefore be as shown: Im(σ ) Re(σ ) -1 1 Region of Stability Re(λh) Im(λh) -1 1 Region of Stability 16.920J/SMA 5212 Numerical Methods for PDEs 19 EXAMPLE 2 Leapfrog Time Discretization In particular, by applying to the 1-D Parabolic PDE 2 2 u u t xυ ? ?= ? ? the central difference scheme for spatial discretization, we obtain which is the tridiagonal matrix Slide 30 EXAMPLE 2 Leapfrog Time Discretization According to analysis of a general triadiagonal matrix B(a,b,c), the eigenvalues of B are 2 2 cos , 1,..., 1 2 2cos j j jb ac j N N j N x piλ pi υλ a0 a1 = + = ?a2 a3a4 a5 a6 a7 a0 a1 = ? + a2 a3a8 a9 ?a4 a5a10 a11 The most “dangerous” mode is that associated with the eigenvalue of largest magnitude max 24xυλ = ? ? i.e. ( )( ) 2 2max 1 max max 2 2max 2 max max 1 1 h h h h h h σ λ λ λ σ λ λ λ = + + = ? + which can be plotted in the absolute stability diagram. 2 2 1 1 2 1 1 1 2 A xυ ?a12 a13 a14 a15 ?a14 a15 a14 a15 = a14 a15? a14 a15 a14 a15 a14 a15 ?a14 a15a16 a17 0 0 16.920J/SMA 5212 Numerical Methods for PDEs 20 One may note that jλ is always real and negative, thereby satisfying the criterion for stability of the space discretization of a parabolic PDE, keeping time continuous. Slide 31 EXAMPLE 2 Leapfrog Time Discretization: Absolute Stability Diagram for σ As applied to the 1-D Parabolic PDE, the absolute stability diagram for σ is In this case, 1σ and 2σ start out being on the unit circle (h = ?t = 0). However, the spurious root (refer to following slide) leaves the unit circle as h starts increasing. Therefore, the spurious root causes the leapfrog time discretization scheme to be unstable, irrespective of how small h = ?t is, although it does not affect the accuracy. The leapfrog time discretization for the 1-D Diffusion Equation is unstable. Slide 32 Im(σ ) Re(σ ) Unit circle σ1 with h increasing σ2 with h increasing Region of Instability Region of Stability σ2 at h = ?t = 0 σ1 at h = ?t = 0 16.920J/SMA 5212 Numerical Methods for PDEs 21 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 h λ λ σ λ → → = but not the accuracy of the scheme.) Slide 33 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 he h λ σ λ= + ( ) ( ) 1 2 2 2 2 4 42 2 2 1 2 2 1 1. 1 2 21 . 2 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 hh he h h λ λ λ λ λ λσ λ λλ λ ? + = + + + = + + + = + + + Slide 34 16.920J/SMA 5212 Numerical Methods for PDEs 22 EXAMPLE 3 Euler-Forward Time Discretization: Stability Analysis Analyze the stability of the explicit Euler-forward time discretization 1n ndu u u dt t + ? = ? as applied to the modal equation du udt λ= 1 1 Substituting where into the modal equation, we obtain (1 ) 0 n n n n duu u h h t dt u h uλ + + = + = ? ? + = Slide 35 EXAMPLE 3 Euler-Forward Time Discretization: Stability Analysis Making use of the shift operator S 1 (1 ) (1 ) [ (1 )] 0n n n n nc h c Sc h c S h cλ λ λ+ ? + = ? + = ? + = Therefore ( ) 1 and n n h h c σ λ λ βσ = + = The Euler-forward time discretization scheme is stable if 1 1hσ λ≡ + < or bounded by 1 s.t. 1 in the -plane.h hλ σ σ λ= ? < Slide 36 characteristic polynomial 16.920J/SMA 5212 Numerical Methods for PDEs 23 EXAMPLE 3 Euler-Forward Time Discretization: Stability Diagram The stability diagram for the Euler-forward time discretization in the λh-plane is Slide 37 EXAMPLE 3 Euler-Forward Time Discretization: Absolute Stability Diagram As applied to the 1-D Parabolic PDE, max 24xυλ λ ?= = ? max 2The stability limit for largest h t λ ?≡ ? = σ leaves the unit circle at σ = ?1, i.e. σ = 1+ λh = ?1 max 2h 2 hλ λ ?= ? a0 = since it is the extreme. Slide 38 Im(λh) Re(λh) -2 0 Region of Stability -1 Unit Circle Im(σ ) -1 1 Re(σ ) σ at h = ?t = 0 σ leaves the unit circle at λh = ?2 σ with h increasing 16.920J/SMA 5212 Numerical Methods for PDEs 24 Predictor-Corrector Time Discretization Consider the numerical stability of the following predictor-corrector time discretization scheme 1 1 1 1 1 ? ?1 ? 2 n n n n n n n n n duu u h dt duu u u h dt + + + + + = + a0 a1 = + +a2 a3 a4 a5 as applied to the typical modal equation taeudtdu μλ += of the parabolic PDE. Substituting dtdu and dtud ? into the predictor-corrector scheme yields ( ) ( ) 1 1 1 1 ( 1) ? where 1 ? ? 2 n n n hn n n n n h n u u h u ae t n t nh u u u h u ae μ μ λ λ + + + + + = + + = ? = a6 a7 = + + +a8 a9 Utilizing the shift operator 1 1? ? n n n n Su u Su u + + = = and rearranging the equations into matrix form, we obtain ( ) ( ) 1 ? 11 11 22 2 n hn n hS h u aehSuh S S μ λ λ a10 a11 ? + a10 a11a10a12a11a13 a14 a13 a14 =a13 a14a13 a14 a13 a14? + ? a15a12a16 a13 a14a13 a14 a15 a16a15 a16 To determine the characteristic polynomial, set ( ) ( ) ( ) ( ) 1 01 11 2 2 S h S h S S λ σ λ ? + Ρ = Ρ = =? + ? 16.920J/SMA 5212 Numerical Methods for PDEs 25 ( ) ( ) 2 2 2 2 11 0 2 0 (trivial root) 11 2 S S S h h h h σ λ λ σ σ λ λ a0 a1 Ρ = Ρ = ? ? ? =a2 a3a4 a5 a6 = = + + i.e. the scheme is a one-root method. Compared to 2 211 ....2he h hλ λ λ= + + + the scheme is second-order accurate. To obtain the particular solution, one can perform a matrix inversion and obtain ( ) 22 2 11 121 hhe heahe p h hhn n λλ λ μ μμ ??? ++ = with the complementary solution being n nn hhc a7a8 a9 a10a11a12 ++== 22 2 11 λλββσ The absolute stability diagram (showing 24x??= υλ ) for the 1-D Parabolic PDE is Im(σ ) Re(σ ) -1 1 0.5 Region of Stability Region of Instability h increasing from 0 h increasing further 16.920J/SMA 5212 Numerical Methods for PDEs 26 When h increases from zero, σ decreases from 1.0. As h continues to increase, σ reaches a minimum of 0.5 with λh = ?1 and then increases. As h increases further, σ returns to 1.0 with λh = ?2. Prior to this point, the scheme is stable. Increasing h and thus σ beyond this point renders the scheme unstable. Hence, this predictor-corrector scheme is stable for small h’s and unstable for large h’s; the limit for stability is λh = ?2 (from above). In general, we can analyze the absolute stability diagram for the predictor-corrector time discretization method in terms of 2( ) : ( ) 1 2hh h λσ σ λ λ= + + or : 1 2 1h hλ λ σ= ? ± ? λ, the eigenvalue(s) of the A matrix can take on complex forms depending on the governing equation (as opposed to negative real values for the 1-D parabolic PDE with central differencing for the spatial derivative). RELATIONSHIP BETWEEN σ AND λh σ = σ(λh) Thus far, we have obtained the stability criterion of the time discretization scheme using a typical modal equation. We can generalize the relationship between σ and λh as follows: ? Starting from the set of coupled ODEs du Au bdt = + a0 a0a1a0 ? Apply a specific time discretization scheme like the leapfrog time discretization as in Example 2 1 1 2 n ndu u u dt h + ?? = Slide 39 16.920J/SMA 5212 Numerical Methods for PDEs 27 RELATIONSHIP BETWEEN σ AND λh σ = σ(λh) ? The above set of ODEs becomes 1 1 2 n n n nu u Au b h + ?? = + a0 a0 a0 a0 ? Introducing the time shift operator S 1 2 2 2 n n n n nn uSu hAu hb S S SA I u b h ? = + + a1 a2 ?? = ?a3 a4 a5 a6 a7 a7 a7 a7 a7 a7 ? 1 1 1 1 1 2 nS SE AE E E E u E b h ? ? ? ? ? a8 a9 ?? = ?a10 a11 a12 a13 a14 a14 Slide 40 RELATIONSHIP BETWEEN σ AND λh σ = σ(λh) ? Putting 1 1, nn n nU E u F E b? ?= = a15a15 a15 a15 we obtain 1 1 2 n nS SE E U F h ? ? a16 a17 ?Λ ? = ?a18 a19 a20 a21 a22 a22 i.e. 1 2 n nS S U F h ? a23 a24 ?Λ ? = ?a25 a26 a27 a28 a29 a29 which is a set of uncoupled equations. 1 1 Premultiplying on the LHS and RHS and introducing operating on n E I EE u ? ?= a30 Λ 1 2 S S h ?? 16.920J/SMA 5212 Numerical Methods for PDEs 28 Hence for each j, j = 1,2,….,N?1, 1 2j j j S S U F hλ ? a0 a1 ?? = ?a2 a3 a4 a5 Slide 41 RELATIONSHIP BETWEEN σ AND λh σ = σ(λh) Note that the analysis performed above is identical to the analysis carried out using the modal equation j dU U F dt λ a6 a7 = +a8 a9a10 a11 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 n j j j U U u EU= a12 a12 Slide 42 RELATIONSHIP BETWEEN σ AND λh σ = σ(λh) Hence, applying any “consistent” numerical technique to each equation in the set of coupled linear ODEs is mathematically equivalent to 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 Slide 43 16.920J/SMA 5212 Numerical Methods for PDEs 29 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. However, use of implicit schemes requires 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. Slide 44 IMPLICIT TIME-MARCHING SCHEME Euler-Backward Consider the Euler-backward scheme for time discretization 1 1n n ndu u u dt h + + ?a0 a1 =a2 a3a4 a5 Applying the above to the modal equation for parabolic PDE tdu u aedt μλ= + yields ( ) ( ) ( ) 1 11 111 n n n hn n hn n u u u ae h h u u ahe μ μ λ λ + ++ ++ ? a6 a7= + a8 a9 ? ? = Slide 45 IMPLICIT TIME-MARCHING SCHEME Euler-Backward Applying the S operator, ( ) ( )11 1 n hnh S u aheμλ + a10 a11 ? ? =a12 a13 16.920J/SMA 5212 Numerical Methods for PDEs 30 the characteristic polynomial becomes ( ) ( ) ( )1 1 0S h Sσ λ a0 a1 Ρ = Ρ = ? ? =a2 a3 The principal root is therefore 2 21 1 ....1 h hhσ λ λλ= = + + +? 2 2 1which, upon comparison with 1 .... , is only 2 first-order accurate. he h hλ λ λ= + + + The solution is ( ) ( ) 11 1 1 1 n u h n h aheU h h e μ μβ λ λ +a4 a5 = +a6 a7? ? ?a8 a9 Slide 46 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. A small time-step size is used for the short time- scaled events, and then a large time-step size used for the longer time-scaled events. There is no constraint on hmax. Slide 47 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 16.920J/SMA 5212 Numerical Methods for PDEs 31 intensive/expensive compared to the multiplication/ addition operations of explicit schemes. Slide 48 SUMMARY ? Stability Analysis of Parabolic PDE a0 Uncoupling the set. a0 Integrating each equation in the uncoupled set → modal equation. a0 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. Slide 49 Reference: Numerical Computation of Internal and External Flows, Vol I & II by C. Hirsch, 1992, Wiley Series.