Tut8_Sol

MATH3240 Numerical Methods for Differential Equations Tutorial 8 (Mar. 30, 2023) 1 Recall: 1. Numerical methods for stiff ODEs Consider the system of ODEs X 0 ( t ) = AX ( t ) + W ( t ) , (1) where X ( t ) = ( x 1 ( t ); x 2 ( t ); . . . ; x n ( t )) T is the solution vector, W(t) is a given vector, and A is a n × n matrix. Assume that A has n linearly independent eigenvectors { v i } n i =1 associated with n distinct eigenvectors λ j , j = 1 , 2 , ..., n and ψ ( t ) is a particular solution to the equation (1). Then the general solution X(t) of the ODE (1) can be represented by X ( t ) = n X i =1 α i e λ i t v i + ψ ( t ) (2) where α 1 , α 2 , . . . , α n are n constants. Consider only the case with asymptotical stability and assume that all eigenvalues satisfy Re ( λ i ) < 0 , i = 1 , 2 , . . . , n. The solution part = ∑ n i =1 α i e λ i t v i is called the transient solution . When t → ∞ , the solution X ( t ) tends to the particular solution ψ ( t ) . And we call ψ ( t ) as the steady-state solution. Assume that σ ≤ Re ( λ i ) ≤ τ < 0 , i = 1 , 2 , . . . , n. Then we call the factor r s = σ τ as the stiffness quotient . And a linear system of ODEs with constant coefficients is stiff if all the eigenvalues of A have negative real parts and r s >> 1 . Definition 1. A-stablility A multistep method is said to be A-stable if all the roots of the polynomial φ ( z ) lie strictly inside the unit disk | z | < 1 when Re ( α ) < 0 . Theorem 1. For a linear multistep method, if it is A-stable, then it must be an implicit method, and its order can not exceed 2. 2. Absolute stability and region of absolute stability Consider a general model ODE: ( x 0 ( t ) = αx ( t ) , x (0) = 1 (3) where α is a complex number. A multistep method for approximating (3) is said to be absolute stable if it holds | x n | → 0 as t n → ∞ . Clearly, the solution x n depends on step size h and parameter α . We call the region D = { any complex z = hα such that | x n | → 0 as t n → ∞} as the region of absolute stability of the multistep method. 1

3. Numerical methods for boundary value problems Consider another type of ODE system: ( x 00 ( t ) = f ( t, x, x 0 ) , a < t < b x ( a ) = α, x ( b ) = β (4) which is called a second order boundary value problem (BVP). Theorem 2. Assume that ∂f ∂x is continuous , nonnegative and bounded in the domain Ω = { ( t, x ); 0 ≤ t ≤ 1 , -∞ < x < ∞} , then the two-point BVP ( x 00 ( t ) = f ( t, x, x 0 ) , 0 < t < 1 x (0) = 0 , x (1) = 0 (5) has a unique solution. 2

Next we try to find eigenvectors by solving ( A - λ i I ) v i = 0 , i = 1 , 2 . We obtain two eiqenvectors v 1 = 1 - 1 , v 2 = 1 1 . Letting P = ( v 1 , v 2 ) = 1 1 - 1 1 , then we have P - 1 = 1 2 1 - 1 1 1 , AP = PD, D = α - β 0 0 α + β , and P - 1 X 0 ( t ) = P - 1 AX ( t ) = P - 1 APY ( t ) = DY ( t ) where Y ( t ) = P - 1 X ( t ) = y 1 ( t ) y 2 ( t ) . Then, we have Y 0 ( t ) = DY ( t ) , or equivalently ( y 0 1 ( t ) = ( α - β ) y 1 ( t ) , y 2 ( t ) 0 = ( α + β ) y 2 ( t ) , Y (0) = P - 1 X (0) = 1 1 . Solving it, we get y 1 = e ( α - β ) t , y 2 = e ( α + β ) t . Hence we obtain the solution to original ODE X ( t ) = PY ( t ) = e ( α - β ) t + e ( α + β ) t - e ( α - β ) t + e ( α + β ) t . Note when α < β < 0 , both x 1 ( t ) , x 2 ( t ) decay to zero. (b) We use the explicit Euler scheme to compute an approximation solution to original ODE X n +1 = X n + hAX n = ( I + hA ) X n , n = 0 , 1 · · · with X 0 = 2 0 . Note X n +1 = X n + hAX n = ( I + hA ) X n = ( PP - 1 + hPDP - 1 ) X n = P ( I + hD ) P - 1 X n Now taking the transformation Z n = P - 1 X n , n = 0 , 1 · · · , we have Z n +1 = ( I + hD ) Z n , where Z 0 = P - 1 X 0 = (1 , 1) T , which gives us Z n = ( I + hD ) n Z 0 = ( 1 + h ( α - β ) ) n ( 1 + h ( α + β ) ) n Hence X n = PZ n = ( 1 + h ( α - β ) ) n + ( 1 + h ( α + β ) ) n - ( 1 + h ( α - β ) ) n + ( 1 + h ( α + β ) ) n . To ensure X n → 0 as n → 0 , we need | 1 + h ( α - β ) | < 1 , | 1 + h ( α + β ) | < 1 , which is equivalent to 0 < h < - 2 α - β , 0 < h < - 2 α + β . Since - 2 α - β > - 2 α + β > 0 , we know X n → 0 as n → ∞ , if and only if 0 < h < - 2 α - β . 4

3. Consider the method x n +1 = x n + h ((1 - θ ) f n + θf n +1 ) , 0 ≤ θ ≤ 1 . Show that it is A-stable if and only if θ ≥ 1 / 2 . Solution. Let f ( t, x ) = αx , with Re ( α ) < 0 , then the numerical method can be written as x n +1 - x n = h ((1 - θ ) αx n + θαx n +1 ) . The characteristic polynomial is φ ( z ) = p ( z ) - hαq ( z ) = z - 1 - hα ( θz + (1 - θ )) , and its root is z = 1 + hα (1 - θ ) 1 - hαθ . By definition, it is A-stable iff (let α = α 1 + iα 2 ) | z | 2 = (1 + α 1 h (1 - θ )) 2 + ( α 2 h (1 - θ )) 2 (1 - α 1 hθ ) 2 + ( α 2 hθ ) 2 < 1 i.e. 2 α 1 h + (1 - 2 θ )( α 2 1 + α 2 2 ) h 2 < 0 . Clearly this inequality holds true for any h > 0 if θ ≥ 1 / 2 . Then the method is A-stable if θ ≥ 1 / 2 . Consider the case θ < 1 / 2 , then we have 2 α 1 h + (1 - 2 θ )( α 2 1 + α 2 2 ) h 2 > 0 i.e. | z | 2 > 1 when h > - 2 α 1 (1 - 2 θ )( α 2 1 + α 2 2 ) . Hence in the case θ < 1 / 2 , this leads to a contradiction with "A-stable". 4. Consider the following 2-step method 3 x n +2 - 4 x n +1 + x n = 2 hf n +2 . Prove that the region of absolute stability contains the negative real line. Solution. Let f ( t, x ) = αx , we have φ ( z ) = 3 z 2 - 4 z + 1 - 2 wz 2 with w = αh = w 1 + iw 2 . So the root is z = 2 ± √ 1 + 2 w 3 - 2 w . Note that the negative real line can be expressed as w 1 < 0 , w 2 = 0 , then we only need to verify this w = w 1 + iw 2 = w 1 < 0 satisfying | z | < 1 . (a) If 1 + 2 w 1 ≥ 0 , i.e. - 1 2 ≤ w 1 ≤ 0 , we have 0 ≤ 1 + 2 w 1 < 1 , 3 < 3 - 2 w 1 ≤ 4 . It’s clear that | z | = | 2 ± √ 1+2 w 1 3 - 2 w 1 | < 1 . (b) If 1 + 2 w 1 < 0 , i.e. w 1 < - 1 2 , we have z = 2 ± i √ - 1 - 2 w 1 3 - 2 w 1 , then | 2 ± i √ - 1 - 2 w 1 | = √ 3 - 2 w 1 , 3 - 2 w 1 > 4 . So | z | < 1 . 5