# On the Application of Homotopy Perturbation Method for Solving Systems of Linear Equations

Authors: S. A. Edalatpanah; M. M. Rashidi
Journal: International Scholarly Research Notices
 (2014)

---

## Abstract

The application of homotopy perturbation method (HPM) for solving systems of linear equations is further discussed and focused on a method for choosing an auxiliary matrix to improve the rate of convergence. Moreover, solving of convection-diffusion equations has been developed by HPM and the convergence properties of the proposed method have been analyzed in detail; the obtained results are compared with some other methods in the frame of HPM. Numerical experiment shows a good improvement on the convergence rate and the efficiency of this method.


---

## Body

## 1. Introduction

Consider the following convection-diffusion equation:(1)
-
(
u
x
x
+
u
y
y
)
︸
Δ
u
+
δ
u
x
+
τ
u
y
=
f
(
x
,
y
)
.
On the unit square domain Ω
=
0,1
×
0,1, with constant coefficients δ, τ, subject to Dirichlet boundary conditions. Discretization by a five-point finite difference operator leads to the following linear system:
(2)
A
X
=
b
,
where X denotes a vector in a finite-dimensional space and A
∈
R
n
2
×
n
2. With discretization on a uniform n
×
n grid, using standard second-order differences for the Laplacian and either centered or upwind differences for the first derivatives, the coefficient matrix has the form
(3)
A
=
tridiagonal
[
b
I
,
tridiagonal
[
c
,
a
,
d
]
,
e
I
]
,
m
m
m
c
b
=
-
1
+
∂
;
c
=
-
1
+
γ
;
a
=
4
;
m
m
m
c
m
m
c
d
=
-
1
-
γ
;
e
=
-
1
-
∂
.
Or with Kronecker product (⊗) we obtain
(4)
A
=
I
⊗
T
x
+
I
⊗
T
y
,
T
x
=
tridiagonal
[
c
,
a
,
d
]
,
T
y
=
tridiagonal
b
I
,
0
,
e
I
,
where ∂
=
τ
h
/
2; γ
=
δ
h
/
2 are Reynolds numbers. Also, the equidistant step-size h
=
1
/
n is used in the discretization and the natural lexicographic ordering is employed to the unknowns and the right hand side satisfies b
i
j
=
h
2
f
i
j
(
x
,
y
). For details, we refer to [1, 2] and the references therein. Convection-diffusion equations appear in many fields of science and engineering and there are some reliable methods for solving this class of problems ([1, 2] and the references therein). Here, we use alternative approach to solve (2) based on homotopy perturbation method (HPM). The HPM was first proposed by the Chinese mathematician He in 1999 [3] and was further developed and improved by him [4–6]. He presented a homotopy perturbation technique based on the introduction of homotopy in topology coupled with the traditional perturbation method for the solution of algebraic equations and ODEs. This technique provides a summation of an infinite series with easily computable terms, which converges to the solution of the problem. In the literature, various authors have successfully applied this method for many kinds of different problems such as nonlinear partial differential equations [7, 8], nonlinear wave equations [5], nonlinear integral and integrodifferential equations [9], fractional IVPs [10], and optimization [11]. Nevertheless, to our knowledge, a few papers have considered this analytical method for system of linear equations. Keramati [12] presented an efficient method for solving system of linear equations based on HPM. Liu [13] by using HPM and splitting methods proposed a new method for linear systems. Saberi Najafi et al. [14] based on HPM proposed a new algorithm for fuzzy linear systems and compared it with Adomian’s decomposition method. They also in [15] show that solving linear systems by using a new method called modified HPM, presented by Noor et al. [16], is impractical.In this paper, by combination of HPM and preconditioning technique, we will present a new model for solving linear systems and apply this method to the above convection-diffusion equation.


## 2. Homotopy Perturbation Method for Linear Systems

Consider linear system (2). HPM procedure for solving of this system is a continuous map from the interval [
0,1
] to a function space where it is as follows:
(5)
H
(
u
,
p
)
=
(
1
-
p
)
F
(
u
)
+
p
L
(
u
)
=
0
,
and p
∈
[
0,1
] is an embedding parameter.Keramati [12] applied a HPM to solve linear systems by the following definitions.LetL
(
u
)
=
A
u
-
b and F
(
u
)
=
u
-
w
0, where w
0 is a known vector. Then homotopy H
(
u
,
p
) is
(6)
H
(
u
,
p
)
=
(
1
-
p
)
(
u
-
w
0
)
+
p
(
A
u
-
b
)
=
0
.
Obviously, we will have
(7)
H
(
u
,
0
)
=
F
(
u
)
,
H
(
u
,
1
)
=
L
(
u
)
.
According to the HPM, we can first use the embedding parameter p as a small parameter and assume that the solution of (2) can be written as a power series in p:
(8)
u
=
u
0
+
p
u
1
+
p
2
u
2
+
⋯
,
and the exact solution is obtained as follows:
(9)
x
=
lim
⁡
p
→
1
u
=
lim
⁡
p
→
1
(
u
0
+
p
u
1
+
p
2
u
2
+
p
3
u
3
+
⋯
)
=
∑
j
=
0
∞
u
j
.
Putting (8) into (5) and comparing the coefficients of identical degrees of p on both sides, we find
(10)
p
0
:
u
0
=
w
0
p
1
:
A
-
I
u
0
+
u
1
-
w
0
-
b
=
0
,
u
1
=
b
-
A
-
I
u
0
+
w
0
,
p
2
:
A
-
I
u
1
+
u
2
=
0
,
u
2
=
-
A
-
I
u
1
,
⋮
.
And in general
(11)
u
n
+
1
=
-
(
A
-
I
)
u
n
,
n
=
1,2
,
…
.
Taking u
0
=
w
0
=
0 yields
(12)
u
1
=
b
,
u
2
=
-
(
A
-
I
)
u
1
=
-
(
A
-
I
)
b
,
u
3
=
-
(
A
-
I
)
u
2
=
A
-
I
2
b
,
⋮
u
n
+
1
=
-
(
A
-
I
)
u
n
=
-
1
n
A
-
I
n
b
.
Therefore, the solution can be of the form
(13)
x
=
∑
i
=
0
∞
-
1
i
A
-
I
i
b
.
This method is efficient for solving system of linear equations but, however, does not converge for some systems when the spectral radius is greater than one. To make the method useful, Liu in [13] added the auxiliary parameter and the auxiliary matrix to the HPM. This method is as follows:
(14)
H
(
u
,
p
)
=
(
1
-
p
)
F
(
u
)
+
(
q
H
)
p
L
(
u
)
=
0
,
where
(15)
L
(
u
)
=
A
u
-
b
,
F
(
u
)
=
Q
u
-
w
0
.
Here, q is an auxiliary parameter, H is the auxiliary matrix, and the operator F
(
u
) is decided by the splitting matrix Q. Then his iterative scheme for u is
(16)
u
n
+
1
=
I
-
q
Q
-
1
H
A
n
(
q
Q
-
1
H
)
b
.
Therefore,
(17)
u
=
∑
i
=
0
∞
u
i
p
i
=
∑
i
=
0
∞
[
I
-
q
Q
-
1
H
A
i
(
q
Q
-
1
H
)
b
]
p
i
.
And finally by setting p
=
1, he obtained the solution as follows:
(18)
u
=
∑
i
=
0
∞
u
i
=
∑
i
=
0
∞
I
-
q
Q
-
1
H
A
i
(
q
Q
-
1
H
b
)
.
This author has adapted the Richardson, Jacobi, and the Gauss-Seidel methods to choose the splitting matrix and obtained that the homotopy series converged rapidly for a large sparse system with a small spectral radius. But to improve the rate of convergence, he does not present on a method for choosing H (see [13], conclusion). The main goal on this paper is studying this problem and doing some modifications for having very effective and straightforward results.Our idea for choosingH is preconditioning technique. Preconditioning methods are the most authoritative techniques to improve poor properties of the basic iterative methods. The main aim of preconditioning method for A
X
=
b is to substitute the original matrix A with an equivalent one A
prec which has better properties concerning the computation of a solution (generally by a certain iterative method). The two matrices are equivalent in the sense that they have the same solution. Simple preconditioners of this type are the left matrix preconditioners. The left preconditioned matrix is a nonsingular matrix M and the preconditioned matrix is defined by A
prec
=
M
A, where M
≈
A
-
1.In this paper, we considerH as a preconditioner and if the auxiliary matrix H is properly chosen, we can be able to achieve the best results. To improve the convergence rate of the basic iterative method, various models of preconditioning systems have been proposed [17–19]. In the literature, various authors have suggested different models of (
I
+
S
)-Type preconditioned [20–28], where all of these models are left preconditioners. We can consider all of the models in [17–28] for H.Here we consider Usui et al. [22] model of (
I
+
S)-Type preconditioners for H, where
(19)
H
=
h
i
j
n
×
n
=
1
,
j
=
i
,
-
a
i
j
a
j
j
,
j
>
i
,
0
,
otherwise
.
Let A
=
D
-
L
-
U, where D is the diagonal matrix and L, U are strictly lower and strictly upper triangular matrices of A, respectively. Then the preconditioned matrix is as follows:
(20)
A
prec
=
H
L
u
=
H
A
u
-
b
prec
,
b
prec
=
H
b
.
We know that H
=
(
I
+
S
), where
(21)
S
=
(
s
i
j
)
=
(
-
a
i
j
a
j
j
)
,
∀
j
>
i
.
Then we have
(22)
H
A
=
(
I
+
S
)
(
D
-
L
-
U
)
=
(
D
-
L
-
U
)
+
S
D
-
S
L
-
S
U
.
Thus,
(23)
H
A
=
D
-
D
1
-
L
+
L
1
-
U
-
S
D
+
U
1
+
S
U
,
where D
1, L
1, and U
1 are, respectively, the diagonal, strictly lower, and strictly upper triangular parts of S
L
=
D
1
+
L
1
+
U
1.Therefore, we have(24)
D
-
=
D
-
D
1
,
L
-
=
L
+
L
1
,
U
-
=
U
-
S
D
+
U
1
+
S
U
.
Furthermore, for more convergence speed in Jacobi-HPM and Gauss-Seidel-HPM we choose Q as follows.In Jacobi-HPM,(25)
Q
=
D
-
.
In Gauss-Seidel-HPM,
(26)
Q
=
D
-
-
L
-
.
Now, we will show that by this choice of H, Q these methods are faster than of the basic methods form point of view of the convergence speed.Definition 1 (see [29, 30]).
(a) A matrixA
=
a
i
j is called a Z-matrix if for any i
≠
j, a
i
j
≤
0.
(b) AZ-matrix is an M-matrix, if A is nonsingular and A
-
1
≥
0.Definition 2 (see [29, 30]).
LetA be a real matrix. The splitting A
=
M
-
N is called(a)
convergent ifρ
(
M
-
1
N
)
<
1 (ρ
(
·
) is the spectral radius of matrix);
(b)
regular ifM
-
1
≥
0 and N
≥
0.Lemma 3 (see [29]).
Let A be aZ-matrix. Then A is M-matrix if and only if there is a positive vector X such that A
X
>
0.Lemma 4 (see [31]).
LetA, B be Z-matrix and A an M-matrix; if A
≤
B, then B is an M-matrix too.Lemma 5 (see [32]).
LetA
=
M
-
N be a regular splitting. Then ρ
(
M
-
1
N
)
<
1 if and only if A is nonsingular and A
-
1 is nonnegative.Lemma 6 (see [24], Theorem 2.7, Remark 2.20).
LetA
=
[
D
-
(
L
+
U
)
]
=
[
(
D
-
L
)
-
U
] and H
A
=
[
D
-
-
(
L
-
+
U
-
)
=
(
D
-
-
L
-
)
-
U
-
] be regular splittings of A and H
A where H is (19). If A is M-matrix, then
(27)
ρ
(
D
-
-
1
(
L
-
+
U
-
)
)
≤
ρ
(
D
-
1
(
L
+
U
)
)
<
1
,
(28)
ρ
D
-
-
L
-
-
1
U
-
≤
ρ
D
-
L
-
1
U
<
1
.
By applying Lemma 3, we can achieve the following results.Theorem 7.
LetA be M-matrix. If we use H in (19) and Q in (25) and (26), then convergence rate of (18) is faster than of the convergence rate in ∑
i
=
0
∞
I
-
q
Q
^
-
1
A
i
(
q
Q
-
1
b
), (
Q
^
=
D or Q
^
=
D
-
L
).Proof.
SinceA is an M-matrix, by Lemma 3 it is easy to see that H
A also is M-matrix, and thus by Lemma 4
Q is M-matrix and Q
-
q
H
A is nonnegative. Then by Lemmas 5 and 6 we can obtain ρ
(
I
-
q
Q
-
1
H
A
)
≤
ρ
(
I
-
q
Q
^
-
1
A
)
<
1.
Therefore by [13, Theorem 2.1] the proof is completed.


## 3. Numerical Experiments

In this section, we give an example to illustrate the obtained results in the previous sections.We test convection-diffusion equation (1) when δ
=
1, τ
=
2.Then, we solved then
2
×
n
2
M-matrix yielded by Jacobi-HPM and Gauss-Seidel-HPM.In this experiment, we choose the right hand side vector, such thatX
=
1,1
,
…
,
1
T is solution of A
X
=
b. The stopping criterion with tolerance ε
=
1
0
-
10 is A
u
-
b
<
ε.In Tables1 and 2 we report theCPU time and number of iterations (Iter) for the Jacobi-HPM and Gauss-Seidel-HPM methods, respectively.Table 1
The results of experiment for Jacobi-HPM.
Method
Algorithm forq
=
1, H
=
I
Algorithm forq
=
1, H
=
(
I
+
S
) in (19)
n
N
Iter
CPU
Iter
CPU
5
25
42
0.015893
24
0.003804
10
100
63
0.501852
37
0.231103
15
225
81
8.018167
47
3.844972
20
400
96
71.971452
56
35.328229Table 2
The results of experiment for Gauss-Seidel-HPM.
Method
Algorithm forq
=
1, H
=
I
Algorithm forq
=
1, H
=
(
I
+
S
) in (19)
n
N
Iter
CPU
Iter
CPU
5
25
24
0.007968
10
0.001203
10
100
38
0.255449
17
0.069256
15
225
51
4.428433
22
1.424873
20
400
62
40.541738
27
13.336331From the tables, we can see that ifH is (19), then the mentioned methods are superior to the basic methods.


## 4. Conclusions

In this paper, we have applied a new method for solving convection-diffusion equations based on the HPM. We can see that how HPM method is influenced for solving system of linear equation, if combined by preconditioning technique. Finally, from theoretical speaking and numerical experiment, it may be concluded that the convergence behaviors of this method is superior to some other methods in the frame of HPM.



---

*Source: 143512-2014-11-16.xml*