Convergent Homotopy Analysis Method for Solving Linear Systems

Abstract

By using homotopy analysis method (HAM), we introduce an iterative method for solving linear systems. This method (HAM) can be used to accelerate the convergence of the basic iterative methods. We also show that by applying HAM to a divergent iterative scheme, it is possible to construct a convergent homotopy-series solution when the iteration matrix G of the iterative scheme has particular properties such as being symmetric, having real eigenvalues. Numerical experiments are given to show the efficiency of the new method.

1. Introduction

Computational simulation of scientific and engineering problems often depend on solving linear system of equations. Such systems frequently arise from discrete approximation to partial differential equations. Systems of linear equations can be solved either by direct or by iterative methods. Iterative methods are ideally suited for solving large and sparse systems. For the numerical solution of a large nonsingular linear system, 𝐴𝑢=𝑏, (1) where 𝐴∈ℝ^𝑛×𝑛 is given, 𝑏∈ℝ^𝑛 is known, and 𝑢∈ℝ^𝑛 is unknown, one class of iterative methods is based on a splitting (𝑀,𝑁) of the matrix 𝐴, that is, 𝐴=𝑀−𝑁, (2) where 𝑀 is taken to be invertible and cheap to invert, which mean that a linear system with matrix coefficient 𝑀 is much more economical to solve than (1). Based on (2), (1) can be written in the fixed-point form 𝑢=𝐺𝑢+𝑐,𝐺=𝑀⁻¹𝑁,𝑐=𝑀⁻¹𝑏, (3) which yields the following iterative scheme for the solution of (1): 𝑢^(𝑘+1)=𝐺𝑢^(𝑘)𝑢+𝑐,𝑘=0,1,2,…,₍₀₎∈ℝ^𝑛isarbitrary. (4)

A sufficient and necessary condition for (4) to converge to the solution of (1) is 𝜌(𝐺)<1, where 𝜌(𝐺) denotes spectral radius. Some effective splitting iterative methods and preconditioning methods were presented for solving the linear system of (1), see [1–9]. Recently, Keramati [10], Yusufoğlu [11], and Liu [12] applied the homotopy perturbation method to obtain the solution of linear systems and deduced the conditions for checking the convergence of homotopy series. In this work, we show how the homotopy analysis method may be regarded as an acceleration procedure based on the iterative method (4). We observe that it is not necessary that the basic method (4) be convergent. When 𝜌(𝐺)>1, it is sufficient that the eigenvalues 𝜆_𝑖=Re(𝜆_𝑖)+𝑖Im(𝜆_𝑖), 𝑖=1,…,𝑛 of iteration matrix 𝐺 satisfy the relation Re(𝜆_𝑖)<1, 𝑖=1,…,𝑛, (or Re(𝜆_𝑖)>1, 𝑖=1,…,𝑛). When 𝜌(𝐺)<1, by applying the homotopy analysis method to the basic iterative method (4), we can improve the rate of convergence of the iterative method (4). This paper is organized as follows. In Section 2, we introduce the basic concept of HAM, derive the conditions for convergence of the homotopy series, and apply the homotopy analysis method to the Jacobi, Richardson, SSOR, and SAOR methods. In Section 3, some numerical examples are presented to show the efficiency of the method. Finally, we make some concluding remarks in Section 4.

2. Basic Idea of HAM

The homotopy analysis method (HAM) [13, 14] was first proposed by S. J. Liao in 1992. The HAM was further developed and improved by S. Liao for nonlinear problem in [15].

Here, we apply the homotopy analysis method (HAM) to the problem (3) for finding the solution of (1) when det(𝐴)≠0. Consider (3), where 𝑢 is unknown vector of (1) and 𝐺 is the iteration matrix of an iterative method. Let 𝑤₀ denote an initial guess of exact solution 𝑢, ℏ≠0 an convergence control parameter. Then, we can apply the homotopy analysis method and define 𝐻(𝑢,𝑝) by 𝐻(𝑢,0)=𝐹(𝑢),𝐻(𝑢,1)=𝐿(𝑢) (5) and a homotopy as follows: 𝐻𝑢,𝑝=1−𝑝𝐹(𝑢)−𝑝ℏ𝐿(𝑢)=0, (6) where 𝐿(𝑢)=(𝐼−𝐺)𝑢−𝑐,𝐹(𝑢)=𝑢−𝑤₀, (7) and 𝑝 is an embedding parameter. Hence, it is obvious that 𝐻(𝑢,0)=𝑢−𝑤₀,𝐻(𝑢,1)=(𝐼−𝐺)𝑢−𝑐. (8) And as the embedding parameter 𝑝∈[0,1] increases from 0 to 1, the solution of 𝐻(𝑢,𝑝) varies continuously from the initial approximation 𝑤₀ to the exact solution 𝑢 of the original equation 𝐿(𝑢)=0. The homotopy analysis method uses the parameter 𝑝 as an expanding parameter (see [16–18]) to obtain 𝑢=𝑢₀+𝑝𝑢₁+𝑝²𝑢₂+⋯, (9) and it gives an approximation to the solution of (3) as 𝑣=lim_𝑝→1𝑢₀+𝑝𝑢₁+𝑝²𝑢₂+⋯=lim^∞_{𝑝→1𝑘=0}𝑝^𝑘𝑢_𝑘=^∞_𝑘=0𝑢_𝑘. (10)

By substituting (9) in (6) and equating the terms with identical powers of 𝑝, we can obtain 𝑝⁰:𝑢₀−𝑤₀𝑝=0,₁:𝑢₁−ℏ(𝐼−𝐺)𝑢₀𝑝−𝑐=0,_𝑖:𝑢_𝑖−𝑢_𝑖−1−ℏ(𝐼−𝐺)𝑢_𝑖−1=0,𝑖=2,3,…. (11)

This implies that 𝑢₀=𝑤₀,𝑢₁=ℏ(𝐼−𝐺)𝑢₀,𝑢−𝑐_𝑖=[(ℏ+1)𝐼−ℏ𝐺]𝑢_𝑖−1,𝑖=2,3,…. (12)

Taking 𝐺_ℏ=(ℏ+1)𝐼−ℏ𝐺 (13) yields that 𝑢₁=ℏ(𝐼−𝐺)𝑢₀,𝑢−𝑐_𝑖=𝐺_ℏ𝑢_𝑖−1,𝑖=2,3,…, (14) where 𝑢₀ is an initial guess of exact solution 𝑢. Therefore, 𝑣=𝑢₀+^∞_𝑖=1𝑝^𝑖𝐺_ℏ^𝑖−1𝑢₁. (15)

By setting 𝑝=1, we obtain 𝑣=𝑢₀+^∞_𝑖=1𝐺_ℏ^𝑖−1𝑢₁. (16)

It is obvious that if 𝜌(𝐺_ℏ)<1, then the series, ∑^∞_𝑖=1𝐺_ℏ^𝑖−1𝑢₁, converges and we have 𝑣=𝑢₀+𝐼−𝐺_ℏ⁻¹𝑢₁=𝑢₀+ℏ⁻¹(𝐼−𝐺)⁻¹𝑢₁=(𝐼−𝐺)⁻¹𝑐, (17) which is the exact solution of (3). A series of vectors can be computed by (14), and our aim is to choose the convergence control parameter ℏ≠0 so that 𝜌(𝐺_ℏ)<1. For improving the rate of convergence of iterative method, we present the following theorem.

Theorem 1. Suppose that 𝜌(𝐺)<1, and let 𝜆_𝑖=Re(𝜆_𝑖)+𝑖Im(𝜆_𝑖), and 𝜇_𝑖, 𝑖=1,2,…,𝑛 be the eigenvalues of 𝐺 and 𝐺_ℏ, respectively. Let 𝛼_𝑖=(1−Re(𝜆_𝑖))²+Im(𝜆_𝑖)², 𝛽_𝑖=2(1−Re(𝜆_𝑖)), and let 𝑔_𝑖(ℏ)=|𝜇_𝑖|²−𝜌(𝐺)²=𝛼_𝑖ℏ²+𝛽_𝑖ℏ+1−𝜌(𝐺)², 𝑖=1,2,…,𝑛. If 𝛼_𝑖≠0 and 𝛽_𝑖−2𝛼_𝑖≠0, 𝑖=1,2,…,𝑛, then (i)the quadratic equation 𝑔_𝑖(ℏ)=0, 𝑖=1,2,…,𝑛, has simple real roots 𝛾₁^(𝑖)<𝛾₂^(𝑖)<0, (ii)ℏ=−1 belongs to the interval [max^𝑛_𝑖=1𝛾₁^(𝑖),min^𝑛_𝑖=1𝛾₂^(𝑖)] and 𝜌(𝐺₋₁)=𝜌(𝐺)<1, (iii)for each ℏ∈[max^𝑛_𝑖=1𝛾₁^(𝑖),min^𝑛_𝑖=1𝛾₂^(𝑖)] and ℏ≠−1, the relation 𝜌(𝐺_ℏ)<𝜌(𝐺)<1 holds.

Proof. (i) We begin by defining two index sets 𝑁₁={𝑖∣|𝜆_𝑖|=𝜌(𝐺)} and 𝑁₂={𝑖∣|𝜆_𝑖|≠𝜌(𝐺)}. Since 𝛼_𝑖>0 and 𝜌(𝐺)<1, for 𝑖∈𝑁₂, we have sign𝑔_𝑖(−∞)=1,𝑔_𝑖(−1)<0,𝑔_𝑖(0)>0. (18)
So, 𝑔_𝑖(ℏ), 𝑖∈𝑁₂, has simple real roots 𝛾₁^(𝑖)<−1,−1<𝛾₂^(𝑖)<0. (19)
For 𝑖∈𝑁₁, we have 𝛼_𝑖−𝛽_𝑖<0 and 𝑔_𝑖(−1)=0,𝑔_𝑖𝛼_𝑖−𝛽_𝑖𝛼_𝑖=0. (20)
So, by using the assumption 𝛽_𝑖−2𝛼_𝑖≠0, 𝑔_𝑖(ℏ), 𝑖∈𝑁₂, has also simple real roots 𝛾₁^(𝑖)<𝛾₂^(𝑖)<0, defined as follows: 𝛾₁^(𝑖)=−1,𝛾₂^(𝑖)=𝛼_𝑖−𝛽_𝑖𝛼_𝑖𝛼if^𝑖−𝛽_𝑖𝛼_𝑖𝛾>−1,_1(𝑖)=𝛼_𝑖−𝛽_𝑖𝛼_𝑖,𝛾₂^(𝑖)𝛼=−1if^𝑖−𝛽_𝑖𝛼_𝑖<−1. (21)
(ii) From part (i), we observe that −1∈[max^𝑛_𝑖=1𝛾₁^(𝑖),min^𝑛_𝑖=1𝛾₂^(𝑖)] and 𝑔_𝑖(−1)=0. This implies that 𝜌(𝐺₋₁)=𝜌(𝐺)<1.
(iii) From part (i), we also observe that for each ℏ∈[max^𝑛_𝑖=1𝛾₁^(𝑖),min^𝑛_𝑖=1𝛾₂^(𝑖)] and ℏ≠−1. We have 𝑔_𝑖(ℏ)<0, and the relation 𝜌(𝐺_ℏ)<𝜌(𝐺)<1 holds.

The following theorem shows that by applying the HAM to a iterative scheme which is divergent, it is possible to construct a convergent homotopy-series vectors when the iteration matrix 𝐺 has particular properties.

Theorem 2. Let 𝜆_𝑖=Re(𝜆_𝑖)+𝑖Im(𝜆_𝑖), and 𝜇_𝑖, 𝑖=1,2,…,𝑛, be the eigenvalues of 𝐺 and 𝐺_ℏ, respectively. Let 𝛼_𝑖=(1−Re(𝜆_𝑖))²+Im(𝜆_𝑖)² and 𝛽_𝑖=2(1−Re(𝜆_𝑖)). Suppose that 𝛼_𝑖≠0, 𝑖=1,2,…,𝑛. (i)If  Re(𝜆_𝑖)>1,  for  𝑖=1,2,…,𝑛, and ℏ∈(0,min_{1≤𝑖≤𝑛}−(𝛽_𝑖/𝛼_𝑖)), then 𝜌(𝐺_ℏ)<1.(ii)If  Re(𝜆_𝑖)<1,   for  𝑖=1,2,…,𝑛, and ℏ∈(max_{1≤𝑖≤𝑛}−(𝛽_𝑖/𝛼_𝑖),0), then 𝜌(𝐺_ℏ)<1.

Proof. From (13), we have 𝜇_𝑖=1+ℏ(1−Re(𝜆_𝑖))−𝑖ℏIm(𝜆_𝑖). So, it is sufficient to have 𝜇_𝑖²−1=𝛼_𝑖ℏ²+𝛽_𝑖ℏ<0,for𝑖=1,2,…,𝑛. (22)
Under the hypothesis of part (i), we have 𝛽_𝑖<0, 𝑖=1,2,…,𝑛, and the relation (22) holds if ℏ∈(0,min_{1≤𝑖≤𝑛}−(𝛽_𝑖/𝛼_𝑖)), and the proof of (i) is complete. A similar argument holds if Re(𝜆_𝑖)<1,  for 𝑖=1,2,…,𝑛, and part (ii) follows.

When the assumption of part (i) (or part (ii)) of Theorem 2 does not hold, the following theorem shows that, for certain cases, instead of (3) we can consider the equivalent equation 𝑢=𝐺𝑢+𝑐,with𝐺=(𝐺−𝐼)²+𝐼,𝑐=(𝐺−𝐼)𝑐, (23) in which the iteration matrix 𝐺 has the eigenvalues with the desired properties.

Theorem 3. Let 𝜆_𝑖=Re(𝜆_𝑖)+𝑖Im(𝜆_𝑖), and 𝜈_𝑖=Re(𝜈_𝑖)+𝑖Im(𝜈_𝑖), 𝑖=1,2,…,𝑛, be the eigenvalues of 𝐺 and 𝐺, respectively.  If (1−Re(𝜆_𝑖))²−(𝐼𝑚(𝜆_𝑖))²>0 (or (1−Re(𝜆_𝑖))²−(Im(𝜆_𝑖))²<0) for 𝑖=1,2,…,𝑛, then Re(𝜈_𝑖)>1 (or Re(𝜈_𝑖)<1) for 𝑖=1,2,…,𝑛.

Proof. The proof immediately follows from the fact that Re(𝜈_𝑖)=((1−Re(𝜆_𝑖))²−(Im(𝜆_𝑖))²)+1.

The following corollary shows that by using the modified linear equation (23) and the homotopy analysis method with the corresponding 𝐺_ℏ=(ℏ+1)𝐼−ℏ𝐺, we can construct a convergent homotopy series for linear system (1).

Corollary 4. If 𝐺 has only real eigenvalues, then there exists ℏ≠0 such that the series of vectors generated by 𝑢₁=ℏ𝐼−𝐺𝑢₀−𝑐,𝑢_𝑖=𝐺_ℏ𝑢_𝑖−1,𝑖=2,3,…, (24) converges to the exact solution of (1).

Proof . The proof immediately follows from Theorems 2 and 3.

This corollary establishes that the series of vectors generated by (24) always converges if the iteration matrix 𝐺 is a symmetric matrix. When 𝐴 is symmetric with diagonal elements positive real numbers, (1) can be written as follows: (𝐷^−1/2𝐴𝐷^−1/2)(𝐷^1/2𝑥)=(𝐷^−1/2𝑏), where 𝐷 is the diagonal of 𝐴. Denoting again by 𝐴, 𝑥, and 𝑏 the expressions (𝐷^−1/2𝐴𝐷^−1/2), (𝐷^1/2𝑥), and (𝐷^−1/2𝑏), respectively, it is obvious the new coefficient 𝐴 is still symmetric and can therefore be written in the form 𝐴=𝐼−𝐿−𝐿^𝑇. An immediate consequence of Corollary 4 and the above discussion is the following results. (i) The series of vectors generated by (24) converges when 𝐴 is a symmetric matrix and the iterative method is the Richardson method 𝑢^(𝑘+1)=(𝐼−𝐴)𝑢^(𝑘)+𝑏. (25)(ii) The series of vectors generated by (24) converges when 𝐴=𝐼−𝐿−𝑈 is a symmetric matrix and the iterative method is the Jacobi method 𝑢^(𝑘+1)=(𝐿+𝑈)𝑢^(𝑘)+𝑏. (26)(iii) If 𝐴=𝐼−𝐿−𝑈 is a symmetric matrix and the iterative method is SAOR method 𝑢^(𝑘+1)=𝐽_𝑟,𝑤𝑢^(𝑘)+𝑐 (27) with 𝐽_𝑟,𝑤=(𝐼−𝑟𝑈)⁻¹)()]×()[(1−𝜔𝐼+𝜔−𝑟𝑈+𝜔𝐿𝐼−𝑟𝐿⁻¹)()],()[(1−𝜔𝐼+𝜔−𝑟𝐿+𝜔𝑈𝑐=𝜔𝐼−𝑟𝑈⁻¹[(2−𝜔)𝐼−(𝜔−𝑟)(𝐿+𝑈)](𝐼−𝑟𝐿)⁻¹𝑏, (28) then the series of vectors generated by (24) converges if 𝜔∈(0,2),𝜔+2−𝜔𝜇_min<𝑟<𝜔+2−𝜔𝜇_max, (29) where 𝜇_min and 𝜇_max stand for the minimum and the maximum eigenvalues of 𝐵=𝐿+𝑈, respectively. This result follows from the fact that in this case the iteration matrix 𝐽_𝑟,𝑤 of SAOR method has real eigenvalues (see Theorem 2 in [19]). (iv) If 𝐴=𝐼−𝐿−𝑈 is a symmetric matrix and the iterative method is SSOR method, then the series of vectors generated by (24) converges if 𝜔∈(0,2). This result immediately follows from the fact that the SAOR method reduces to the SSOR method for 𝑟=𝜔.

3. Numerical Examples

For numerical comparison, we use some matrices from the University of Florida sparse matrix Collection [20]. These matrices with their properties are shown in Table 1. We determined the spectral radii of iteration matrices of the classical SOR, AOR, SSOR, SAOR, Jacobi, and Richardson methods (𝜌(𝐺)) as well as those of the corresponding 𝐺_ℏ after the application HAM method with the experimentally computed optimal value of ℏ(𝜌(𝐺_ℏopt)). In Tables 2–4, we list 𝜌(𝐺), 𝜌(𝐺_ℏopt), the interval of convergence which introduced in Theorems 1 and 2, the experimentally computed optimal value of ℏ(ℏ_opt), and the spectral radius of iteration matrix 𝜌(𝐺_ℏopt) which introduced in Corollary 4.

Table 1 

Test problem information.

Table 2 

The basic method is convergent (Theorem 1).

Table 3 

The basic method is divergent (Theorem 2).

Table 4 

The basic method is divergent (Theorem 3).

In Table 2, we consider the convergent classical methods (𝜌(𝐺)<1). It is easy to verify that the numerical results are consistent with Theorem 1, and we observe that for the convergent classical methods by choosing suitable convergence control parameter ℏ, the rate of the convergence of the HAM method is faster than that of corresponding classical method.

In Table 3, we consider the divergent classical methods when Re(𝜆_𝑖)>1 (or Re(𝜆_𝑖)<1) for 𝑖=1,2,…,𝑛. We can see that the numerical results are consistent with Theorem 2. These results show that by applying the HAM to a iterative scheme which is divergent, it is possible to construct a convergent homotopy-series vectors when the iteration matrix 𝐺 has mentioned properties.

In Table 4, we report the results obtained for the symmetric matrix Si2 which has positive diagonal elements. For this example, the classical methods diverge and there exist 𝑖,𝑗 such that Re(𝜆_𝑖)<1 and Re(𝜆_𝑗)>1. We observe that the results are consistent with Theorem 3 and Corollary 4. The results show that the HAM method is convergent but the rate of the convergence is slow.

Finally, Tables 3 and 4 show that it is not necessary to choose the the parameters 𝑟 and 𝜔 in the convergence interval of the classical methods. In the case of divergence, under the assumptions of Theorems 2 and 3, the application of the HAM can generate the convergent homotopy-series vectors for linear system (1).

4. Conclusion

In this paper, we proposed to apply the homotopy analysis method to the classical iterative methods for solving the linear system of equations. The theoretical results show the HAM can be used to accelerate the convergence of the basic iterative methods. In addition, we show that by applying the HAM to a divergent iterative scheme, it is possible to construct a convergent homotopy-series solution when the iteration matrix 𝐺 of the iterative scheme has particular properties. The numerical experiments confirm the theoretical results and show the efficiency of the new method.