Survey of Numerical Stability Issues in Convergence Acceleration - PowerPoint PPT Presentation

Survey of Numerical Stability Issues in Convergence Acceleration AVRAM SIDI Computer Science Department Technion–Israel Institute of Technology Haifa, ISRAEL 1

Introduction Notation: { A m } ; A limit or antilimit of { A m } . { A ( j ) n } : Obtained by applying an extrapolation method to { A m } . When computed numerically (i.e., in floating-point arithmetic), the A i and hence the A ( j ) are in error. This may cause numerical instabilities, which eventually n may destroy the accuracy of the computed A ( j ) completely. n For almost all methods, K n K n A ( j ) γ ( j ) γ ( j ) � � = ni A j + i ; ni = 1 . n i =0 i =0 The γ ( j ) ni are functions of the ∆ A i . γ ( j ) = γ ( j ) ni + δ ( j ) If ¯ A i = A i + ǫ i , ¯ are the computed quantities, then the ni ni computed A ( j ) is n K n A ( j ) γ ( j ) ¯ ¯ � = ¯ A j + i . n ni i =0 2

Then K n K n A ( j ) − A | ≤ | A ( j ) | γ ( j ) | δ ( j ) | ¯ ni || ¯ � � − A | + ni || ǫ j + i | + A j + i | . n n i =0 i =0 If | ǫ i | = | A i || ρ i | and | δ ( j ) ni | = | γ ( j ) ni || η ( j ) ni | , where | ρ i | , | η ( j ) ni | ≤ u , then � K n K n � A ( j ) − A | ≤ | A ( j ) | γ ( j ) | γ ( j ) | ¯ � � ni || ¯ − A | + u ni || A j + i | + A j + i | . n n i =0 i =0 A i | = u | A i | + O ( u 2 ) , we have By A i ≈ ¯ A i and u | ¯ K n A ( j ) − A | � | A ( j ) | γ ( j ) | ¯ � − A | + 2 u ni || A j + i | . n n i =0 In case { A m } converges, A i ≈ A for all large i . Therefore, K n A ( j ) − A | � | A ( j ) | γ ( j ) | ¯ � − A | + 2 u | A | ni | . n n i =0 A ( j ) � | A ( j ) K n | ¯ − A | − A | n n | γ ( j ) � + 2 u ni | . | A | | A | i =0 3

Now, the theoretical relative error | A ( j ) − A | / | A | → 0 as j → ∞ or n → ∞ . n Consequently, for all large j or n , A ( j ) K n | ¯ − A | n | γ ( j ) � � 2 u ni | . | A | i =0 Thus, the quantity K n Γ ( j ) | γ ( j ) � = ni | ≥ 1 n i =0 plays a very important role when applying extrapolation methods in floating- point arithmetic. In case Γ ( j ) is of order 10 p and u is of order 10 − s , at most n A ( j ) s − p decimal digits of ¯ can be trusted. n Note: u is of order 10 − 16 for double precision, and of order 10 − 35 for quadru- A ( j ) ple precision. Therefore, better accuracy can be obtained from ¯ by in- n creasing the precision of the floating-point arithmetic, provided the A i are also computed in this arithmetic. We want sup j Γ ( j ) or sup n Γ ( j ) to be finite and small. In case Γ ( j ) is un- n n n bounded in j or n , we want to be able to force it to grow as slowly as possible. 4

Example. Levin u transformation Defined via the linear system n − 1 β i A l = A ( j ) � + ( l + 1)∆ A l ( l + 1) i , j ≤ l ≤ j + n. n i =0 Thus, i =0 ( − 1) i � n � � n ( j + i + 1) n − 2 ( A j + i /a j + i ) A ( j ) i = ; a i = A i − A i − 1 . n � n � � n i =0 ( − 1) i ( j + i + 1) n − 2 (1 /a j + i ) i ( − 1) i � n � ( j + i + 1) n − 2 (1 /a j + i ) γ ( j ) i ni = , i = 0 , 1 , . . . , n. � n � � n r =0 ( − 1) r ( j + r + 1) n − 2 (1 /a j + r ) r Best results are obtained from { A ( j ) n } ∞ n =0 with j fixed (such as j = 0 ). 5

Example. u transformation on A n = � n i =1 1 /i 2 , n = 1 , 2 , . . . , lim n →∞ A n = π 2 / 6 . E (0) E (0) Γ (0) ¯ ¯ n n (d) n (q) n 0 3 . 92 D − 01 3 . 92 D − 01 1 . 00 D + 00 2 1 . 21 D − 02 1 . 21 D − 02 9 . 00 D + 00 4 1 . 90 D − 05 1 . 90 D − 05 9 . 17 D + 01 6 6 . 80 D − 07 6 . 80 D − 07 1 . 01 D + 03 8 1 . 56 D − 08 1 . 56 D − 08 1 . 15 D + 04 10 1 . 85 D − 10 1 . 83 D − 10 1 . 35 D + 05 12 1 . 09 D − 11 6 . 38 D − 13 1 . 60 D + 06 14 2 . 11 D − 10 2 . 38 D − 14 1 . 92 D + 07 16 7 . 99 D − 09 6 . 18 D − 16 2 . 33 D + 08 18 6 . 10 D − 08 7 . 78 D − 18 2 . 85 D + 09 20 1 . 06 D − 07 3 . 05 D − 20 3 . 50 D + 10 22 1 . 24 D − 05 1 . 03 D − 21 4 . 31 D + 11 24 3 . 10 D − 04 1 . 62 D − 22 5 . 33 D + 12 26 3 . 54 D − 03 4 . 33 D − 21 6 . 62 D + 13 28 1 . 80 D − 02 5 . 44 D − 20 8 . 24 D + 14 30 1 . 15 D − 01 4 . 74 D − 19 1 . 03 D + 16 E (0) A (0) ¯ (d) : relative error in ¯ in double precision. ( ≈ 16 decimal digits) n n E (0) A (0) ¯ (q) : relative error in ¯ in quadruple precision. ( ≈ 35 decimal digits) n n 6

Improving stability Example. Levin–Sidi d (1) transformation Defined via the linear system, with 1 ≤ R 0 < R 1 < R 2 < · · · , n − 1 β i A R l = A ( j ) � + R l a R l , j ≤ l ≤ j + n ; a i = A i − A i − 1 . n R i i =0 l The A ( j ) and Γ ( j ) can be computed via the W algorithm: n n A Rj M ( j ) R j a Rj , N ( j ) R j a Rj , H ( j ) = ( − 1) j | N ( j ) 1 = = | , j ≥ 0 . 0 0 0 0 For n = 1 , 2 , . . . do M ( j +1) − M ( j ) N ( j +1) − N ( j ) H ( j +1) − H ( j ) M ( j ) , N ( j ) , H ( j ) n − 1 n − 1 n − 1 n − 1 n − 1 n − 1 = = = , j ≥ 0 . n n n R − 1 j + n − R − 1 R − 1 j + n − R − 1 R − 1 j + n − R − 1 j j j � � = M ( j ) H ( j ) A ( j ) , Γ ( j ) � � n n = � , j ≥ 0 . � � n n N ( j ) N ( j ) � � n n � end do Best results are obtained from { A ( j ) n } ∞ n =0 with j fixed (such as j = 0 ). 7

Example: Logarithmic convergence d (1) transformation on A n = � n i =1 1 /i 2 , n = 1 , 2 , . . . , lim n →∞ A n = π 2 / 6 . R 0 = 1 , R l = max { R l − 1 + 1 , ⌊ σR l − 1 ⌋} , with σ = 1 . 3 (GPS). E (0) E (0) Γ (0) ¯ ¯ n R n n (d) n (q) n 0 1 3 . 92 D − 01 3 . 92 D − 01 1 . 00 D + 00 2 3 1 . 21 D − 02 1 . 21 D − 02 9 . 00 D + 00 4 5 1 . 90 D − 05 1 . 90 D − 05 9 . 17 D + 01 6 7 6 . 80 D − 07 6 . 80 D − 07 1 . 01 D + 03 8 11 1 . 14 D − 08 1 . 14 D − 08 3 . 04 D + 03 10 18 6 . 58 D − 11 6 . 59 D − 11 3 . 75 D + 03 12 29 1 . 58 D − 13 1 . 20 D − 13 3 . 36 D + 03 14 48 1 . 55 D − 15 4 . 05 D − 17 3 . 24 D + 03 16 80 7 . 11 D − 15 2 . 35 D − 19 2 . 76 D + 03 18 135 5 . 46 D − 14 1 . 43 D − 22 2 . 32 D + 03 20 227 8 . 22 D − 14 2 . 80 D − 26 2 . 09 D + 03 22 383 1 . 91 D − 13 2 . 02 D − 30 1 . 97 D + 03 24 646 1 . 00 D − 13 4 . 43 D − 32 1 . 90 D + 03 26 1090 4 . 21 D − 14 7 . 24 D − 32 1 . 86 D + 03 28 1842 6 . 07 D − 14 3 . 27 D − 31 1 . 82 D + 03 30 3112 1 . 24 D − 13 2 . 52 D − 31 1 . 79 D + 03 E (0) A (0) ¯ (d) : relative error in ¯ in double precision. ( ≈ 16 decimal digits) n n E (0) A (0) ¯ (q) : relative error in ¯ in quadruple precision. ( ≈ 35 decimal digits) n n 8

Example: Linear convergence u and d (1) transformations on A n = � n i =1 z i − 1 /i , n = 1 , 2 , . . . , lim n →∞ A n = − log(1 − z ) /z ≡ f ( z ) ; z = 0 . 95 . R l = κ ( l + 1) , with κ > 0 integer (APS). ( κ = 1 gives the u transformation.) E (0) Γ (0) E (0) Γ (0) ¯ ¯ n n ( κ = 1) n ( κ = 1) n ( κ = 10) n ( κ = 10) 0 6 . 83 D − 01 1 . 00 D + 00 1 . 72 D − 01 1 . 00 D + 00 4 1 . 70 D − 01 3 . 92 D + 02 1 . 09 D − 04 1 . 49 D + 01 8 9 . 65 D − 03 1 . 64 D + 04 2 . 36 D − 08 9 . 92 D + 01 12 6 . 69 D − 04 7 . 86 D + 05 5 . 45 D − 12 6 . 71 D + 02 16 4 . 88 D − 05 3 . 87 D + 07 1 . 28 D − 15 4 . 55 D + 03 20 3 . 63 D − 06 1 . 93 D + 09 3 . 04 D − 19 3 . 09 D + 04 24 2 . 73 D − 07 9 . 66 D + 10 7 . 24 D − 23 2 . 10 D + 05 28 2 . 06 D − 08 4 . 85 D + 12 1 . 73 D − 26 1 . 43 D + 06 32 1 . 56 D − 09 2 . 44 D + 14 1 . 92 D − 28 9 . 73 D + 06 36 1 . 19 D − 10 1 . 23 D + 16 3 . 97 D − 27 6 . 62 D + 07 40 9 . 02 D − 12 6 . 17 D + 17 1 . 34 D − 26 4 . 51 D + 08 44 6 . 87 D − 13 3 . 11 D + 19 3 . 30 D − 26 3 . 07 D + 09 48 4 . 47 D − 14 1 . 57 D + 21 3 . 05 D − 25 2 . 09 D + 10 52 1 . 45 D − 12 7 . 89 D + 22 5 . 19 D − 25 1 . 42 D + 11 56 1 . 45 D − 11 3 . 98 D + 24 2 . 43 D − 23 9 . 67 D + 11 60 1 . 75 D − 09 2 . 01 D + 26 1 . 39 D − 22 6 . 58 D + 12 E (0) A (0) A ( j ) n : the computed A ( j ) ¯ n : relative error in ¯ ¯ n ( z ) . n ( z ) . E (0) (Computations in quadruple-precision.) Note: ¯ 24 ( κ = 20) = 2 . 21 D − 33 . 9

Example: Linear convergence (cont’d). u and d (1) transformations on A n = � n i =1 z i − 1 /i , n = 1 , 2 , . . . , lim n →∞ A n = − log(1 − z ) /z ≡ f ( z ) . R l = κ ( l + 1) , with κ > 0 integer (APS). ( κ = 1 gives the u transformation.) E (0) E (0) E (0) ¯ ¯ ¯ s z n ( κ = 1) “smallest” error 28 ( κ = 1) 28 ( κ = s ) 1 0 . 5 4 . 97 D − 30 ( n = 28) 4 . 97 D − 30 4 . 97 D − 30 0 . 5 1 / 2 ≈ 0 . 707 2 1 . 39 D − 25 ( n = 35) 5 . 11 D − 21 3 . 26 D − 31 0 . 5 1 / 3 ≈ 0 . 794 3 5 . 42 D − 23 ( n = 38) 4 . 70 D − 17 4 . 79 D − 31 0 . 5 1 / 4 ≈ 0 . 841 4 1 . 41 D − 21 ( n = 41) 1 . 02 D − 14 5 . 74 D − 30 0 . 5 1 / 5 ≈ 0 . 871 5 1 . 31 D − 19 ( n = 43) 3 . 91 D − 13 1 . 59 D − 30 0 . 5 1 / 6 ≈ 0 . 891 6 1 . 69 D − 18 ( n = 43) 5 . 72 D − 12 2 . 60 D − 30 0 . 5 1 / 7 ≈ 0 . 906 7 1 . 94 D − 17 ( n = 44) 4 . 58 D − 11 4 . 72 D − 30 E (0) A (0) A ( j ) n ( z ) are the computed A ( j ) ¯ = | ¯ ( z ) − f ( z ) | , and the ¯ n ( z ) . n n (Computations in quadruple-precision arithmetic.) Best results for z κ fixed and away from 1 , the point of singularity of f ( z ) . In our case, z κ = 0 . 5 maintained for every z > 0 in the table. Note that R l = ⌊ κ ( l + 1) ⌋ with non-integer κ > 1 also gives excellent results. 10