Numerical Solutions to Partial Differential Equations Zhiping Li - PowerPoint PPT Presentation

Numerical Solutions to Partial Differential Equations Zhiping Li LMAM and School of Mathematical Sciences Peking University

A Model Problem in a 2D Box Region Let us consider a model problem of parabolic equation: u t = a ( u xx + u yy ) , ( x , y ) ∈ Ω , t > 0 , u ( x , y , 0) = u 0 ( x , y , 0) , ( x , y ) ∈ Ω , u ( x , y , t ) = 0 , ( x , y ) ∈ ∂ Ω , t > 0 , where a > 0 is a constant, Ω = (0 , X ) × (0 , Y ) ⊂ R 2 . 1 For integers N x ≥ 1 and N y ≥ 1, let h x = △ x = XN − 1 and x h y = △ y = YN − 1 be the grid sizes in the x and y directions; y 2 a uniform parallelopiped grid with the set of grid nodes J Ω × R + = { ( x j , y k , t m ) : 0 ≤ j ≤ N x , 0 ≤ k ≤ N y , m ≥ 0 } , where x j = j h x , y k = k h y , t m = m τ ( τ > 0 time step size). 3 the space of grid functions U = { U m j , k = U ( x j , y k , t m ) : 0 ≤ j ≤ N x , 0 ≤ k ≤ N y , m ≥ 0 } .

Finite Difference Methods for Parabolic Equations Explicit and Implicit Schemes in Box Regions The Explicit Scheme in a 2D Box Region The Forward Explicit Scheme in a 2D Box Region The forward explicit scheme and its error equation U m +1 − U m � U m j +1 , k − 2 U m j , k + U m U m j , k +1 − 2 U m j , k + U m � j , k j , k j − 1 , k j , k − 1 = a + . h 2 h 2 τ x y e m +1 = [1 − 2( µ x + µ y )] e m � e m j +1 , k + e m � � e m j , k +1 + e m � − T m j , k + µ x + µ y j , k τ, j − 1 , k j , k − 1 j , k where µ x = a τ x , µ y = a τ y are the grid ratios in x and y directions. h 2 h 2 The truncation error is Tu ( x , y , t ) = 1 2 u tt ( x , y , t ) τ − a ∂ 4 x u ( x , y , t ) h 2 x + ∂ 4 y u ( x , y , t ) h 2 � � y 12 + O ( τ 2 + h 4 x + h 4 y ) . 3 / 35

Finite Difference Methods for Parabolic Equations Explicit and Implicit Schemes in Box Regions The Explicit Scheme in a 2D Box Region The Forward Explicit Scheme in a 2D Box Region 1 the condition for the maximum principle: µ x + µ y ≤ 1 2 ; 2 Fourier modes and amplification factors λ l ( l = ( l x , l y )): l e i ( α x x j + α y y k ) = λ m l e i ( l x j π N − 1 + l y k π N − 1 ) , U m j , k = λ m x y X , − N x + 1 ≤ l x ≤ N x , α y = l y π α x = l x π Y , − N y + 1 ≤ l y ≤ N y ; ( l x , l y represent the frequency or wave number in x , y direction.) � � + µ y sin 2 α y h y µ x sin 2 α x h x 3 amplification factor λ l = 1 − 4 2 2 � � 2 N x + µ y sin 2 l y π µ x sin 2 l x π = 1 − 4 ; 2 N x 4 L 2 stable if and only if µ x + µ y ≤ 1 / 2; 5 Convergence rate is O ( τ + h 2 x + h 2 y ). 4 / 35

The θ -Scheme in a 2D Box Region U m +1 − U m � � � � + δ 2 + δ 2 δ 2 δ 2 j , k j , k y y x U m x U m +1 = (1 − θ ) a j , k + θ a j , k τ h 2 h 2 h 2 h 2 x y x y � U m j +1 , k − 2 U m j , k + U m U m j , k +1 − 2 U m j , k + U m � j − 1 , k j , k − 1 = (1 − θ ) a + h 2 h 2 x y � U m +1 j +1 , k − 2 U m +1 + U m +1 U m +1 j , k +1 − 2 U m +1 + U m +1 � j − 1 , k j , k − 1 j , k j , k + θ a + . h 2 h 2 x y The error equation: � � � � �� [(1 + 2 θ )( µ x + µ y )] e m +1 e m +1 j +1 , k + e m +1 e m +1 j , k +1 + e m +1 = θ µ x + µ y j − 1 , k j , k − 1 j , k + [1 − 2(1 − θ )( µ x + µ y )] e m j , k m + 1 e m j +1 , k + e m e m j , k +1 + e m � � � � �� + (1 − θ ) µ x + µ y − τ T 2 , j − 1 , k j , k − 1 j , k

Finite Difference Methods for Parabolic Equations Explicit and Implicit Schemes in Box Regions The θ -Scheme in a 2D Box Region The θ -Scheme in a 2D Box Region The truncation error O ( τ 2 + h 2 � x + h 2 if θ = 1 y ) , 2 , T m + ∗ = j if θ � = 1 O ( τ + h 2 x + h 2 y ) , 2 , 1 the condition for the maximum principle: 2( µ x + µ y ) (1 − θ ) ≤ 1; Nx + ly k π i ( lx j π l e i ( α x x j + α y y k ) = λ m Ny ) , 2 Fourier modes: U m j , k = λ m l e X , α y = ly π Nx , α y y k = ly k π l = ( l x , l y ), α x = lx π Y , α x x j = lx j π Ny ; 3 amplification factor � � 2 N x + µ y sin 2 l y π µ x sin 2 l x π 1 − 4(1 − θ ) 2 N x λ l = ; � � 2 N x + µ y sin 2 l y π µ x sin 2 l x π 1 + 4 θ 2 N x 6 / 35

Finite Difference Methods for Parabolic Equations Explicit and Implicit Schemes in Box Regions The θ -Scheme in a 2D Box Region The θ -Scheme in a 2D Box Region 4 for θ ≥ 1 / 2, unconditionally L 2 stable; 5 for 0 ≤ θ < 1 / 2, L 2 stable iff 2(1 − 2 θ )( µ x + µ y ) ≤ 1; 6 the matrix of the linear system is still symmetric positive definite and diagonal dominant, however, each row has now up to 5 nonzero elements with a band width of the order O ( h − 1 ); 7 if solved by the Thompson method, the cost is O ( h − 1 ) times of that of the explicit scheme; 8 in 3D, µ x + µ y ⇒ µ x + µ y + µ z , and O ( h − 1 ) ⇒ O ( h − 2 ), if solved by the Thompson method, the cost is O ( h − 2 ) times of that of the explicit scheme. 7 / 35

Finite Difference Methods for Parabolic Equations 2D and 3D ADI and LOD Schemes Alternative Approaches for Solving n -D Parabolic Equations To reduce the computational cost, we may consider to apply highly efficient iterative methods to solve the linear algebraic equations, for example, the preconditioned conjugate gradient method; the multi-grid method; etc.. Alternatively, to avoid the shortcoming of the implicit difference schemes for high space dimensions, we may develop the alternating direction implicit (ADI) schemes; the locally one dimensional (LOD) schemes. 8 / 35

Finite Difference Methods for Parabolic Equations 2D and 3D ADI and LOD Schemes A 2D Alternating Direction Implicit (ADI) Scheme A Fractional Steps 2D ADI Scheme by Peaceman and Rachford 1 in the odd fractional steps implicit in x and explicit in y , and in even fractional steps implicit in y and explicit in x : � 1 − 1 � � 1 + 1 � m + 1 2 µ x δ 2 2 µ y δ 2 U m U 2 = j , k , x j , k y � 1 − 1 � � 1 + 1 � m + 1 2 µ y δ 2 U m +1 2 µ x δ 2 = U 2 , y j , k x j , k 2 numerical boundary conditions are easily imposed directly by m + 1 m + 1 those of the original problem, since U 2 ∼ u 2 ; j , k j , k 9 / 35

Finite Difference Methods for Parabolic Equations 2D and 3D ADI and LOD Schemes A 2D Alternating Direction Implicit (ADI) Scheme A Fractional Steps 2D ADI Scheme by Peaceman and Rachford 3 one step equivalent scheme � 1 − 1 � � 1 − 1 � � 1 + 1 � � 1 + 1 � 2 µ x δ 2 2 µ y δ 2 U m +1 2 µ x δ 2 2 µ y δ 2 U m = j , k . x y x y j , k 4 Crank-Nicolson scheme � 1 − 1 x − 1 � � 1 + 1 x + 1 � 2 µ x δ 2 2 µ y δ 2 U m +1 2 µ x δ 2 2 µ y δ 2 U m = j , k . y j , k y m + 1 m + 1 = a 2 τ 3 [ u xxyyt ] + O ( τ 5 + τ 3 ( h 2 5 Since µ x µ y δ 2 x δ 2 x + h 2 y )) , y δ t u 2 2 j , k j , k the truncation error of the scheme is O ( τ 2 + h 2 x + h 2 y ). 10 / 35

Finite Difference Methods for Parabolic Equations 2D and 3D ADI and LOD Schemes A 2D Alternating Direction Implicit (ADI) Scheme Stability of the 2D ADI Scheme by Peaceman and Rachford l e i ( l x j π N − 1 + l y k π N − 1 1 For the Fourier mode U m j , k = λ m ) , x y 2 N x )(1 − 2 µ y sin 2 l y π (1 − 2 µ x sin 2 l x π 2 N y ) λ l = , (1 + 2 µ x sin 2 l x π 2 N x )(1 + 2 µ y sin 2 l y π 2 N y ) the scheme is unconditionally L 2 stable; 2 the two fractional steps can be equivalently written as � � m + 1 m + 1 m + 1 j , k + µ y + µ x = (1 − µ y ) U m � U m j , k − 1 + U m � (1 + µ x ) U 2 j − 1 , k + U 2 2 U , j , k +1 j , k j +1 , k 2 2 � � m + 1 m + 1 m + 1 + µ x + µ y � � (1 + µ y ) U m +1 U m +1 j , k − 1 + U m +1 = (1 − µ x ) U 2 U j − 1 , k + U 2 2 , j , k j , k +1 j , k j +1 , k 2 2 thus, the maximum principle holds if max { µ x , µ y } ≤ 1; 11 / 35

Finite Difference Methods for Parabolic Equations 2D and 3D ADI and LOD Schemes A 2D Alternating Direction Implicit (ADI) Scheme Cost of the 2D ADI Scheme by Peaceman and Rachford 1 order the linear system by ( N x − 1) k + j and ( N y − 1) j + k in odd and even steps, the corresponding matrixes are tridiagonal; 2 the computational cost: 3 times of that of the explicit scheme. 12 / 35

Finite Difference Methods for Parabolic Equations 2D and 3D ADI and LOD Schemes A 2D Alternating Direction Implicit (ADI) Scheme The Idea of Peaceman and Rachford Doesn’t Work for 3D In 3D, we can still construct fractional step scheme by 1 dividing each time step into 3 fractional steps; 2 introducing U m + 1 3 and U m + 2 3 at t m + 1 3 and t m + 2 3 ; 3 in the 3 fractional time steps, applying schemes which are in turn implicit in x -, y - and z -direction and explicit in the other 2 directions respectively. However, the equivalent one step scheme is, up to a higher order term, the same as the θ -scheme with θ = 1 / 3, which has a local truncation error O ( τ + h 2 x + h 2 y ) instead of what we expect to have O ( τ 2 + h 2 x + h 2 y ) for an ADI method. 13 / 35