Von Neumann stability analysis In numerical algorithms for - PDF document

Von Neumann stability analysis In numerical algorithms for differential equations the concern is the growth of round-off errors and/or initially small fluctuations in initial data which might cause a large deviation of final answers from the exact solution. The method of stability analysis shown next was developed by the mid-twentieth century Hungarian mathematician — and father of the electronic computer — John von Neumann. The von Neuman stability analysis is based on the decomposition of numerical errors of numerical approximations into Fourier series. Fourier series decomposes any periodic function or periodic signal into the sum of a (possibly infinite) set of simple oscillating functions, namely sines and cosines (or, equivalently, complex exponentials). We have chosen complex exponentials to represent errors as they are much easier to work with than real trigonometric functions. Consider the one-dimensional heat equation 𝜖𝑢 = 𝛽 𝜖 2 𝑉 𝜖𝑉 𝜖𝑦 2 [1] Stability Analysis of FTCS scheme. Approximate the numerical error. Using the FTCS method the discretized form of equation [1] is: 𝑜+1 = 𝑉 𝑜 + 𝜇(𝑉 𝑜 + 𝑉 𝑜 ) [2] 𝑜 𝑉 − 2𝑉 𝑘 𝑘 𝑘+1 𝑘 𝑘−1 where 𝜇 = 𝛽∆𝑢/∆𝑦 2 (called the mesh ratio), 𝑘 is the space index, and 𝑜 is the time index. Define the error in the numerical approximation as 𝑜 = 𝑉 𝑜 − 𝑣 𝑘 𝑜 [3] 𝜗 𝑘 𝑘 𝑜 is the exact solution at grid point (𝑘, 𝑜) . Both the numerical solution 𝑉 𝑜 and the exact solution where 𝑣 𝑘 𝑘 𝑜 satisfy equation [2], therefore, the error 𝜗 𝑘 𝑜 also follows the discretized ODE equation [2] 𝑣 𝑘 𝑜+1 = 𝜗 𝑘 𝑜 + 𝜇(𝜗 𝑘+1 𝑜 + 𝜗 𝑘−1 𝑜 ) [4] 𝑜 𝜗 𝑘 − 2𝜗 𝑘 The equations [2] and [4] show that both the error and the numerical solution have the same growth or decay behavior with respect to time. For linear differential equations with periodic boundary condition, the spatial variation of error may be expanded in a finite Fourier series, in the interval L, as 𝑁 𝜗(𝑦) = ∑ 𝐵 𝑛 𝑓 𝑗𝑙 𝑛 𝑦 [5] 𝑛=1 𝜌𝑛 where the wavenumber 𝑙 𝑛 = 𝑀 where 𝑛 = 1,2, … , 𝑁 and 𝑁 = 𝑀/∆𝑦 , 𝑗 = √−1 , and 𝑦 is the independent space variable. 𝑓 𝑗𝑙 𝑛 𝑦 is the complex exponential (recall the Euler’s formula: 𝑓 𝑗𝜒 = cos 𝜒 + 𝑗 sin 𝜒 , where 𝑓 is the base of natural logarithm and 𝜒 is a real number). [Von Neumann stability analysis-Presentation-V5] Page 1 of 6

The time dependence of the error is included by assuming the amplitude of error 𝐵 𝑛 is a function of time. Since the error tends to grow or decay exponentially with time, it is reasonable to assume that the amplitude varies exponentially with time; hence 𝑁 𝑜 = ∑ 𝑓 𝑏𝑢 𝑓 𝑗𝑙 𝑛 𝑦 𝜗(𝑦, 𝑢) = 𝜗 𝑘 [6] 𝑛=1 where 𝑏 is a constant. Since the difference equation for error is linear (the behavior of each term of the series is the same as the series itself), it is enough to consider the growth of error of a typical term: 𝑜 = 𝑓 𝑏𝑢 𝑓 𝑗𝑙 𝑛 𝑦 [7] 𝜗 𝑛 (𝑦, 𝑢) = 𝜗 𝑘 The stability characteristics can be studied using just this form for the error with no loss in generality. Define a Stability Criterion . The goal is to show that the error incurred by a particular numerical scheme doesn’t grow in the evolving steps of time. Define the amplification factor 𝑜+1 𝜗 𝑘 𝐻 ≡ 𝑜 [8] 𝜗 𝑘 Then, the necessary and sufficient condition for the error to remain bounded, maintaining numerical stability is |𝐻| ≤ 1 𝑝𝑠 [−1 ≤ 𝐻 ≤ 1 ] [9] How the FTCS scheme complies with the Stability Criterion . The goal now is to show that equation [9] holds for the particular FTCS numerical scheme. To find out how the error varies in steps of time, substitute equation [7] into each term of equation [4], as shown below 𝑜 = 𝑓 𝑏𝑢 𝑓 𝑗𝑙 𝑛 𝑦 [10.𝑏] 𝜗 𝑘 𝑜+1 = 𝑓 𝑏(𝑢+∆𝑢) 𝑓 𝑗𝑙 𝑛 𝑦 [10.𝑐] 𝜗 𝑘 𝑜 = 𝑓 𝑏𝑢 𝑓 𝑗𝑙 𝑛 (𝑦+∆𝑦) [10. 𝑑] 𝜗 𝑘+1 𝑜 = 𝑓 𝑏𝑢 𝑓 𝑗𝑙 𝑛 (𝑦−∆𝑦) [10.𝑒] 𝜗 𝑘−1 By eqs. [10] in equation [4], the last yields 𝑓 𝑏(𝑢+∆𝑢) 𝑓 𝑗𝑙 𝑛 𝑦 = 𝑓 𝑏𝑢 𝑓 𝑗𝑙 𝑛 𝑦 + 𝜇(𝑓 𝑏𝑢 𝑓 𝑗𝑙 𝑛 (𝑦+∆𝑦) − 2𝑓 𝑏𝑢 𝑓 𝑗𝑙 𝑛 𝑦 + 𝑓 𝑏𝑢 𝑓 𝑗𝑙 𝑛 (𝑦−∆𝑦) ) [11] Divide by 𝑓 𝑏𝑢 𝑓 𝑗𝑙 𝑛 𝑦 to yield after simplification 𝑓 𝑏∆𝑢 = 1 + 𝜇(𝑓 𝑗𝑙 𝑛 ∆𝑦 + 𝑓 −𝑗𝑙 𝑛 ∆𝑦 − 2) [12] Using the identities (see Appendix) [Von Neumann stability analysis-Presentation-V5] Page 2 of 6

𝑗𝑙𝑛∆𝑦 −𝑓 − 𝑗𝑙𝑛∆𝑦 𝑙 𝑛 ∆𝑦 𝑓 2 2 𝑡𝑗𝑜 ( 2 ) = [13] , and 2𝑗 ) = − [𝑓 𝑗𝑙 𝑛 ∆𝑦 + 𝑓 −𝑗𝑙 𝑛 ∆𝑦 − 2] sin 2 (𝑙 𝑛 ∆𝑦 [14] 2 4 By eq. [14] in eq. [12], the last may be written as 𝑓 𝑏∆𝑢 = 1 − 4λ sin 2 (𝑙 𝑛 ∆𝑦 ) [15] 2 However, by eqs [10] in eq [8], the last yields 𝑜+1 𝑜 = 𝑓 𝑏(𝑢+∆𝑢) 𝑓 𝑗𝑙 𝑛 𝑦 𝜗 𝑘 = 𝑓 𝑏∆𝑢 [16] 𝐻 = 𝑓 𝑏𝑢 𝑓 𝑗𝑙 𝑛 𝑦 𝜗 𝑘 Thus, by equations [15], [16] plugged into [9], the condition for stability is given by (1 − 4𝜇 sin 2 (𝑙 𝑛 ∆𝑦 )) ≥ −1 [17.𝑏] 2 [17] |𝐻| = |1 − 4𝜇 sin 2 (𝑙 𝑛 ∆𝑦/2)| ≤ 1 { (1 − 4𝜇 sin 2 (𝑙 𝑛 ∆𝑦 )) ≤ 1 [17.𝑐] 2 For the above conditions, equations [17.a] and [17.b], to hold for all m (and therefore all sin 2 (𝑙 𝑛∆𝑦 /2) ), we need to consider [17.a] and [17.b], separately. Equation [17.a] yields [18] 4𝜇 sin 2 (𝑙 𝑛 ∆𝑦/2) ≤ 2 Note that the term 4𝜇 sin 2 (𝑙 𝑛 ∆𝑦/2) is always positive, i.e., always in the range [0,1]. The worst case is when sin 2 (𝑙 𝑛 ∆𝑦/2) = 1 , therefore, equation [18] yields 1 𝛽∆𝑢 1 [19] 𝜇 ≤ 2 or 𝜇 = ∆𝑦 2 ≤ 2 Equation [17.b] yields [20] 4𝜇 sin 2 (𝑙 𝑛 ∆𝑦/2) ≥ 0 Since sin 2 (𝑙 𝑛 ∆𝑦/2) range is [0,1] equation [20] yields [21] 𝜇 ≥ 0 1 which always holds. Hence, combining [19] and [21]: 0 ≤ 𝜇 ≤ 2 Equation [19] gives the stability requirement for the FTCS scheme as applied to the one-dimensional heat equation. We say that the method is conditionally stable ; for a given ∆𝑦 , the allowed value of ∆𝑢 must be small enough to satisfy equation [19] or [Von Neumann stability analysis-Presentation-V5] Page 3 of 6

[22] ∆𝑢 ≤ ∆𝑦 2 2𝛽 Usually we establish the space discretization by assuming ∆𝑦 , then we compute the allowed ∆𝑢 with equation [22 ] so to be sure our method won’t ex ceed the stability requirement. For instance, if we have ∆𝑦 = 0.01 , and 𝛽 = 1 , then we can only use a time step size ∆𝑢 ≤ 0.00005 , which is minuscule. It would take an inordinately large number of time steps to compute the value of the solution at even moderate final times, e.g., t = 1. Moreover, owing to the limited accuracy of computers, the propagation of round-off errors might then cause a significant reduction in the overall accuracy of the final solution values. Since not all choices of space and time steps lead to a convergent scheme, that’s the reason the explicit scheme FTCS is called conditionally stable , capisce my friend? Stability Criterion for the Crank-Nicolson Scheme By the CN scheme, the discretized form of equation [1] 𝑙 = 𝜇 𝑙+1 ] + 𝜇 𝑙+1 − 𝑉 𝑗 𝑙+1 − 2𝑉 𝑗 𝑙+1 + 𝑉 𝑗+1 𝑙 + 𝑉 𝑗+1 𝑙 ] [23] 𝑙 𝑉 𝑗 2 [𝑉 𝑗−1 2 [𝑉 𝑗−1 − 2𝑉 𝑗 𝑛 also follows the discretized ODE equation Using similar arguments as in the FTCS method, the error 𝜗 𝑘 of the CN scheme 𝑙 = 𝜇 𝑙+1 ] + 𝜇 𝑙+1 − 𝜗 𝑗 𝑙+1 − 2𝜗 𝑗 𝑙+1 + 𝜗 𝑗+1 𝑙 + 𝜗 𝑗+1 𝑙 ] [24] 𝑙 𝜗 𝑗 2 [𝜗 𝑗−1 2 [𝜗 𝑗−1 − 2𝜗 𝑗 𝑙 Divide [24] by 𝜗 𝑗 𝑙+1 𝑙 𝑙+1 𝑙+1 𝑙+1 𝑙 𝑙 𝑙 𝜗 𝑗 𝑙 − 𝜗 𝑗 𝑙 = 𝜇 2 [𝜗 𝑗−1 𝑙 − 2𝜗 𝑗 + 𝜗 𝑗+1 𝑙 ] + 𝜇 2 [𝜗 𝑗−1 𝑙 − 2𝜗 𝑗 𝑙 + 𝜗 𝑗+1 𝑙 ] [25] 𝑙 𝜗 𝑗 𝜗 𝑗 𝜗 𝑗 𝜗 𝑗 𝜗 𝑗 𝜗 𝑗 𝜗 𝑗 𝜗 𝑗 Simplify each term in equation [25]: 𝑙+1 𝑙 = 𝐻 = 𝑓 𝑏(𝑢+∆𝑢) 𝜗 𝑗 [26. 𝑏] 𝑓 𝑏𝑢 𝜗 𝑗 𝑙+1 = 𝑓 𝑏(𝑢+∆𝑢) 𝑓 𝑗𝑙 𝑛 (𝑦−∆𝑦) = 𝐻 𝑓 𝑗𝑙 𝑛 𝑦 𝑓 𝑗𝑙 𝑛 (−∆𝑦) 𝜗 𝑗−1 = 𝐻𝑓 −𝑗𝑙 𝑛 ∆𝑦 [26. 𝑐] 𝑙 𝑓 𝑏𝑢 𝑓 𝑗𝑙 𝑛 𝑦 𝑓 𝑗𝑙 𝑛 𝑦 𝜗 𝑗 𝑙+1 𝑙 = 𝑓 𝑏(𝑢+∆𝑢) 𝑓 𝑗𝑙 𝑛 (𝑦+∆𝑦) = 𝐻 𝑓 𝑗𝑙 𝑛 𝑦 𝑓 𝑗𝑙 𝑛 ∆𝑦 𝜗 𝑗+1 = 𝐻𝑓 𝑗𝑙 𝑛 ∆𝑦 [26. 𝑑] 𝑓 𝑏𝑢 𝑓 𝑗𝑙 𝑛 𝑦 𝑓 𝑗𝑙 𝑛 𝑦 𝜗 𝑗 𝑙 𝑙 = 𝑓 𝑏𝑢 𝑓 𝑗𝑙 𝑛 (𝑦−∆𝑦) = 𝑓 𝑗𝑙 𝑛 𝑦 𝑓 𝑗𝑙 𝑛 (−∆𝑦) 𝜗 𝑗−1 = 𝑓 −𝑗𝑙 𝑛 ∆𝑦 [26. 𝑒] 𝑓 𝑏𝑢 𝑓 𝑗𝑙 𝑛 𝑦 𝑓 𝑗𝑙 𝑛 𝑦 𝜗 𝑗 [Von Neumann stability analysis-Presentation-V5] Page 4 of 6