- Email: [email protected]
Application of homotopy perturbation methods for solving systems of linear equations
Applied Mathematics and Computation 217 (2011) 5259–5264
Contents lists available at ScienceDirect
Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc
Application of homotopy perturbation methods for solving systems of linear equations Hsuan-Ku Liu Department of Mathematics and Information Education, National Taipei University of Education, Taiwan
a r t i c l e
i n f o
a b s t r a c t In this paper, homotopy perturbation methods (HPMs) are applied to obtain the solution of linear systems, and conditions are deduced to check the convergence of the homotopy series. Moreover, we have adapted the Richardson method, the Jacobi method, and the Gauss– Seidel method to choose the splitting matrix. The numerical results indicate that the homotopy series converges much more rapidly than the direct methods for large sparse linear systems with a small spectrum radius. Ó 2010 Elsevier Inc. All rights reserved.
Keywords: HPM Linear systems Richardson method Jacobi method Gauss–Seidel method
1. Introduction Many numerical methods, such as iterative methods [5] and HPMs [4], are suggested to search for the solution of linear systems. Keramati [4] applied a HPM to solve linear systems. The splitting matrix of this method is only the identity matrix. However, this method does not converge for some systems when the spectrum radius is greater than 1. To make the method available, the auxiliary parameter and the auxiliary matrix are added to the homotopy method; in other words, the HPMs are applied to consider the linear systems. Nowadays, HPMs are utilized to determine the solution of many problems, such as nonlinear heat transfer equations [2], nonlinear wave equations [2,3,4] and van der Pol equations [6]. In addition to HPMs, many analytic techniques have been proposed using homotopy, a fundamental concept of topology. The idea of homotopy is a continuous map from the interval [0, 1] to a function space. To solve a nonlinear equation f(u) = 0, the continuation procedure is done to track the set of zeros of
Hðu; pÞ ¼ ð1 pÞf ðuÞ þ pgðuÞ ¼ 0; where g(u) = 0 is a starting problem whose solution u0 is known. Watson [7] pointed out that continuation can fail because the curve of zeros of H(u, p) emanating from (u0, 0) may (i) have turning points, (ii) bifurcate, (iii) fail to exist at some p values, or (iv) wander off to infinity without reaching p = 1. In 1976, Chow et al. [1] proposed the probability-one homotopies that was able to overcome all the four shortcomings of the continuation and homotopy methods. In this study, the homotopy was always continuous for linear systems. The goal of this paper was to propose a new iterative method for solving large sparse linear systems
Ax ¼ b;
ð1Þ
where
A ¼ ½aij ;
x ¼ ½xj ;
b ¼ ½bj ;
i ¼ 1; 2; . . . ; n; j ¼ 1; 2; . . . ; n:
E-mail address: [email protected] 0096-3003/$ - see front matter Ó 2010 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2010.11.024
5260
H.-K. Liu / Applied Mathematics and Computation 217 (2011) 5259–5264
Instead of using the traditional homotopy, HPMs introduce a nonzero auxiliary parameter and a nonzero auxiliary matrix. Hence, the advantage of HPM is the freedom to choose the auxiliary parameter, where the auxiliary parameters may play the same role as overrelaxation parameters in iteration methods [5]. In summary, the main contributions of this paper are (1) to give an approximation for (1) by applying HPMs, (2) to derive conditions for checking the convergence of the homotopy series, (3) to adapt the Richardson method, Jacobi method and Gauss–Seidel method for choosing the splitting matrix. The numerical studies show that the homotopy series converges rapidly when the spectrum radius of the linear system is small. This paper is organized as follows. In Section 2, we have introduced the basic concept of HPM, and have derived the conditions for the convergence of the homotopy series. In Section 3, we have introduced the homotopy Richardson method, homotopy Jacobi method, and homotopy Gauss–Seidel method. In Section 4, we have presented a numerical example which is a large sparse linear system. The concise conclusion is provided in Section 5. 2. Basic concept in HPMs We have considered the HPMs in a general mathematical setting. A general type of homotopy method for solving (1) can be described as follows. By applying a certain matrix Q, we can define homotopy H(u, p) by
Hðu; 0Þ ¼ FðuÞ and Hðu; 1Þ ¼ LðuÞ and a convex homotopy as follows
Hðu; pÞ ¼ ð1 pÞFðuÞ þ ðqHÞpLðuÞ ¼ 0;
ð2Þ
where
LðuÞ ¼ Au b and FðuÞ ¼ Qu w0 : Here, we have the freedom to choose the auxiliary parameter q, the auxiliary matrix H, the initial approximation w0, and the auxiliary operator F(u). Note that the operator F(u) is decided by the splitting matrix Q. HPM uses the homotopy parameter p as an expanding parameter to obtain
u ¼ u0 þ u1 p þ u2 p2 þ
ð3Þ
and it gives an approximation to the solution of (1) as
v ¼ lim ðu0 þ u1 p þ u2 p2 þ Þ: p!1 By substituting (3) in (2), and by equating the terms with the identical power of p, we can obtain
p0 : Qu0 w0 ¼ 0 p1 : Qu1 þ ðqHA QÞu0 þ w0 qHb ¼ 0 pi : Qui þ ðqHA Q Þui1 ¼ 0; i ¼ 2; 3; . . . This implies
u0 ¼ Q 1 w0 u1 ¼ ðI qQ 1 HAÞu0 þ Q 1 ðqHb w0 Þ ui ¼ ðI qQ 1 HAÞui1 ; i ¼ 2; 3; . . . Taking w0 = qHb yields
u0 ¼ qðQ 1 HÞb ui ¼ ðI qQ 1 HAÞi ðqQ 1 HÞb; i ¼ 1; 2; 3; . . .
ð4Þ
Therefore,
u¼
1 X i¼0
ui pi ¼
1 h X i¼0
i ðI qQ 1 HAÞi ðqQ 1 HÞb pi :
ð5Þ
H.-K. Liu / Applied Mathematics and Computation 217 (2011) 5259–5264
5261
By setting p = 1, the solution can be of the form
v¼
1 X
ui ¼
i¼0
1 X ðI qQ 1 HAÞi ðqQ 1 HÞb: i¼0
A series of vectors can be computed from (4), and our objective is to choose the splitting matrix Q, the auxiliary parameter q and the auxiliary matrix H so that (1) the sequence {ui} is computed easily, (2) the series u in (5) converges rapidly to the solution. In our analysis, we tried to determine whether the sequence {ui} can be easily computed if the system
Qxi ¼ ðQ qHAÞxi1 can be solved easily and quickly. Hence the splitting matrix was chosen so that the above mentioned equation can be solved easily and quickly. To verify whether the series u converges, we proposed the following theorems. Let q(S) denote the spectral radius, that is,
qðSÞ ¼ max jkj; k2rðSÞ
where r(S) denotes the set of all eigenvalues of S and is called the spectrum of S. P k Saad [5] showed that the series 1 k¼0 S converges if and only if q(S) < 1. Consequently, we obtained the following results. Theorem 2.1. The series
u¼
1 X ðI qQ 1 HAÞi ðqQ 1 HÞb i¼0
converges if and only if
qðI qQ 1 HAÞ < 1: Moreover, the following corollary holds for qðI qQ 1 HAÞ 6 kI qQ 1 HAk for any matrix norm. Corollary 2.2. The series
u¼
1 X ðI qQ 1 HAÞi ðqQ 1 HÞb i¼0
converges if
kI qQ 1 HAk < 1: Definition 2.3. Let A, M, N be the three given matrices satisfying A = M N. The pair of matrices M, N is a regular splitting of A if M is nonsingular and M1 and N are nonnegative. Lemma 2.4 [5]. Let M, N be a regular splitting of a matrix A. Then q(M1N) < 1 if and only if A is nonsingular and A1 is nonnegative. By using Lemma 2.4 we can obtain the following results. Theorem 2.5. Let Q be a nonsingular matrix such that Q1 and Q qHA are nonnegative. Then the sequence
u¼
1 X ðI qQ 1 HAÞi ðqQ 1 HÞb i¼0
converges if qHA is nonsingular and (qHA)1 is nonnegative. Proof. Suppose that Q is a nonsingular matrix such that Q1 and Q qHA are nonnegative. By using Lemma 2.4 we can obtain q(I qQ1HA) < 1 as both Q and Q qHA are a regular splitting of qHA. This implies that
u½m ¼
m X ðI qQ 1 HAÞi ðqQ 1 HÞb i¼0
converges by using Theorem 2.1. h
5262
H.-K. Liu / Applied Mathematics and Computation 217 (2011) 5259–5264
In summary, the auxiliary parameter q and the auxiliary matrix H are properly chosen such that the power series (5) converges at p = 1. In other words, we choose the auxiliary parameter and the auxiliary matrix such that one of the following three conditions holds: (1) q(I qQ1HA) < 1; (2) kI qQ 1 HAk < 1; (3) The matrices Q and qHA are nonsingular, and the matrices Q1, (Q qHA) and (qHA)1 are nonnegative. 3. Homotopy perturbation methods Following the idea of iteration methods, a considerable number of iteration methods could be adapted to the homotopy iteration approximations, such as the Richardson method, Jacobi method and Gauss–Seidel method. In this section, we shall introduce these approximations and call them homotopy Richardson method, homotopy Jacobi method and homotopy Gauss–Seidel method. 3.1. Homotopy Richardson method Let QRich = In, we have
u0 ¼ qHb; ui ¼ ðI qHAÞi u0 : It can be noted that if q = 1 and H = In, then we get the same results as proposed by Karamati [4]. 3.2. Homotopy Jacobi method Another illustration is provided by Jacobi approximation, in which Q is the diagonal matrix whose diagonal entries are the same as those in the matrix A, namely QJac = D = diag (aii). Using (4), we can obtain
u0 ¼ qD1 Hb; ui ¼ ðI qD1 HAÞi u0 : 3.3. Homotopy Gauss–Seidel method The Gauss–Seidel method is defined by letting Q to be the lower triangular part of S, including its diagonal. Hence, we have
0
Q GS
a11 B a B 21 B B . ¼ D þ L ¼ B .. B B @ an1;1 an;1
0 a22 an1;2 an;2
0
an1;n1
an;n1
0
1
C 0 C C C C: C C 0 A an;n
Using (4), we have
Q GS u0 ¼ qHb; Q GS ui ¼ ðQ GS qHAÞi u0 :
ð6Þ
As Q is a lower triangular matrix, the solution of (6) can be computed easily and quickly. Thus, in this section, the three homotopy iterative approximations, the homotopy Richardson approximation, homotopy Jacobi approximation and homotopy Gauss–Seidel approximation, have been presented. 4. Comparison with the other methods In this section, we will compare HPMs with the direct methods and iteration methods [8]; in other words, the HPMs, the iteration methods, and the direct methods are applied to find the solution of the system
Ax ¼ b; where A is a 1000 1000 matrix with q(I A) = 0.4997 and b is a 1000 1 vector. The errors, computation times and spectrum radius for the direct methods and the HPMs are displayed in Table 1. As the real solution is not known, errors
5263
H.-K. Liu / Applied Mathematics and Computation 217 (2011) 5259–5264 Table 1 The direct methods and HPMs. Methods
Direct methods GJ
Spectrum Iterations MAD Time (s)
3.52e15 1.42e4
Homotopy perturbation methods LU
Jacobi
Richardson
Gauss-Seidel
3.44e15 1.19e4
0.3414 30 2.33e15 0.0345
0.4997 30 2.25e10 0.0215
0.6336 30 3.4e7 0.0869
50 2.44e15 0.0334
50 2.65e11 0.1375
Table 2 The iteration methods. Methods
Jacobi
Richardson
Gauss-Seidel
Iterations MAD Time (s)
30 1.44e15 0.0343
30 2.26e10 0.0206
30 2.53e7 0.42
are computed by the maximum absolute difference between Ax and b. As we can see the errors for the methods are all less than 1010. The computation times for the direct methods, such as Gauss–Jordan elimination and LU decomposition, are greater than 104 but the computation times for HPMs are less than 102. This implies that HPMs converge rapidly for large sparse systems with small spectrum radius. In this example, the spectrum radius for the Jacobi, Richardson, and Gauss–Seidel methods are 0.3414, 0.4997, and 0.6336, respectively. Thus, the Jacobi method converges to 1015 after 30 iterations but the Gauss–Seidel method converges only to 1011 after 50 iterations. This implies that the spectrum radius provides a method for selecting the splitting matrices. For the 30 iterations we also observed that the computation time of the Gauss–Seidel method is greater than that of the other two homotopy methods. This implies that the homotopy Gauss–Seidel method may not be a good method for solving large sparse systems. When HPMs were compared with iteration methods, we found that the errors and the computation times were very close (Table 2). However, we also observed some advantages of HPMs for solving a sequence of linear systems
Ax ¼ bi ; where fbi gni¼1 is a sequence of distinct vectors. We could compute the homotopy series
E¼
1 X
ðI Q 1 HAÞi ðQ 1 HÞ
ð7Þ
i¼0
and define
xi ¼ Ebi : When the change is only bi, we do not compute the homotopy series E again. Hence, HPMs can be considered better than the iterative methods for solving a sequence of linear systems. The multiplications in (7) are matrix multiplications, but the multiplications in (4) are vector multiplications; in other words the computation time for (7) may be greater than that for (4). 5. Conclusion This paper applied HPMs to solve the linear systems and derive conditions to check the convergence of the homotopy series. Moreover, the Richardson method, Jacobi method, and Gauss–Seidel method were adapted to choose the splitting matrix. From the numerical results, we observed that the homotopy series converged rapidly for a large sparse system with a small spectrum radius. To improve the rate of convergence, future studies should focus on a method for choosing an auxiliary parameter and an auxiliary matrix. References [1] S.N. Chow, J. Mallet-Paret, J.A. Yorke, Finding zeros of maps: homotopy methods that are constructive with probability one, Math. Comp. 32 (1978) 887– 899. [2] D.D. Ganji, A. Sadighi, Application of homotopy-perturbation and variational iteration methods to nonlinear heat transfer and porous media equations, Appl. Math. Comput. 27 (2007) 24–34. [3] J.-H. He, HPM: a new nonlinear analytical technique, Appl. Math. Comput. 135 (2003) 73–79. [4] B. Keramati, An approach to the solution of linear system of equations by He’s HPM, Chaos Solitons Fract. doi:10.1016/j.chaos.2007.11.020. [5] Y. Saad, Iterative Methods for Sparse Linear Systems, second ed., SIAM, 2003.
5264
H.-K. Liu / Applied Mathematics and Computation 217 (2011) 5259–5264
[6] T. Ozis, A. Yildirim, A note on He’s homotopy perturbation method for van der Pol oscillator with very strong nonlinearity, Chaos Solitons Fract. 34 (2007) 989–991. [7] L.T. Watson, Probability-one homotopies in computational science, J. Comput. Appl. Math. 140 (2002) 785–807. [8] W.Y. Yang, W. Cao, T.-S. Chung, J. Morries, Applied Numerical Methods Using Matlab, John Wiley & Sons, Inc., 2005.
Contents lists available at ScienceDirect
Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc
Application of homotopy perturbation methods for solving systems of linear equations Hsuan-Ku Liu Department of Mathematics and Information Education, National Taipei University of Education, Taiwan
a r t i c l e
i n f o
a b s t r a c t In this paper, homotopy perturbation methods (HPMs) are applied to obtain the solution of linear systems, and conditions are deduced to check the convergence of the homotopy series. Moreover, we have adapted the Richardson method, the Jacobi method, and the Gauss– Seidel method to choose the splitting matrix. The numerical results indicate that the homotopy series converges much more rapidly than the direct methods for large sparse linear systems with a small spectrum radius. Ó 2010 Elsevier Inc. All rights reserved.
Keywords: HPM Linear systems Richardson method Jacobi method Gauss–Seidel method
1. Introduction Many numerical methods, such as iterative methods [5] and HPMs [4], are suggested to search for the solution of linear systems. Keramati [4] applied a HPM to solve linear systems. The splitting matrix of this method is only the identity matrix. However, this method does not converge for some systems when the spectrum radius is greater than 1. To make the method available, the auxiliary parameter and the auxiliary matrix are added to the homotopy method; in other words, the HPMs are applied to consider the linear systems. Nowadays, HPMs are utilized to determine the solution of many problems, such as nonlinear heat transfer equations [2], nonlinear wave equations [2,3,4] and van der Pol equations [6]. In addition to HPMs, many analytic techniques have been proposed using homotopy, a fundamental concept of topology. The idea of homotopy is a continuous map from the interval [0, 1] to a function space. To solve a nonlinear equation f(u) = 0, the continuation procedure is done to track the set of zeros of
Hðu; pÞ ¼ ð1 pÞf ðuÞ þ pgðuÞ ¼ 0; where g(u) = 0 is a starting problem whose solution u0 is known. Watson [7] pointed out that continuation can fail because the curve of zeros of H(u, p) emanating from (u0, 0) may (i) have turning points, (ii) bifurcate, (iii) fail to exist at some p values, or (iv) wander off to infinity without reaching p = 1. In 1976, Chow et al. [1] proposed the probability-one homotopies that was able to overcome all the four shortcomings of the continuation and homotopy methods. In this study, the homotopy was always continuous for linear systems. The goal of this paper was to propose a new iterative method for solving large sparse linear systems
Ax ¼ b;
ð1Þ
where
A ¼ ½aij ;
x ¼ ½xj ;
b ¼ ½bj ;
i ¼ 1; 2; . . . ; n; j ¼ 1; 2; . . . ; n:
E-mail address: [email protected] 0096-3003/$ - see front matter Ó 2010 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2010.11.024
5260
H.-K. Liu / Applied Mathematics and Computation 217 (2011) 5259–5264
Instead of using the traditional homotopy, HPMs introduce a nonzero auxiliary parameter and a nonzero auxiliary matrix. Hence, the advantage of HPM is the freedom to choose the auxiliary parameter, where the auxiliary parameters may play the same role as overrelaxation parameters in iteration methods [5]. In summary, the main contributions of this paper are (1) to give an approximation for (1) by applying HPMs, (2) to derive conditions for checking the convergence of the homotopy series, (3) to adapt the Richardson method, Jacobi method and Gauss–Seidel method for choosing the splitting matrix. The numerical studies show that the homotopy series converges rapidly when the spectrum radius of the linear system is small. This paper is organized as follows. In Section 2, we have introduced the basic concept of HPM, and have derived the conditions for the convergence of the homotopy series. In Section 3, we have introduced the homotopy Richardson method, homotopy Jacobi method, and homotopy Gauss–Seidel method. In Section 4, we have presented a numerical example which is a large sparse linear system. The concise conclusion is provided in Section 5. 2. Basic concept in HPMs We have considered the HPMs in a general mathematical setting. A general type of homotopy method for solving (1) can be described as follows. By applying a certain matrix Q, we can define homotopy H(u, p) by
Hðu; 0Þ ¼ FðuÞ and Hðu; 1Þ ¼ LðuÞ and a convex homotopy as follows
Hðu; pÞ ¼ ð1 pÞFðuÞ þ ðqHÞpLðuÞ ¼ 0;
ð2Þ
where
LðuÞ ¼ Au b and FðuÞ ¼ Qu w0 : Here, we have the freedom to choose the auxiliary parameter q, the auxiliary matrix H, the initial approximation w0, and the auxiliary operator F(u). Note that the operator F(u) is decided by the splitting matrix Q. HPM uses the homotopy parameter p as an expanding parameter to obtain
u ¼ u0 þ u1 p þ u2 p2 þ
ð3Þ
and it gives an approximation to the solution of (1) as
v ¼ lim ðu0 þ u1 p þ u2 p2 þ Þ: p!1 By substituting (3) in (2), and by equating the terms with the identical power of p, we can obtain
p0 : Qu0 w0 ¼ 0 p1 : Qu1 þ ðqHA QÞu0 þ w0 qHb ¼ 0 pi : Qui þ ðqHA Q Þui1 ¼ 0; i ¼ 2; 3; . . . This implies
u0 ¼ Q 1 w0 u1 ¼ ðI qQ 1 HAÞu0 þ Q 1 ðqHb w0 Þ ui ¼ ðI qQ 1 HAÞui1 ; i ¼ 2; 3; . . . Taking w0 = qHb yields
u0 ¼ qðQ 1 HÞb ui ¼ ðI qQ 1 HAÞi ðqQ 1 HÞb; i ¼ 1; 2; 3; . . .
ð4Þ
Therefore,
u¼
1 X i¼0
ui pi ¼
1 h X i¼0
i ðI qQ 1 HAÞi ðqQ 1 HÞb pi :
ð5Þ
H.-K. Liu / Applied Mathematics and Computation 217 (2011) 5259–5264
5261
By setting p = 1, the solution can be of the form
v¼
1 X
ui ¼
i¼0
1 X ðI qQ 1 HAÞi ðqQ 1 HÞb: i¼0
A series of vectors can be computed from (4), and our objective is to choose the splitting matrix Q, the auxiliary parameter q and the auxiliary matrix H so that (1) the sequence {ui} is computed easily, (2) the series u in (5) converges rapidly to the solution. In our analysis, we tried to determine whether the sequence {ui} can be easily computed if the system
Qxi ¼ ðQ qHAÞxi1 can be solved easily and quickly. Hence the splitting matrix was chosen so that the above mentioned equation can be solved easily and quickly. To verify whether the series u converges, we proposed the following theorems. Let q(S) denote the spectral radius, that is,
qðSÞ ¼ max jkj; k2rðSÞ
where r(S) denotes the set of all eigenvalues of S and is called the spectrum of S. P k Saad [5] showed that the series 1 k¼0 S converges if and only if q(S) < 1. Consequently, we obtained the following results. Theorem 2.1. The series
u¼
1 X ðI qQ 1 HAÞi ðqQ 1 HÞb i¼0
converges if and only if
qðI qQ 1 HAÞ < 1: Moreover, the following corollary holds for qðI qQ 1 HAÞ 6 kI qQ 1 HAk for any matrix norm. Corollary 2.2. The series
u¼
1 X ðI qQ 1 HAÞi ðqQ 1 HÞb i¼0
converges if
kI qQ 1 HAk < 1: Definition 2.3. Let A, M, N be the three given matrices satisfying A = M N. The pair of matrices M, N is a regular splitting of A if M is nonsingular and M1 and N are nonnegative. Lemma 2.4 [5]. Let M, N be a regular splitting of a matrix A. Then q(M1N) < 1 if and only if A is nonsingular and A1 is nonnegative. By using Lemma 2.4 we can obtain the following results. Theorem 2.5. Let Q be a nonsingular matrix such that Q1 and Q qHA are nonnegative. Then the sequence
u¼
1 X ðI qQ 1 HAÞi ðqQ 1 HÞb i¼0
converges if qHA is nonsingular and (qHA)1 is nonnegative. Proof. Suppose that Q is a nonsingular matrix such that Q1 and Q qHA are nonnegative. By using Lemma 2.4 we can obtain q(I qQ1HA) < 1 as both Q and Q qHA are a regular splitting of qHA. This implies that
u½m ¼
m X ðI qQ 1 HAÞi ðqQ 1 HÞb i¼0
converges by using Theorem 2.1. h
5262
H.-K. Liu / Applied Mathematics and Computation 217 (2011) 5259–5264
In summary, the auxiliary parameter q and the auxiliary matrix H are properly chosen such that the power series (5) converges at p = 1. In other words, we choose the auxiliary parameter and the auxiliary matrix such that one of the following three conditions holds: (1) q(I qQ1HA) < 1; (2) kI qQ 1 HAk < 1; (3) The matrices Q and qHA are nonsingular, and the matrices Q1, (Q qHA) and (qHA)1 are nonnegative. 3. Homotopy perturbation methods Following the idea of iteration methods, a considerable number of iteration methods could be adapted to the homotopy iteration approximations, such as the Richardson method, Jacobi method and Gauss–Seidel method. In this section, we shall introduce these approximations and call them homotopy Richardson method, homotopy Jacobi method and homotopy Gauss–Seidel method. 3.1. Homotopy Richardson method Let QRich = In, we have
u0 ¼ qHb; ui ¼ ðI qHAÞi u0 : It can be noted that if q = 1 and H = In, then we get the same results as proposed by Karamati [4]. 3.2. Homotopy Jacobi method Another illustration is provided by Jacobi approximation, in which Q is the diagonal matrix whose diagonal entries are the same as those in the matrix A, namely QJac = D = diag (aii). Using (4), we can obtain
u0 ¼ qD1 Hb; ui ¼ ðI qD1 HAÞi u0 : 3.3. Homotopy Gauss–Seidel method The Gauss–Seidel method is defined by letting Q to be the lower triangular part of S, including its diagonal. Hence, we have
0
Q GS
a11 B a B 21 B B . ¼ D þ L ¼ B .. B B @ an1;1 an;1
0 a22 an1;2 an;2
0
an1;n1
an;n1
0
1
C 0 C C C C: C C 0 A an;n
Using (4), we have
Q GS u0 ¼ qHb; Q GS ui ¼ ðQ GS qHAÞi u0 :
ð6Þ
As Q is a lower triangular matrix, the solution of (6) can be computed easily and quickly. Thus, in this section, the three homotopy iterative approximations, the homotopy Richardson approximation, homotopy Jacobi approximation and homotopy Gauss–Seidel approximation, have been presented. 4. Comparison with the other methods In this section, we will compare HPMs with the direct methods and iteration methods [8]; in other words, the HPMs, the iteration methods, and the direct methods are applied to find the solution of the system
Ax ¼ b; where A is a 1000 1000 matrix with q(I A) = 0.4997 and b is a 1000 1 vector. The errors, computation times and spectrum radius for the direct methods and the HPMs are displayed in Table 1. As the real solution is not known, errors
5263
H.-K. Liu / Applied Mathematics and Computation 217 (2011) 5259–5264 Table 1 The direct methods and HPMs. Methods
Direct methods GJ
Spectrum Iterations MAD Time (s)
3.52e15 1.42e4
Homotopy perturbation methods LU
Jacobi
Richardson
Gauss-Seidel
3.44e15 1.19e4
0.3414 30 2.33e15 0.0345
0.4997 30 2.25e10 0.0215
0.6336 30 3.4e7 0.0869
50 2.44e15 0.0334
50 2.65e11 0.1375
Table 2 The iteration methods. Methods
Jacobi
Richardson
Gauss-Seidel
Iterations MAD Time (s)
30 1.44e15 0.0343
30 2.26e10 0.0206
30 2.53e7 0.42
are computed by the maximum absolute difference between Ax and b. As we can see the errors for the methods are all less than 1010. The computation times for the direct methods, such as Gauss–Jordan elimination and LU decomposition, are greater than 104 but the computation times for HPMs are less than 102. This implies that HPMs converge rapidly for large sparse systems with small spectrum radius. In this example, the spectrum radius for the Jacobi, Richardson, and Gauss–Seidel methods are 0.3414, 0.4997, and 0.6336, respectively. Thus, the Jacobi method converges to 1015 after 30 iterations but the Gauss–Seidel method converges only to 1011 after 50 iterations. This implies that the spectrum radius provides a method for selecting the splitting matrices. For the 30 iterations we also observed that the computation time of the Gauss–Seidel method is greater than that of the other two homotopy methods. This implies that the homotopy Gauss–Seidel method may not be a good method for solving large sparse systems. When HPMs were compared with iteration methods, we found that the errors and the computation times were very close (Table 2). However, we also observed some advantages of HPMs for solving a sequence of linear systems
Ax ¼ bi ; where fbi gni¼1 is a sequence of distinct vectors. We could compute the homotopy series
E¼
1 X
ðI Q 1 HAÞi ðQ 1 HÞ
ð7Þ
i¼0
and define
xi ¼ Ebi : When the change is only bi, we do not compute the homotopy series E again. Hence, HPMs can be considered better than the iterative methods for solving a sequence of linear systems. The multiplications in (7) are matrix multiplications, but the multiplications in (4) are vector multiplications; in other words the computation time for (7) may be greater than that for (4). 5. Conclusion This paper applied HPMs to solve the linear systems and derive conditions to check the convergence of the homotopy series. Moreover, the Richardson method, Jacobi method, and Gauss–Seidel method were adapted to choose the splitting matrix. From the numerical results, we observed that the homotopy series converged rapidly for a large sparse system with a small spectrum radius. To improve the rate of convergence, future studies should focus on a method for choosing an auxiliary parameter and an auxiliary matrix. References [1] S.N. Chow, J. Mallet-Paret, J.A. Yorke, Finding zeros of maps: homotopy methods that are constructive with probability one, Math. Comp. 32 (1978) 887– 899. [2] D.D. Ganji, A. Sadighi, Application of homotopy-perturbation and variational iteration methods to nonlinear heat transfer and porous media equations, Appl. Math. Comput. 27 (2007) 24–34. [3] J.-H. He, HPM: a new nonlinear analytical technique, Appl. Math. Comput. 135 (2003) 73–79. [4] B. Keramati, An approach to the solution of linear system of equations by He’s HPM, Chaos Solitons Fract. doi:10.1016/j.chaos.2007.11.020. [5] Y. Saad, Iterative Methods for Sparse Linear Systems, second ed., SIAM, 2003.
5264
H.-K. Liu / Applied Mathematics and Computation 217 (2011) 5259–5264
[6] T. Ozis, A. Yildirim, A note on He’s homotopy perturbation method for van der Pol oscillator with very strong nonlinearity, Chaos Solitons Fract. 34 (2007) 989–991. [7] L.T. Watson, Probability-one homotopies in computational science, J. Comput. Appl. Math. 140 (2002) 785–807. [8] W.Y. Yang, W. Cao, T.-S. Chung, J. Morries, Applied Numerical Methods Using Matlab, John Wiley & Sons, Inc., 2005.