 Journal of Applied Mathematics and Physics, 2013, 1, 128-143 Published Online November 2013 (http://www.scirp.org/journal/jamp) http://dx.doi.org/10.4236/jamp.2013.15019 Open Access JAMP Higher Order Aitken Extrapolation with Application to Convergin g and Div e r g i ng Gauss-Seidel Iterations Ababu Teklemariam Tiruneh Department of Environmental Health Science, University of Swaziland, Mbabane, Swaziland Email: ababute@yahoo.com Received October 13, 2013; revised November 13, 2013; accepted November 20, 2013 Copyright © 2013 Ababu Teklemariam Tiruneh. This is an open access article distributed under the Creative Commons Attribu- tion License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. ABSTRACT In this paper, Aitken’s extrapolation normally applied to convergent fixed point iteration is extended to extrapolate the solution of a divergent iteration. In addition, higher order Aitken extrapolation is introduced that enables successive decomposition of high Eigen values of the iteration matrix to enable convergence. While extrapolation of a convergent fixed point iteration using a geometric series sum is a known form of Aitken acceleration, it is shown that in this paper, the same formula can be used to estimate the solution of sets of linear equations from diverging Gauss-Seidel iterations. In both convergent and divergent iterations, the ratios of differences among the consecutive values of iteration eventu-ally form a convergent (divergent) series with a factor equal to the largest Eigen value of the iteration matrix. Higher order Aitken extrapolation is shown to eliminate the influence of dominant Eigen values of the iteration matrix in suc-cessive order until the iteration is determined by the lowest possible Eigen values. For the convergent part of the Gauss-Seidel iteration, further acceleration is made possible by coupling of the extrapolation technique with the succes-sive over relaxation (SOR) method. Application examples from both convergent and divergent iterations have been provided. Coupling of the extrapolation with the SOR technique is also illustrated for a steady state two dimensional heat flow problem which was solved using MATLAB programming. Keywords: Linear Equations; Gauss-Seidel Iteration; Aitken Extrapolation; Acceleration Technique; Iteration Matrix; Fixed Point Iteration 1. Introduction Iterative solutions to systems of equations are widely employed for solving scientific problems. They have several advantages over direct methods. For equations involving large number of unknowns iterative solutions involve operations of order (n2) compared to direct solu-tions of order (n3). In computer applications, iterative solutions require much less memory and are quite simple to program . In addition iterative solutions are in many cases applicable to non-linear sets of equations. One set of iterative methods that are in wide use is centered on generation of Krylov sub space such as the method of conjugate gradients developed by Hestens and Stiefel . Such methods are guaranteed to converge in at most N steps for problems involving N unknowns. However, iterative solutions based on stationary methods such as Gauss-Seidel and others that are simpler to pro- gram and do not generate new iteration matrix are re- ceiving greater application . Iterative methods that mimic the physical processes involved such as the for-ward-backward method have been shown to produce solutions for some problems faster than Krylov methods . However, fixed point iterative methods are known to be less robust than Krylov methods and convergence is not guaranteed for ill-conditioned systems of equations. Iterative processes involving fixed point iterations that are convergent such as Jacobi and Gauss-Seidel iterations are known to form error series that are diminishing in the form of geometric series [5,14,16]. The geometric series sum based form of Aitken extrapolation for fixed point iteration involving the system of equations AXB is based on the formula : 1kkEXX where X is the Aitken extrapolation of the solution, Xk is the approximation to the solution at the kth iteration, Ek is the difference in X at the kth iteration 1–kkXX and  A. T. TIRUNEH 129is the ratio of the difference in X values 1kkEE. The geometric series extrapolation requires at least three points of the iteration. The extrapolation in terms of the three points at the k, k + 1 and k + 2 iterations takes the form . 21212kk kkkkXX XXX1kXXX which is a known form of Aitken extrapolation. Irons and Shrive  made a modification to Aitken’s method for scalars (sequences) which can also be applied to individual x values of fixed point iterations such as the Gauss-Seidel iteration. The extrapolation to the fifth point uses four points. After the four point extrapolation the fixed point iteration formula is used to obtain the fifth point. In this way by joint application of both extrapola-tion and fixed point iteration, the method is said to be dynamic with better convergence properties. A dynamic model in such a form does not require restarting the it-eration procedure as the method generates all necessary iterates from the latest extrapolation . The reduced rank extrapolation method  extends the scalar form of Aitken extrapolation into vectors of a given dimension. The extrapolation formulation for the iteration vector is a vector parallel of Aitken’s extrapola-tion: 1KKXXMIE where M is the iteration matrix of the fixed point itera-tion, I is the identity matrix of rank N, X are the iteration vector and Ek is the vector difference in X at the kth itera-tion. The method involves computing the generalized inverse involving the second order difference vector and is time consuming for solving non-linear systems of equations of larger size. Gianola and Schaffer  applied geometric series ex-trapolation for Jacobi and Gauss-Seidel iterations in ani-mal models. The optimal relaxation factor was lower when solutions were extrapolated, but its value was not as critical in the case of extrapolation. Fast Eigen vector computations that require matric in-version or decomposition are unsuitable for large size matrix problems as many of them involve operations of the order O(n3). Kamvar, et al.  applied Aitken ex-trapolation for accelerating page rank computations. They showed that Aitken acceleration computes the prin-cipal eigenvector of a Markov matrix under the assump-tion that the power-iteration estimate xk can be ex-pressed as a linear combination of the first two eigen-vectors. Calude Breziniski and Michela Redivo Zaglia  proposed extension of Aitken’s extrapolation into a gen-eral form involving transforming the sequence of itera-tion into a different form using known sequences which can lead to stabilization and convergence of the original iteration. The transformation, however, is not simple and straight forward and required further refinement. Chebyshev acceleration is also a way of transforming the iteration sequence which, for iteration matrix of known upper and lower bound Eigen values, the trans-formed sequence using Chebyshev polynomials leads to convergence of the fixed point iteration. In Chebyshev acceleration, the sequence of iteration values is modified by multiplication with Chebyshev polynomials which are constructed from known or estimated ranges of Eigen values of the iteration matrix. The choice of the form of Chebyshev polynomials is such that the procedure leads to progressive reduction of the norm of the error vector through a min-max property which minimizes the maxi-mum value that the polynomial has for the range of Ei-gen values specified . However, Chebyshev accelera-tion has the drawback of the need to accurately estimate the bounds of Eigen values of the iteration matrix, be-cause outside the domain of Eigen values the polynomial shows divergence and the min-max property does not hold. The following sections in this paper will give the basis of Aitken extrapolation for convergent fixed point itera-tions and the derivation for the extension of the extrapo-lation to divergent fixed point iterations. The paper con-tinues by introducing higher order Aitken extrapolation and shows how each order of extrapolation successively reduces the rank of Eigen values which helps in stabiliz-ing the extrapolation of a divergent iteration. Thereafter, application examples of both convergent and divergent fixed point iterations follow that include finite difference solution of Laplace equation to a two dimensional heat flow problem which was solved using MATLAB pro-gramming. 2. Method Development For solving a system of linear equations using fixed point iteration, the Aitken extrapolation formula can be written in the form: 1kkkexx where: x = The estimate of the solution at the limit of itera-tion; ek = The difference in consecutive x values, i.e., 1kkxx; k = The ratio of differences in x values, i.e., 121kkkkkkkexxexx1  It will now be shown that the above formula is appli-cable to both convergent as well as divergent Gauss- Seidel iteration. In addition it will also be shown that successive application of the Aitken extrapolation for- Open Access JAMP A. T. TIRUNEH Open Access JAMP 130 The differences in the solution vector Xk among con-secutive steps of iteration are formulated as follows: mula to a higher order will result in deflation of the dominant Eigen values one by one there by transforming a divergent iteration to a convergent form. Let the system of linearized equations for a given problem be represented in the matrix form: AXB where A is the coefficient matrix, B is the right hand side vector and X is the solution vector. Writing the matrix A further in terms of the components L, U and D matrices gives AULD where U and L are the upper and low triangular matrices respectively and D is the diagonal matrix. For the Gauss-Seidel iteration, the system of equations now can be written as: ULDX B LDXUX B Using the (k + 1)th and kth iteration X-vectors for the left and right hand sides of equation above respectively; 1kkLDXUX B Solving for 1kX;  111kkXLDUXLD B (1) This is the Gauss-Seidel iteration and the matrix; 1LD U  is called the iteration matrix. The actual iteration in terms of the x values (scalars) takes the form; 1111111inkkiiijjjjiii iiiibkijjxax axaa a    kNUBBk (2)  111kkXLDUXLD   1121kkXLDUXLD  121 1kk kXXLDUXX   In other words, 11kkELDUE  (4) For a convergent Gauss-Seidel iteration, the differ-ences in x values Ek written in terms of the difference between consecutive vectors of x values in Equation (4) above converges to the zero difference vectors when the solution X vector is reached. Denoting by M the iteration matrix; 1MLD U  and the difference vector of X among consecutive steps of the iteration by: 012 1,, ,, ,,kkEEE EE 1kkEME Assuming that the initial difference in x values (here after termed the error vector) E0 can be written as a linear combination of the Eigen vectors 12,,Nvvv of the itera-tion matrix M of size N with corresponding Eigen values 12,,,N in decreasing order of magnitude where 1 is the dominant Eigen value: 01 1122331011223311112223331112223 331111111222333oNNNNNNkkk kkNkkk kkNEXXcvcvcvcvEMEc MvcMvcMvcMvEcvcvcvc vEcvcvcvc vEcvcvcvc v        1NNNNN For the Gauss-Seidel iteration the vector of x values is successively computed from the fixed-point iteration process by:  111kkXLDUXLDBMX(3) where the matrix is the iteration matrix and . 1LDM1DBNL Taking the ratio of the Euclidean norms of Ek and Ek+1; 12222 2111 1111222333112222 21112223 331222112132111223311 1kkk kNN Nkkkk kkNN NkkkNNNcv c vcvcvEEcvc vcvcvcvcvcvc v       21k1222232111 223311 1kk kkNNNcvcvcvcv        2 A. T. TIRUNEH 131 If 1 is the dominant Eigen value, then the ratios 2121 1,,,N  tend to zero at higher k values so that 122111111122111lim limkkkkkkcvEEcv  (5) Therefore, for both converging and diverging Gauss- Seidel iterations, the error vector ratios are determined by the dominant Eigen value of the iteration matrix, 1. 2.1. Case I: Convergent Gauss-Seidel Iteration For a convergent iteration in which the dominant Eigen value is less than one, the difference vector Ek converges to zero. Denoting the estimate of the largest Eigen value of the iteration matrix M at the kth iteration by 1; and taking the difference among respective values of x as; 11kkee Starting with the kth difference in x values among con-secutive steps of iteration, a geometric series is formed in terms of the e values; 2311 1,, ,kkkkee e e  The sum of this diminishing geometric series is given by (since 1 < 1); 11ikiikee It is possible to write the e values so that; 11223kkkkkkkkkexxexxexx12 So that; 11ikikikeexx (6) Therefore, the solution x at convergence is extrapo-lated from Equation (6) above; 11kkexx (7) 2.2. Case II: Divergent Gauss-Seidel Iteration For a divergent Gauss-Seidel iteration, the error propaga-tion is studied as the iteration progresses. Let X0 be the initial estimate of the solution while Xs is the true solu-tion of the system of equation. The initial error vector Ero can then be written as: 0ro sEXX The difference between consecutive iteration values X is as defined before, i.e., 1kkEX Xk The Gauss-Seidel iteration formula given in Equation (3), i.e.,  111kkXLDUXLD B can also be rewritten as: 1kkXMX N where 1LDMU is the iteration matrix as defined before and 1.DBNLThe successive values of X of the iteration can now be written as follows; 10 0srsro sroXMXNM XENMXNME XME  The above expression is true because at the true solu-tion Xs, the Gauss-Seidel iteration satisfies the relation: ssXMX N Similarly for X2; 21 sroXMXNMXME N  222sro sroXMXNMEXM E Proceeding similarly at the kth iteration, the X-values can be written as: kksXXME ro (8) Interms of the difference between consecutive x-esti-mates, recalling the formula derived earlier in Equation (4), i.e., 1123 100 00001000kkkkkkikiELDUEMEXXEMEMEME MEXX ME   (9) It is seen from Equations (8) and (9) above that the X values increase in proportion to the iteration matrix M. Representing this increase in proportion to M by the largest Eigen value of the iteration matrix, 1, Equations (8) and (9) can now be written as: 110kkks rossXXEX XX (10) 10101000111kkikiEXX EX  (11) Open Access JAMP A. T. TIRUNEH 132 The expression in Equation (11) above is derived us-ing the general geometric series sum formula for an ex-panding geometric series. Equating the expressions for Xk in Equations (10) and (11): 0110 0111kkssEXXXX 01101111kksEXX 0001111sEEXX  0011sEXX In general for iteration estimation made from the kth it-eration vales of xk and individual ekvalues the formula can be written as: 11kskexx (12) It can be seen, therefore, that the same formula used for extrapolating a convergent Gauss-Seidel iteration can be used to extrapolate the solution from the divergent Gauss-Seidel iteration. Such a procedure which is un-conventional works well as the examples that follow il-lustrate. 2.3. Higher Order Aitken Extrapolation For the first-order Aitken extrapolation, the ratio of the norms of the error vector was shown in Equation (5) to be equal to the dominant Eigen value of the iteration ma-trix, i.e., 12211111122111kkkkcvEEcv For the second order Aitken extrapolation, the ex-trapolation made at the kth and (k + 1)th iterations are con-sidered:  11111kkskEXX   1111111kkskEXX The second order error in terms of the extrapolated X vectors can be written as:   2111111111121111Δ1ksk skkkkkkkkEX XEEXXEEE  Similarly for the (k + 1)th iteration;  121 1111Δ1kkkEEE At the kth iteration, the ratio of the norm of the error vector becomes;    1112111111Δ1Δ1kkkkkkEEEEEE  (13) In terms of the Eigen values  and Eigen vectors v of the iteration matrix M, the terms in the norm expression of Equation (13) above are given by: 11112 223 33kk kkNEcvcvcvcv  kNN1NNv 11111111222333kkk kkNEcvcvcvc     11111111222233 33Δ1111kk kkkkkNn NNEEEcvc vcvc   v   111121111111222211333 3Δ1111kkkkkkkNnN NEEEcvc vcvc   v Collecting the terms for both the numerator and de-nominator yields,   1111111111111Δ11111Δ111kNkikii iikNikii ikiEEcvEcvE(14) Examination of the terms on the right hand side of Equation (14) above reveals that the first terms of the summation (i.e. i = 1) in both the numerator and de-nominator vanish. Therefore, the second-order Aitken extrapolation reduces the Eigenvalue so that 2 becomes the dominant Eigen value. As will be shown below, the second dominant Eigen value of the iteration matrix, 2, will be equal to the ratio of the error vectors for the sec-ond order Aitken extrapolation. Factoring out the 2 term in Equation (14) above will Open Access JAMP A. T. TIRUNEH Open Access JAMP 133 result in the following expression:   11121122231121112222 3112111Δ11111Δ1111111kNkiikiiikkNkkiikiiiEcvc vEEEcvc v1     Once 1 has been eliminated and 2 is the dominant Eigen value, the ratios: 122,3,4,kkiii 5, being less than one, will vanish at higher k values so that the error norm ratios become    111122122211111222211212Δ11111Δ111kkkkkkkkkkEEcvEEEcvEEE    Similarly, it is easy to show that for the third and higher order Aitken extrapolations the error vector ratios correspond to the ith Eigen value of the iteration matrix. The higher order decomposition of Eigen values through higher order Aitken extrapolation can be generalized as: 1 for 1,2,3,,ikikEiENk (15) In effect, higher order Aitken extrapolation succes-sively decomposes the dominant Eigen values so that the error terms are determined eventually by the lowest Ei-gen value of the iteration matrix. This works for both the convergent and divergent Gauss-Seidel iteration. How-ever, the procedure is best in decomposing the first two dominant Eigen values beyond which the decomposition might be slow or inexact due to the successively small difference in the error vectors. The examples that follow later for diverging Gauss-Seidel iteration illustrate this fact. 2.4. Coupling of SOR Technique with Geometric Series Extrapolation Theextrapolation to the Gauss-Seidel iteration can well be extended to the successive over relaxation (SOR) method. In matrix form, the SOR iteration process is:  1111kXDL UDXDL B  (16) where  is the relaxation factor and the other terms are as deifned above. The iteration matrix is the coefficthe Xk term in Equation (16) above and is given byient of :  11MDL UD (17) The iteration formula for the successive over relaxa-tion technique in terms of individual x values is given by;   k111kkkxxbaxax(18) The acceleration factor  cannot be easily determinedin advance. It depends on the coefficient matrix A. If the coefficient matrix A is symmetric as well as posdefinite, the spectral radius of the iteration matrixw mentioned ea ext 1, 2, ,iiiijjijjji jiiiain itive M ill be less than one—ensuring convergence of the proc-ess—when the  value lies between 0 and 2. The procedure for extrapolation of the SOR process using geometric series sum, based on the dominant Eigen value of the iteration matrix as a ratio of the geometric series, follows a similar process to the onerlier. The only change is in the iteration matrix which is modified by the relaxation factor  while the condition for convergence (i.e. the dominant Eigen value of the iteration matrix being less than one) remains the same. However, it should be noted that the optimum relaxa-tion factor  is not necessarily the same as the SOR- optimum when the SOR technique is combined with Aitken extrapolation. This is illustrated in the applicationample of the heat flow problem presented in this paper. The opt for the SOR techniques is so chosen that the two dominant Eigen-values become equal in magnitude and these Eigen values at the optimum  can be complex numbers leading possibly to the failure of extrapolation methods. The fact that, at the optimum value of the ac-celeration factor , the dominant Eigen values turn out to be complex numbers is also shown in the heat flow ex-ample presented in this paper. Coupling the SOR tech-nique with the Aitken extrapolation at the exact optimum  is not necessary as the result is not very sensitive to the  value as will be shown in the heat flow example that follows. However, failure is not necessarily always the case for coupling at the optimum  value. The applica-ion example suggests that in the case of coupling of A. T. TIRUNEH 134 SThe first example below is a simple 2 × 2 equation in x eidel iteration starts at values distant from the solution, i.e. the starting values of x = 10,000 and y = 4250. The x and y values of the eration, differences in consecutive steps of the iteration e and the ratios, , are coue of the iteration ma- tri to −0.5). Therefore, the second x value, i.e., x OR with Aitken’s extrapolation, the  value can be chosen so that it is slightly less than the SOR optimum value enabling the coupling to be made at real Eigen values without the method failing to lead to convergence to the solution. 3. Examples from Convergent Iterations 3.1. Example 1 and y values. 27xy 2xyThe Gauss-Sitimputed and shown in Table 1. It is clear that the ratios of differences in x and y val-ues converge almost immediately to −0.50 for both the x and y values of the iteration. It will be later shown that, this value is the largest Eigen-valx, M. For calculating the extrapolated x value of the itera- tion, x = 10,000 cannot be used as the x0 value since the ratio at this level (−0.26) is not sufficiently convergent (compared 1= −2125.5 is taken. For this value the 0xe value is listed in Table 1 as 3186.75. For the y iteration y0 = 4250 can be taken as the ratio converged immediately with the first iteration. The value is also taken to be −6373.5. 0yThe x and y values are calculated as; e003186.752121.5 3.000110.5xxexx  006373.54250 1.000110.5yyeyy These values as expected are exactly equal to the solu-tion of the equations. To see that the ratios of alternative differences are the same as the largest Eigen value of the iteration matrix, the equation: 272xyxy is written as AXB; 21 711 2xy Using the relation A = L + D + U 21002001;; ;11 100100ALDU        ion of the Gauss-Seidel iteration for the 2 by 2 equation. Consecutive Difference ei (y) Ratio ei+1/ei (x) Ratio ei+1/ei (y) It can be shown that using the L, D and U matrices; Table 1. Computation for geometric series extrapolatIteration Step x y Consecutive Difference ei (x) 0 10000.00 4250.00−12121.50 −6373.50 −0.26 −0.50 1 50 1025 1065−2121.50 −2123.503186.75 3186.75 −0.50 −0.2 65.3.2−1593.38 −1593.38 −0.50 −0.50 3 −528.13 −530.13796.69 796.69 −0.50 −0.50 4 268.56 266.56−398.34 −398.34 −0.50 −0.50 5 −129.78 −131.78199.17 199.17 −0.50 −0.50 6 69.39 67.39 −99.59 −99.59 −0.50 −0.50 7 −30.20 −32.2049.79 49.79 −0.50 −0.50 8 19.60 17.60 −24.90 −24.90 −0.50 −0.50 9 −5.30 −7.30 12.45 12.45 −0.50 −0.50 10 7.15 5.15 −6.22 −6.22 −0.50 −0.50 11 0.93 −1.07 3.11 3.11 −0.50 −0.50 12 4.04 2.04 −1.56 −1.56 −0.50 −0.50 13 2.48 0.48 0.78 0.78 −0.50 −0.50 14 3.26 1.26 −0.39 −0.39 15 2.87 0.87 −2.87 −0.87 Open Access JAMP A. T. TIRUNEH 135 1012012LDU and that 1012012MLDU  The Eigen values of the iteration matrix M matrix are solving the determinant equation; found by 012det 0102MI  10  210 or 0.2 The largest Eigen value is −0.5 and it can be seen from Table 1 that the ratios of consecutive differences for both x and y variables converge to the largest Eigen value of the iteration matrix. 3.2. Exn vector is The iteration starts with the uss-Seidel iteration was carried out 13 times at which level four decimal-digit accuracy was obtained for thrences in x, y and z values. Ts and ratios of consecutive e values. can be shown that the iteration matrix M is given by; ample 2 The 2nd example is a 3 × 3 systems of linear equations given below 3523 442xyzxzxyz   The solutioT, , 1, 1, 1xyz vector: TT ,,10,8,5xyz  The Gae ratio of consecutive diffeable 2 shows the corresponding x, e valueIt is shown in Table 2 that, with the geometric series extrapolation a 5 digit accuracy has been obtained for the solution with just 9 iterations. The normal Gauss-Seidel iteration requires more than 60 iterations to arrive at 5 digit accuracy. For the coefficient matrix A of the given equation, it 11103311006682MLDU  The Eigen values are found by solving the determinant; 11011033110det 006611082MI  221038 0,0.1525793240.81924599or Therefore, the largest Eigen value of the iteration ma-trix M is  = −0.81924599 and all the ratios used above were approaching towards the largest Eigen value of the iteration matrix accurate to 5 digits as shon in Table 2. 3.3. Applicatsteel plate with dimension 10 × 20 cm has one of its 10 e wion Example (Heat Flow Problem) Example of application of the Aitken extrapolation me- thod for a steady state heat flow problem involving Laplace equation is given below . A rectangular thin cm edges held at 100˚C and the other three edges arheld at 0˚C. The thermal conductivity is given as k = 0.16 cal/sec·cm2·˚C/cm. Figure 1 shows the steel plate steady state conditions temperatures. The steady state heat flow problem is described by the Laplace equation: 22220uuxy with the boundary conditions, ,0,100,0 and 20,100Cuxuxu yuy . The nine-point finite difference formula of the Laplace equation is used for computation. is formula is sym-bolically represented by: ss- equation after 9 iterations. RatTh Table 2. Results geometric series extrapolation of the GauValues after 9 iterations Differences in values Seidel iteration for the 3 by 3ios of differences Extrapolated values x0 = 1.44846653 0xe= −0.81586443 x = −0.81923923 x = 1.000001908 y0 = 1.79206166 y = −0.81924816 y = 0.999998918 0ye = −1.44095868 z0 = 1.31013205 0ze = −0.56420578 z = −0.81924493 z = 1.000000209 Open Access JAMP A. T. TIRUNEH 136 2,6ijuh ,0iju 21414204141The ahe differen1lgebraic form of tce equation is: 1,1 1,1 1,11,121,1,, 1, 1,1644 44 20ij ij ijiji ji jijijijuuuuhuu u uu   0 wpera-tures at the intersection of the rectangular grid lines. Using a grid size of 2.5 cm, the 21 interior grid points shown in Figure 2 are generated. The coefficient matrix A of the linearized form of the Laplace equation AU = B by finite differencing is given near syitially. In addition the geh ions for the ndel itera-ton ed. In addition torms of ts for each step oere com- peometric series addition tctors and normctors, the consecutive difference ratios, as appn of the between 2.1 an here h is the grid size and the u values are the temin Table 3, followed by the right hand side vector B. A MATLAB program was written to solve the listem of equations for the 21 unknown temperatures using the Gauss-Seidel iteration inometric series extrapolation of the Gauss-Seidel was computed by writing the corresponding program in the MATLAB environment. The solution vector X for eacof the 100 iteratormal Gauss-Seiiwas computhe error vectorhe Euclidian nf the iteration wuted. For the gextrapolation, ino the solution ve of the error veroximatiomaximum Eigen-value were also computed. Figure 3 shows a comparison of the number of itera-tions required to arrive at more or less the same magni-tude of the norm of the error vector for the normal Gauss-Seidel iteration and for the extrapolation. It is ob-served from Figure 3 that the extrapolation procedure will give an acceleration factor in the range d 2.2. For example whereas 16 iterations are required for the normal Gauss-Seidel iteration giving error norm of 0.06+ , the geometric series extrapolation required only 7 iterations. For the subsequently smaller error norms, the numbers of iterations required are in the ratio of 16/7, 18/8, 20/9, 22/10, 24/11. It is clear that as the error norm gets smaller, the acceleration factor approaches a value of 2. The right hand side vector of the finite difference equation is given as: T,0, 600,0,0,0,0,0,0, 550  0,0,0,0,0,0, 550,0,0,0,0B,0 Figure 1. Rectangular plate for the steady state heat flow problem with boundary conditions. Figure 2. Rectangular plate for the heat flow problem with interior points and boundar y conditions shown. Open Access JAMP A. T. TIRUNEH 137Table 3. Coefficient matrix A of the finite difference form of the heat flow problem. −20 4 0 0 0 0 0 4 1 0 0 0 0 0 0 0 0 0 0 0 0 4 −20 4 0 0 0 0 1 4 1 0 0 0 0 0 0 0 0 0 0 0 0 4 −20 4 0 0 0 0 1 4 1 0 0 0 0 0 0 0 0 0 0 0 0 4 −20 4 0 0 0 0 1 4 1 0 0 0 0 0 0 0 0 0 0 0 0 4 −20 4 0 0 0 0 1 4 1 0 0 0 0 0 0 0 0 0 0 0 0 4 −20 4 0 0 0 0 1 4 1 0 0 0 0 0 0 0 0 0 0 0 0 4 −20 0 0 0 0 0 4 1 0 0 0 0 0 0 0 4 1 0 0 0 0 0 −20 4 0 0 0 0 0 4 1 0 0 0 0 0 1 4 1 0 0 0 0 4 −204 0 0 0 0 1 4 1 0 0 0 0 0 1 4 1 0 0 0 0 4 −204 0 0 0 0 1 4 1 0 0 0 0 0 1 4 1 0 0 0 0 4 −204 0 0 0 0 1 4 1 0 0 0 0 0 1 4 1 0 0 0 0 4 −204 0 0 0 0 1 4 1 0 0 1 0 0 0 0 1 4 −−− − − −−0 0 0 1 4 1 0 0 0 0 4 −204 0 0 0 0 1 4 0 0 1 4 0 0 0 0 0 4 −200 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 4 1 4 1 0 1 0 0 0 0 0 0 0 0 4 204 200 4 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 4 1 0 0 0 0 4 204 0 0 0 0 0 0 0 0 0 0 0 0 1 4 1 0 0 0 0 4 20 4 0 0 0 0 0 0 0 0 0 0 0 0 1 4 1 0 0 0 0 4 20 4 0 0 0 0 0 0 0 0 0 0 0 0 1 4 1 0 0 0 0 4 204 0 0 0 0 0 0 0 0 0 0 0 0 1 4 0 0 0 0 0 4 20 Figure 3. Comparison of number of iterations required for given norm of error vector between the normal Gauss-Seidel it-eration and the geometric series extrapolation. 3.4. Application Example—Double Acceleration Combined with the SOR Technique For the rectangular plate heat flow problem given above with 32 mesh divisions of interval cm in both x and y directions forming 21 interior points for the Gauss-Seidel iteration, the successive over relaxa- tion technique (SOR) was applied in conjunction with the geometric series sum based Aitken’s extrapolation. The optimum acceleration factor  to be used in Equa-tions (17) and (18) for rectangular regions with uniform boundary conditions (Drichlet conditions) is given by . 2.5h 2424opt c (19) The value of c in Equation (19) above is given by; πcos coscpq where p and q are the number of mesh divisions in the x and y directions. For the given problem p = 8 and q = 4, Open Access JAMP A. T. TIRUNEH 138 so that; ππcos cos1.6309884c 241.267241.63098opt Therefore, the minimum number of iterations required is achieved by using the optimum acceleration factor  of 1.267. This is shown in Table 4 for the SOR column where the minimum number of iteration of 35 was re-quired to reach to the solution vector within 10−15 accu-racy. The normal Gauss-Seidel process required 80 itera-tions whereas the geometric series extrapolation needed 47 iterations. Coupling of the SOR technique with geo-4, a reduc-by c extrapo-d on thetween SOR and the coupled iteration is insignifi-cant suggesting that a value  slightly less than the SOR optimum should be used if coupling interesting to note  that the largest Eigen values of the iteration matrix for the SOR technique, i.e., metric series extrapolation resulted in further reduction in the number of iterations required from 35 to 2tio by about 31% from the SOR result. The rate of reduction in number of iterations required oupling SOR with Aitken geometric seriesnlation is displayed in Figure 4 below for different values of the SOR acceleration factor,  basee values given in Table 4. It can be seen from the figure that sig-nificant reduction in the number of iterations required is achieved for  values between 1 and opt. In fact, the optimum value of  for the coupled iteration (SOR + Aitken geometric series extrapolation) lies below the SOR optimum for, i.e., at  = 1.23. Beyond the optimum acceleration factor, the differ-ence bis to be made. It is  11MDL UD  turn out to be a complex numbers at the optimum  value and beyond (refer to Table 4 for  = 1.267 and above). For all the  values above the SOR optimum, the corresponding largest Eigen values are complex numbers and the spectra radii are increasing. The progressive re-duction in largest Eigen values of the iteration matrix is also evident as the  value increases towards the SOR optimum. However, the extrapolation did not fail even if the dominant Eigen values were complex numbers at end beyond the optimum  values. The advantage of coupling the SOR technique with Gauss-Seidel iteration and as such does not involve extra 2 × 2 system of eq Aitken extrapolation is evident from example. In addition, the Aitken extrapolation is based on the the abovecalculation except generating a geometric series sum. The optimum SOR value is not always predictable. However, Aitken’s geometric series extrapolation can still work with or without the use of optimum  value. The above example shows that for acceleration factors slightly greater than one, significant reduction in the number of iterations required were obtained when the SOR was coupled with extrapolation. 4. Examples from a Divergent Gauss-Seidel Iteration 4.1. A 2 × 2 Equation with Diagonally Non-Dominant Coefficient Matrix The first example below is a simpleuations with diagonally non-dominant coefficient ma- Figure 4. Comparison of number of iterations required forgeometric series extrapolation coupled with SOR. In the figuRelaxation. given norm of error vector between the SOR method and the re, GE = Geometric series extrapolation. SOR = Successive Over Open Access JAMP A. T. TIRUNEH 139Table 4. Comparison of number of iterations required betweenwith the SOR technique. w SOR + GE SOR Extra acceleration Normal SOR and Aitken’s geometric series extrapolation coupled Largest Eigen value of the SOR iteration matrix 0.8 78 127 1.63 80 0.859905 0.9 59 102 1.73 80 1 47 80 1.70 80 1.05 40 70 1.75 80 1.1 34 62 1.82 80 1.15 26 54 2.08 80 1.2 25 45 0.827175 0.788581 1.23 24 40 1.67 80 1.25 26 35 1.35 80 1.267 27 35 1.30 80 0.760252 0.729593 0.691062 0.637938 1.80 80 Complex number 1.24 25 37 1.48 80 Complex number 1.3 32 37 1.16 80 Complex number 1.4 41 48 1.17 80 Complex number 1.6 72 76 1.06 80 Complex number 1.8 298 168 0.56 80 Complex number 0.59002 0.532422 trix The Gauss-Seidel iteration starts at values of x = 8 and y = 10. The x and y values of the iteration, differences in consecutive steps of the iteration ei and the ratios, , were computed and are shown in Table 5. As Table 5 shows, the Gauss-Seidel based iteration diverges as the atrix. However, the Aitken ex-trapolatios (shons otab innvo tx =y ng the es of econd iteration in Tab l e 5, the extra y val, xs and ysculated using Eqtion (12; 210 1215 510xyxy coefficient matrix is also diagonally non-dominant. This is also evident from the ratio of differences in x and y values shown in Tabl e 5. This ratio is −15 for both x and y iterations and it corresponds to the dominant Eigen value of the iteration mn iterationvariably cown inerge t the last columhe true solutions f the 1 and le)= 1.Usipolated valux and the sues are calua) as0010676 1.00011xsxexx 800150032 1.00011ysyeyy  se valus exped are exactly equal to the stion of the equations. 4.2. A 4 × 4 Equations with Diagonally Non-Dominant Coefficient Matrix A second example of a four by four system of equations in which the coefficient matrix is once again diagonally non-dominant is presented below. 4002026 15The es acteolu-12320234 1232078313656120 125346123 1207625127xxx420 125 14520481x  The normal Gauss-Seidel iteration quickly diverges as such. Higher order Aitken extrapolation has to be carried out sinmatrixfifth or-der Aitken extrapolatissfully achieves conver-gence whereas the firstinant Eigen values were successfully decomposed with and second order Aitken extrapolations rely. The solution of the s equations together with the Eigen values of the iteration matrix for the Gauss- Seidel iteration as obta program are given in Table 6. Figure 5 shows the sition of error ratios for the five-order Aitken extrtion. Comparing the ra- tio of successive errors t each level of extrapola-tion with the Eigen values of the iteration matrix given expected and it is not possible to reach to the solution as ce the two dominant Eigen values of the iteration are large (16.7 and 5.77 respectively). A on succe two domthe first spectiveystem ofined using MATLABdecompoapolagiven aOpen Access JAMP A. T. TIRUNEH 140 Table 5. Gauss-Seidel iteration results for a diagonally non-dominant and diverging 2 × 2 system of equations. Iteration Consecutive Difference) Consecutive Differenceei (y) Ratio ei + 1/ei(x) RExtrapolation (X) Extrapolation(Y) Step x y ei (xatio ei + 1/ei (y) 0 8 10 −52 −144 −13.846 4.49 1 −15 1 − − 2160 −15 1 1 62026 00 −32400 −15 1 1 −10124 −30374 162000 486000 −15 1 1 15145 000 −7290000 −15 1 1 −22 −64 0 1.09 × 108 −15 1 1 1876 102515626 −5.5E+08 −1.6 × 109 1 −5.1 × 108 −1.538 × 109 8.2 × 1009 2.46 × 10108 7.69 ×109 2.3066 × 10101.2 × 1011 −3.7 54 5.2 × 1012 44 134 720−15 2 76 −108 −15 3 −15 4 876 5626−2430 −15 5 7812483437 3645000−15 6 3417−15 −15 1 7 −15 −15 1 1 × 1011 −15 −15 1 1 × 1012 −9 −1.2 × 1011 −3.46 × 1011 1.85 × 1012 5.10 1.73 × 1012 5.1899 × 1012 −1.7 × 1012 −Table 6. Values of solution to the four by four system of equaprogram. Solutions of the systems of equation tions and Eigen values of the iteration matrix using MATLAB Eigen values of the iteration matrix X1 = 3.054225004761563 −16.700300071436452 X2 = −2.904223059942874 X3 = −0.661832433353327 X4 = −4.154545738306979 5.774221605390342 0.085523954767930 0 Figure 5. Variation of ratos of the error vectors at 1st to 5th oenapolationin ble 6ve shoat the first dominant Ei-gen valuesdecomd exactly −16.7003 and 2 = 5.7742). How the 3igher order extolasition slot eventually reduces stly g convef the itera-tion. F ur ss of rf the normof the error vector computed from irder Aitk extr. Ta abows thtwo are poseever, for(i.e. 1 = rd and hraption the decompows buufficien enablinrgence o ige 6 shows the progreeduction o EBAX As can be seenom Ta , with e fifth ordextrapolatthe no the erroector eveally down 0−9. Itld be recalled that tnor-mal Gauss-Seidel iteration is rapidly diverging for this m of equs. On tther hanigher order Ait-ken extrapolation as applied in this example successfully frble 7ther Ait-ken ion, rm ofr vntureduces to 1 shouhe systeationhe od hOpen Access JAMP A. T. TIRUNEH 141 F duction of the error vector for a 5th ord. Table 7. Progress of the 5th order Aitken extrapolation for the 4 ns. Iteration number X2 X3 rm of the error vector igure 6. Progress of reer Aitken extrapolation× 4 system of equatioX1 X4 No1 1 1 1 1 1689.086 2 3.044218606 −2.904617307 0.622362121 −4.153534798 8.14532 3 3.05422335 −2.904224136 −0.66182978 −4.15448631 0.00781 4 3.054225007 −2.904223059 −0.66183244 −4.154547438 0.000223 5 3.054225005 −2.90422306 −0.661832433 −4.15454574 6.37 × 10−6 6 3.054225005 −2.90422306 −0.661832433 −4.154545738 1.83 × 10−7 7 3.054225005 2.90422306 −0.661832433 −4.154545738 5.32 × 10−9 − converges to the solution of the system of equations. 4.3. A 6 × 6 System of Equations with a Diagonally Non -D ominant Coefficient Matrix A further example of diagonally non-dominant six by six systems of linear eThe dominant Eigen value of the Gauss-Seidel itera-tion matrix is −75.7966. Table 8 shows the exact solu-tions of the system of equations together with the Eigen values of the iteration matrix which were obtained from a MATLAB program. With the combination of dominant Eigen values shown in Table 8, the normal Gauss-Seidel iteration quickly diverges. However, a 4th order Aitken extrapolation was enough to bring the iteration to con-vergence. Figure 7 shows the decomposition of the ratio of error vectors at each order of Aitken extrapolation. At the first and second order extrapolation the ratio of the errors are exactly equal to the first two dominant Eigen values of the iteration matrix (i.e., 75.7966 and 11.7041). However, the third and higher order extolations decompose ling convergence The variation of the norm of the error vector with the number of fifth order Aitken iteratio is given in Figure 8.series making it suitable for extrapolation through ex-amination of the ratios of consecutive differences in the solution vector at each step of the iteration. The ratio belongs to the largest Eigen value of the iteration matrix. When sufficient digits of accuracy are obtained for the ratio, the process can be extrapolated towards the solu-tion using a geometric series sum. This procedure is the equivalent of Aitken extrapolation for a convergent itera-tion. A significant reduction in the number of iteration required is obtained through such extrapolation. Cou- quations is given below: slowly to successively lower ratio enabof the Aitken extrapolation. 6T02000 80008746577830 3460 1270 4810x1234540035432 104820035600485 30202000196 10545483497403048 6312034748202034 5457680xxxx 7623.8 2690x rapsn As can be seen from the figure, the norm of the error vector reduces quickly with the first few iterations. 5. Conclusions The Gauss-Seidel iteration in the case of a convergent iteration is known to follow a diminishing geometric Open Access JAMP A. T. TIRUNEH 142 Table 8. Values of solution to the six by six system of equations and Eigen values of the iteration matrix. Solutions of the systems of equation Eigen values of the iteration matrix X1 = −0.563147393304287 −75.796638515975403 X2 = −0.731832157439239 −11.704130560629785 X3 = 1.857839885254053 −0.368124943597328 X4 = −6.666186315796603 −0.027943199473836 X5 =−1.725114498846410 −0.000000000000001 X6 = −1.833877505834098 0 Figuion of raror vecth orderolation. re 7. Variattios of the ertors at 1st to 5 Aitken extrap ctor for a 4th order Aitken extrapolation. less than the optimum as the heat flow example presFigure 8. Progress of reduction of the error ve pling of the successive over relaxation technique with Aitken’s extrapolation is possible with further reduction in iteration while employing relaxation factors not nec-escoupling with SOR technique is done at  value typically ented . ergent Gauss-Seidel iterations the ap-sarily restricted the optimum value which may be dif-ficult to predict in advance for some types of equations. Coupling of extrapolation with SOR technique is nor-mally not always possible at the optimum acceleration factor w because the largest Eigen value at this optimum value can turn out to be a complex number. Therefore, in this paper showedIn the case of divplication of Aitken extrapolation formula is made possi-ble and in many cases the extrapolation at each level of the Gauss-Seidel iteration indicates convergence towards the solution. Higher order Aitken extrapolation succes-sively decomposes the dominant Eigen values of the it-eration matrix. In doing so, the iteration is successively transformed from an expanding (divergent) form to a Open Access JAMP A. T. TIRUNEH 143stable convergent iteration. At each stage of the applica-tion of higher order Aitken extrapolation, the ratio of the error vectors (differences in successive x values) ap-proaches the dominant Eigen value for that order of ex-trapolationan interesting possibility of stabilizing, i.e., converting a divergent fixd convergent iteration. In genelate the solution from divergent as an interesting possibility that expe of application of fixed point iteratiohampered by problems of di H. L. Stone, “Iterative Solution of Implicit Approxima- tions of Multidimensional Partial Differential Equations,” SIAM Journal on Numerical Analysis, Vol. 5, No. 3, 1968, pp. 530-558. http://dx.doi.org/10.1137/0705044. Higher order Aitken extrapolation provides e point iteration to a stable,ral, the ability to extrapoGuss-Seidel iteration iands further the scopn which in some cases is vergence. REFERENCES  M. R. Hestenes and E. Stiefel, “Methods of Conjugate Gradients for Solving Linear Systems,” Journal of Re- search at the National Bureau of Standards, Vol. 49, 1952, pp. 409-436. http://dx.doi.org/10.6028/jres.049.044  G. H. Golub and C. F. Van Loan, “Matrix Computations,” 3 Edition, Johns Hopkins University Press, Baltimore, 1996.  J. C. West, “Preconditioned Iterative Solution of Scatter- ing from Rough Surfaces,” IEEE Transactions on Anten- nas and Propagation, Vol. 48, No. 6, 2000, pp. 1001- 1002. http://dx.doi.org/10.1109/8.865236  A. C. Aitken, “Studies in Practical Mathematics. The Eva- luation of LatProceedings of the Royal Society of Edinburgh, Vol. 57, 1937, p. 269. 70, No. 12, 1987, pp. 2577-2584. s in Engi- un,” John Wiley S. R. Capehelerating Iterative Me- thods for thical Problems,” Ph.D. DissertationOrsity, 1989.  P. Henrici, l Analysis,” John Wi- ley & Sons, S. D. Kamv,. D. Manning and G. H. Golub, “Extrapolatiothods for Accelerating Page Rank Computations,” gs of the Twelfth Interna- tional World Wide Web Conference, 2003, pp. 261-270. http://dx.doi.org/10.1145/775152.775190ent Roots and Latent Vectors of a Matrix,”  M. D. Gianola and L. R. Schaffer, “Extrapolation and Convergence Criteria with Jacobi and Gauss-Seidel Itera- tion in Animal Models,” ournal of Diary Science, Vol. J B. Irons and N. G. Shrive, “Numerical Methodneering and Applied Science: Numbers are F & Sons, New York, 1987. art, “Techniques for Acce Solution of Mathemat, klahoma State Unive“Elements of Numerica New York, 1964. ar T. H. Haveliwala, Cn MeProceedin  C. Breziniski and M. R. Zaglia, “Generalizations of Ait- ken’s Process for Accelerating the Convergence of Se- quences,” Journal of Computational and Applied Mathe- matics, Vol. 26, No. 2, 2007, pp. 171-189.  A. Bjorck, “Numerical Methods in Scientific Comput- ing,” Linkoping University Royal Institute of Technology, Vol. 2, 2009.  D. M. Young, “Iterative Solution of Large Linear Sys- tems,” Academic Press, New York, 1971.  G. Dahlquist and A. Bjorck, “Numerical Methods,” Pren- tice Hall, Englewood Cliffs, 1974.  C. F. Gerald and P. O. Wheatley, “Applied Numerical Analysis,” 5th Edition, 1994. plied Iterative Me- thods,” Academic Press, New York, 1981.  L. A. Hageman and D. M. Young, “ApOpen Access JAMP