16920J/SMA5212 Numerical Methods for Partial Differential equations Lecture 5 Finite Differences Parabolic problems B.C. Khe Thanks to franklin tan 19 February 2003
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 OUTLINE Governing equation Stability Analysis 3 Examples Relationship between oand /h Implicit Time-Marching Scheme Summary Side 2 GOVERNING EQUATION Consider the parabolic pde in 1-D du a2u ∈[0.x] subject to u=lo atx=0, u=u atx=Tt Ifu≡ viscosity→ Diffusion Equation If u=thermal conductivity Heat Conduction Equation Side 3 STABILITY ANALYSIS Discretization Keeping time continuous, we carry out a spatial discretization of the rhs of du a at dx 0 X=丌 2
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 If υ ≡ viscosity → Diffusion Equation 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 x t x υ π ∂ ∂ = ∈ ∂ ∂ 0 subject to u u at x 0, u u at x = = = π = π x = 0 x = π 0 u uπ u ( x,t) = ? 2 2 u u t x υ ∂ ∂ = ∂ ∂ x = 0 x = π 0 x 1 x 2 x N 1 x − N x
16.920J/SMA 5212 Numerical Methods for PDEs There is a total of N+I grid points such that x=jAx, j=0,1,2 Side 4 STABILITY ANALYSIS Discretization Use the Central Difference scheme for a-u which is second-order accurate Schemes of other orders of accuracy may be constructed 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 The order of accuracy, p of a spatial difference scheme is represented as O(ArP) Generally, to represent the spatial operator to a higher order of accuracy, more nodal j2-1 户+1j+2 Consider the following procedure of determining the spatial operator dx cyo(△x2)
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 ( ) p O ∆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 ✁ ✂ ✄ ☎ ✆ up to the order of accuracy ( ) 2 O ∆x : There is a total of 1 grid points such that , 0,1,2,...., N j x j x j N + = ∆ = 2 U 2 se the Central Difference Scheme for u x ∂ ∂ 2 1 1 2 2 2 2 ( ) j j j j u u u u O x x x + − + − ✝ ✞ ∂ = + ∆ ✟ ✠ ∂ ∆ ✡ ☛ j−2 • j−1 • j • j+1 • j+2 • j du dx ☞ ✌ ✍ ✎ ✏ ✑
16.920J/SMA 5212 Numerical Methods for PDEs Let d be represented by u at the nodes j-1, J, and j*l with a-1, d and a being the coefficients to be determined ie 11-+a0u+a1u O(△xP d Seek Taylor Expansions for u;-,u, and u +r about u, and present them in a (Note that p is not known a priori but is determined at the end of the analysis when the as are made known. This column consists of all the terms on the LHS of (1) 0 a1|-△xa △x3a 0 0 a1 x·C1 △x-a1 a S S2 Each cell in this row comprises the sum of its corresponding column
16.920J/SMA 5212 Numerical Methods for PDEs 4 1. Let j du dx ✁ ✂ ✄ ☎ ✆ 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 ( ) p j j j j du u u u O x dx α− − α α + ✝ ✞ + + + = ∆ ✟ ✠ ✡ ☛ 2. Seek Taylor Expansions for j 1 u − , j u and j 1 u + about j u 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 ′′′ uj ′ 0 1 0 0 α−1uj−1 α−1 1 −∆x α− ⋅ 2 1 1 2 ∆x α− ⋅ 3 1 1 6 − ∆x α− ⋅ α0uj α0 0 0 0 1 j 1 α u + α1 1 ∆x ⋅α 2 1 1 2 ∆x ⋅α 3 1 1 6 ∆x ⋅α 1 1 k j k j k k u α u = + =− ′ + ☞ S1 S2 S3 S4 ( 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 where S=(a1+a+a)u S2=(1-△xa1+△xa1) +-△ S4=-△x2a1+△xa1| S1+S2+S3+S4+ 3. Make as many S,'s as possible vanish by choosing appropriate as In this instance, since we have three unknowns a_, a and a, we can therefore set (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 s pertaining to the lower-order terms to zero, thereafter followed by ncreasingly higher-order terms) 0 0 0 Solving the 2△
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 = + =− ′ + = + + + + 3. Make as many Si ’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 i S pertaining to the lower-order terms to zero, thereafter followed by increasingly higher-order terms.) Hence, 1 0 1 0 1 1 1 1 1 0 1 1 0 1 0 x α α α − ✁ ✂ ✁ ✂✄✁ ✂ ☎ ✆ ☎ ✆ ☎ ✆ ☎ ✆ − = − ☎ ✆ ☎ ✆ ∆ ☎ ✆ ☎ ✆ ☎ ✆ ✝ ✞✄✝ ✞ ☎ ✆ ✝ ✞ 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 α α α α α α α α α − − − − = + + ′ = −∆ ⋅ + ∆ ⋅ ✟ ✠ ′′ = ∆ ⋅ + ∆ ⋅ ✡ ☛ ☞ ✌ ✟ ✠ ′′′ = − ∆ ⋅ + ∆ ⋅ ✡ ☛ ☞ ✌
16.920J/SMA 5212 Numerical Methods for PDEs 4. Substituting the ax's into u, +2au +=S+S2+S,+S+.yields 2△x(1“ △x2·l,+ higher- order terms In other words, l--+0(△x2)+ dx 0(△x2),p i.e. the above representation is accurate up to O(4r) Some useful points to note These 4 steps are the general procedure used to obtain the representation of the spatial operator up to the order of accuracy O(ArP) d2u For other spatial operators, say we simply replace/du in(1)with dx d the said spatial operator 3. For one-sided representations, one can choose nodal points u,k,k>0. This may be important especially for representations on a boundary. For example +a+am+a2+…=0(△) dt One possibility is dx 2Ax 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) 6
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 = + =− ′ + = + + + + yields ( ) 2 1 1 1 1 2 6 j j j j u u u x u x + − ′ ′′′ − − = ∆ ⋅ + ∆ higher-order terms In other words, ( ) 1 1 2 .... 2 j j j j du u u u O x dx x + − − ✁ ✂ ′ = = + ∆ + ✄ ☎ ∆ ✆ ✝ i.e. the above representation is accurate up to ( ) 2 O ∆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 ( ) p O ∆x . 2. For other spatial operators, say 2 2 j d u dx ✞ ✟ ✠ ✡ ☛ ☞ , we simply replace j du dx ✌ ✍ ✎ ✏ ✑ ✒ in (1) with the said spatial operator. 3. For one-sided representations, one can choose nodal points uj+k , k ≥ 0 . This may be important especially for representations on a boundary. For example 0 1 1 2 2 .... ( ) p j j j j du u u u O x dx α α + α + ✓ ✔ + + + + = ∆ ✕ ✖ ✗ ✘ One possibility is ( ) 1 2 2 3 4 2 j j j j du u u u O x dx x − + + + ✙ ✚ + = ∆ ✛ ✜ ∆ ✢ ✣ 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). ( ), 2 p O ∆x p =
16.920J/SMA 5212 Numerical Methods for PDEs STABILITY ANALYSIS Discretization We obtain at x dt Ar(uo-2u+"2) dax2(4-2n2+n) 动=a2(1-21+1 N-1 (a2-2ux=1+lx) dt Note that we need not evaluate u at x=x, and x=XN since u, and u, are given as boundary conditions Side 6 STABILITY ANALYSIS Matrix Formulation Assembling the system of equations, we obtain dt 0 dhu,△x2 dt 0 1-2 DUN Slide 7 7
16.920J/SMA 5212 Numerical Methods for PDEs 7 STABILITY ANALYSIS Discretization We obtain at 1 1 2 1 2 : ( 2 ) o du x u u u dt x υ = − + ∆ 2 2 2 1 2 3 : ( 2 ) du x u u u dt x υ = − + ∆ 2 1 1 : ( 2 ) j j j j j du x u u u dt x υ = − − + + ∆ 1 1 2 2 1 : ( 2 ) N N N N N du x u u u dt 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 0 1 2 1 0 1 2 1 1 0 1 2 o j j N N N du u u dt x du u dt du x u dt u du u x dt υ υ υ − − ✁ ✁ − ✂ ✄ ✁ ✁ ✂ ✄ ∆ ✂ ✄ ✂ ✄ ✂ ✄ ✂ ✄ ✂ ✄ ✂ ✄ ✂ ✄ ✂ ✄ ✂ ✄ ✂ ✄ ✂ ✄ − ✂ ✄ ✂ ✄ ✂ ✄ ✂ ✄ ✂ ✄ ✂ ✄ ✂ ✄ ✂ ✄ ✂ ✄ ✂ ✄ ✂ ✄ ✂ ✄ ✂ ✄ = + ✂ ✄ ✂ ✄ ✂ ✄ ✂ ✄ ∆ − ✂ ✄ ✂ ✄ ✂ ✄ ✂ ✄ ✂ ✄ ✂ ✄ ✂ ✄ ✂ ✄ ✂ ✄ ✂ ✄ ✂ ✄ ✂ ✄ ✂ ✄ ✂ ✄ ✂ ✄ ✂ ✄ ✂ ✄ ✂ ✄ ✂ ✄ ✂ ✂ ✄ ✂ ✄ ✂ ✄ ✂ ✂ ✄ ✂ ✄ − ✂ ✄ ☎ ✆✝☎ ✆ ✂ ✂ ∆☎ ✆ ✂☎ ✆✄ ✄ ✄ ✄ 0 ✄ 0 A
16.920J/SMA 5212 Numerical Methods for PDEs STABILITY ANALYSIS PDE to Coupled ODEs Or in compact form du lu+b where u=[a1u2…… △x2 We have reduced the I-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 元={4,2,…,…,不- A-11=0 For each eigenvalue 1. we can evaluate the eigenvector p consisting of a set of mesh point values v/,i.e y=Y vy STABILITY ANALYSIS Eigenvalue and Eigenvector of Matrix A The(N-1)x(N-1)matrix E formed by the(N-1)columns r diagonalizes the matrix A by E-AE=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 ,....,λ j ,....,λN −1} from det( A− λI ) = 0. 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 by j N N E N V A − × − − 1 E AE − = Λ where [ 1 2 1 ] T N u u u u = − 2 2 0 0 0 T o N u u b x x υ υ ✁ ✂ = ✄ ☎ ∆ ∆ ✆ ✝ ✞ du Au b dt = + ✟ ✟✠✟ 1 2 1 Tj j j j V N v v v − ✡ ☛ = ☞ ✌
16.920J/SMA 5212 Numerical Methods for PDEs 0 0 Slide 10 STABILITY ANALYSIS Coupled oDEs to Uncoupled ODEs Starting fre Au+b Premultiplication by E yield e-ld e-l4u+E-b AEE)b+E-b E=(EAEE-b+E-b STABILITY ANALYSIS Coupled oDEs to Uncoupled ODEs Continuing from di AE-u+E-6 Letu=E-i and F=E-B,we have
16.920J/SMA 5212 Numerical Methods for PDEs 9 Slide 10 STABILITY ANALYSIS Coupled ODEs to Uncoupled ODEs Starting from du Au b dt = + ✟ ✟✠✟ 1 Premultiplication by E yields − 1 du 1 1 E E Au E b dt − − − = + Slide 11 STABILITY ANALYSIS Coupled ODEs to Uncoupled ODEs Continuing from 1 du 1 1 E E u E b dt − − − = Λ + 1 1 Let U E u and F E b, we have − − = = ✁ ✁ ✁ ✁ 1 2 1 where N λ λ λ − ✂ ✄ ☎ ✆ ☎ ✆ ☎ ✆ Λ = ☎ ✆ ☎ ✆ ☎ ✆ ✝ ✞ 0 0 ( ) 1 du 1 1 1 E E A EE u E b dt − − − − = + ✟ ✟ ✟ I Λ ( ) 1 du 1 1 1 E E AE E u E b dt − − − − = + ✠ ✠ ✠
16.920J/SMA 5212 Numerical Methods for PDEs di=au+F which is a set of Uncoupled ODEs STABILITY ANALYSIS Coupled oDEs to Uncoupled ODEs Expanding yields dU=U1+F due=au, +F2 d=4+F Since the equations are independent of one another, they be solved The idea then is to solve for u and determine i= eU STABILITY ANALYSIS Coupled oDEs to Uncoupled ODEs Considering the case of b independent of time, for the general j"equation is the solution forj=1, 2, . N-l
16.920J/SMA 5212 Numerical Methods for PDEs 10 d U U F dt = Λ + ✁ ✁ ✁ which is a set of Uncoupled ODEs! Slide 12 STABILITY ANALYSIS Coupled ODEs to Uncoupled ODEs Expanding yields 1 1 1 1 dU U F dt = λ + 2 2 2 2 dU U F dt = λ + j j j j dU U F dt = λ + 1 1 1 1 N N N N dU U F dt λ − = − − + − Since the equations are independent of one another, they can be solved separately. The idea then is to solve for U and determine u = EU ✂ ✂ ✂ Slide 13 STABILITY ANALYSIS Coupled ODEs to Uncoupled ODEs Considering the case of independent of time, for the general equation, th b j ✄ jt 1 j j j j U c e F λ λ = − is the solution for j = 1,2,….,N−1