- Email: [email protected]
A fast algorithm for solving Toeplitz penta-diagonal systems
Applied Mathematics and Computation 215 (2010) 3830–3838
Contents lists available at ScienceDirect
Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc
A fast algorithm for solving Toeplitz penta-diagonal systems q S.S. Nemani Department of Mathematics and Computer Science, Georgian Court University, Lakewood, NJ 08701, United States
a r t i c l e
i n f o
a b s t r a c t In this paper a fast algorithm for solving a large system with a symmetric Toeplitz pentadiagonal coefficient matrix is presented. This efficient method is based on the idea of a system perturbation followed by corrections and is competitive with standard methods. The error analysis is also given. Ó 2009 Elsevier Inc. All rights reserved.
Keywords: Factorization Toeplitz matrix Penta-diagonal matrix Perturbed system
1. Introduction Consider an n n system of linear equations given by
Ax ¼ b; where
2
d
a
1
3
7 6 a 1 7 6a d 7 6 7 61 a d a 1 7 6 7 6 . . . . . .. .. .. .. .. 7 6 7 6 A¼6 7 . . . . . 7 6 . . . . . . . . . . 7 6 7 6 7 6 1 a d a 1 7 6 7 6 1 a d a5 4 1 a d and jdj > 2jaj þ 2. This means the matrix A is strictly diagonally dominant. This special kind of system appears in many applications: spline approximation [1], solution of certain differential equations [3,5,7], prediction [8], etc. The authors in [4,6,9] consider the solution of symmetric circulant tridiagonal and strictly diagonally dominant matrices. They perturb the matrix so that it consists of a near-Toeplitz matrix and a correction matrix. In this paper, a fast algorithm for solving a large system with a symmetric Toeplitz penta-diagonal coefficient matrix is presented. The operation counts for the method considered take into account the basic arithmetic operations and weights them equally. The error analysis is given using the max norm. Some numerical examples are presented to illustrate the efficiency and stability of the method described in this paper. 2. Computation of perturbation matrix Consider an n n penta-diagonal Toeplitz system of linear equations, given by q
Supported by GCU Faculty Summer Research grant. E-mail address: [email protected]
0096-3003/$ - see front matter Ó 2009 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2009.11.026
S.S. Nemani / Applied Mathematics and Computation 215 (2010) 3830–3838
Ax ¼ b;
3831
ð2:1Þ
where
3 d a 1 7 6 a 1 7 6a d 7 6 7 61 a d a 1 7 6 7 6 . . . . . .. .. .. .. .. 7 6 7 6 A¼6 7: . . . . . 7 6 . . . . . . . . . . 7 6 7 6 7 6 1 a d a 1 7 6 7 6 1 a d a5 4 1 a d 2
Conditions on the scalars a and d are given in Theorem 2.1. In order to solve Eq. (2.1), we need to perturb the system as follows:
A ¼ A0 þ C so that A0 has the Toeplitz LU factorization and the matrix C has very few non-zero elements. Let
2
3
1
6 6 l1 6 6 L ¼ 6 l2 6 6 4
7 7 7 7 7 7 7 5
1 l1 1 .. .. .. . . . l2
l1
ð2:2Þ
1
and
2 6 6 6 6 U¼6 6 6 4
u1
u2 1 .. .. .. . . . u1
u2 u1
3 7 7 7 7 ; 17 7 7 u2 5 u1
ð2:3Þ
where l1 ; l2 ; u1 and u2 form the solution of the following non-linear equations:
l2 u1 ¼ 1;
ð2:4Þ
l2 u2 þ l1 u1 ¼ a;
ð2:5Þ
l2 þ l1 u2 þ u1 ¼ d;
ð2:6Þ
l1 þ u2 ¼ a:
ð2:7Þ 0
Then it is easy to show that LU ¼ A where
2
u1 6 6 l1 u1 6 6 1 6 6 6 6 0 A ¼6 6 6 6 6 6 6 4
u2 l1 u2 þ u1 a ..
.
3
1
7 7 7 7 d a 1 7 7 .. .. .. .. 7 . . . . 7 7: .. .. .. .. .. 7 . . . . . 7 7 1 a d a 17 7 7 1 a d a5 1 a d a
1
ð2:8Þ
Thus,
C ¼ ½ðl2 þ l1 u2 Þe1 þ ðl2 u2 Þe2 eT1 þ ½l1 e1 þ l2 e2 eT2 ; where ei ¼ ð0; . . . ; 1; 0; . . . ; 0ÞT , 1 at the ith co-ordinate. The superscript T corresponds to the transpose operation. Note that in order to know the matrix A0 , we need to know l1 ; l2 ; u1 and u2 . Instead of solving the non-linear equations 2 (2.4)–(2.7) directly, we introduced an auxiliary polynomial defined as f ðxÞ ¼ x4 þ ax3 þ dx þ ax þ 1.
3832
S.S. Nemani / Applied Mathematics and Computation 215 (2010) 3830–3838
Theorem 2.1. Let the scalar a and d be such that
jdj > 2jaj þ 2
ð2:9Þ
and for d > 0,
d<
a2 þ 8 : 4
ð2:10Þ 2
If r1 ; r 2 ; r3 and r 4 are the roots of f ðxÞ ¼ x4 þ ax3 þ dx þ ax þ 1 with jr 1 j 6 jr 2 j 6 jr3 j 6 jr 4 j, then jr 1 j < jr 2 j < 1 and 1 < jr 3 j < jr 4 j. Proof. Note that f ð0Þ ¼ 1 means 0 is not a solution to the equation f ðxÞ ¼ 0. Divide the equation by x2 to obtain
x2 þ
1 1 þ a x þ þ d ¼ 0: x2 x
After completing the square and setting y ¼ x þ 1x, the above equation reduces to
y2 þ ay þ ðd 2Þ ¼ 0:
ð2:11Þ
Let y1 ; y2 be solutions of Eq. (2.11). Then
y1 ¼
a þ
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi a2 4d þ 8 2
a
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi a2 4d þ 8 : 2
and
y2 ¼
The condition (2.10) guarantees that y1 and y2 are both real and distinct. Next, we show that jyi j > 2 for i ¼ 1; 2. Case (i) d < 0. pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi From Eq. (2.9), we get a2 4d þ 8 > jaj2 þ 8jaj þ 16 ¼ ðjaj þ 4Þ2 . Therefore, a2 4d þ 8 > jaj þ 4. (This follows from the fact that for any two positive real numbers x and y with x2 > y2 implies x > y.) Subcase (i) a > 0. Since a þ jaj ¼ 0, we get
y1 ¼
a þ
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi a2 4d þ 8 a þ ðjaj þ 4Þ ¼ 2: > 2 2
Since a jaj ¼ 2a < 0, we get
y2 ¼
a
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi a2 4d þ 8 a ðjaj þ 4Þ < 2: < 2 2
Subcase (ii) a < 0. Since a þ jaj ¼ 2a > 0, we get
y1 ¼
a þ
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi a2 4d þ 8 a þ ðjaj þ 4Þ > >2 2 2
and since a jaj ¼ 0, we get
y2 ¼
a
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi a2 4d þ 8 a ðjaj þ 4Þ < ¼ 2: 2 2
Case (ii) d > 0. 2 Using conditions (2.9) and (2.10), we get 2jaj þ 2 < d < a 4þ8. Therefore, 8 < jaj. Thus a2 4d þ 8 < jaj2 8jaj pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi < ðjaj 4Þ2 . This means that a2 4d þ 8 < jðjaj 4Þj ¼ jaj 4. Subcase (i) a > 0.
y1 ¼
a þ
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi a2 4d þ 8 a þ ðjaj 4Þ < ¼ 2: 2 2
Since y2 < y1 , we get y2 < 2. Subcase (ii) a < 0.
y2 ¼
a
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi a2 4d þ 8 a ðjaj 4Þ ¼ 2: > 2 2
Since y2 < y1 , we get y1 > 2.
S.S. Nemani / Applied Mathematics and Computation 215 (2010) 3830–3838
3833
Thus combine all cases, we have jyi j > 2 for i ¼ 1; 2. Since x2 yi x þ 1 ¼ 0 for i ¼ 1; 2 applying quadratic formula gives us
xi; ¼
yi
qffiffiffiffiffiffiffiffiffiffiffiffiffiffi y2i 4 2
:
Since jyi j > 2 for i ¼ 1; 2, the roots xi; are real and distinct. For yi > 2; xi;þ > 1 and jxi; j < 1. For yi < 2; xi; < 1 and jxi;þ j < 1. h From Theorem 2.1 we know that xi; xi;þ ¼ 1. This gives us that r 4 ¼ r 1 and r 3 ¼ r 1 1 2 . Let us assume that r; s 2 fri : 1 6 i 6 4g be such that jr þ sj þ jrsj < 1. (Existence of such r and s for some cases and more information about the roots are described in Appendix.) This means that rs – 1. Therefore, fr i : 1 6 i 6 4g ¼ fr; s; r 1 ; s1 g. Hence we have
f ðxÞ ¼ ðx rÞðx sÞðx s1 Þðx r 1 Þ; which is equal to
x4 ðr þ s þ s1 þ r 1 Þx3 þ ðrs þ rs1 þ 2 þ sr1 þ s1 r 1 Þx2 ðr þ s þ s1 þ r 1 Þx þ 1: Hence comparing the coefficients we get
r þ s þ s1 þ r 1 ¼ a
ð2:12Þ
rs þ rs1 þ 2 þ sr 1 þ s1 r1 ¼ d:
ð2:13Þ
and
Thus, if we set
l1 ¼ ðr þ sÞ;
l2 ¼ rs;
u2 ¼ ðr 1 þ s1 Þ and u1 ¼ r1 s1 :
Then using Eqs. (2.12) and (2.13), it can be shown that l1 ; l2 ; u1 and u2 satisfy Eqs. (2.4)–(2.7). Note that jl1 j þ jl2 j ¼ jr þ sj þ jrsj < 1. Dividing by jrsj both sides in this inequality, we get ju1 j > ju2 j þ 1 and hence both matrices L and U are strictly diagonally dominant. The computation of li and ui for i ¼ 1; 2 gives rise to a Toeplitz LU factorization of the matrix A0 in Oð1Þ time. To solve Eq. (2.1), we first solve A0 x0 ¼ b. To find x0 , we first solve Ly ¼ b and then solve Ux0 ¼ y by forward and backward substitution, respectively. To obtain the vector x0 , it takes 8n operations. The diagonally dominant property of matrices L and U guarantees the numerical stability of the process of solving equation A0 x0 ¼ b. To estimate kx0 k we need the followings. Lemma 2.1. Let L be the matrix as given in Eq. (2.2) and jl1 j þ jl2 j < 1. If Ly ¼ b, then
kyk 6
kbk : ð1 jl1 j jl2 jÞ
Proof Case (i) kyk ¼ jyi j for some i; 3 6 i 6 n. Since Ly ¼ b we get l2 yi2 þ l1 yi1 þ yi ¼ bi . Therefore,
kyk ¼ jyi j 6 jbi j þ jl2 jjyi2 j þ jl1 jjyi1 j 6 kbk þ jl2 jkyk þ jl1 jkyk: Thus
kyk 6
kbk : ð1 jl1 j jl2 jÞ
Case (ii) kyk ¼ jy1 j. Therefore,
kyk ¼ jy1 j ¼ jb1 j 6 kbk: Case (iii) kyk ¼ jy2 j, we have
kyk ¼ jy2 j ¼ jb2 jl1 y1 j 6 kbk þ jl1 jkyk: Thus
kyk 6
kbk : ð1 jl1 jÞ
The result follows from the fact that 1 jl1 j jl2 j 6 1 jl1 j 6 1.
h
Similarly, using the fact that ju1 j ju2 j 1 6 ju1 j ju2 j 6 ju1 j, we get the following lemma.
3834
S.S. Nemani / Applied Mathematics and Computation 215 (2010) 3830–3838
Lemma 2.2. Let U be the matrix as given in Eq. (2.3) and ju1 j > ju2 j þ 1. If Ux0 ¼ y, then
kx0 k 6
kyk : ðju1 j ju2 j 1Þ
Now we are ready to calculate a bound for kx0 k. It is given in the following theorem. Theorem 2.2. Let A0 be the matrix as given in Eq. (2.8) . If A0 x0 ¼ b, then
kx0 k 6
kbk : ð1 jl1 j jl2 jÞðju1 j ju2 j 1Þ
Proof. Since x0 is solution of A0 x0 ¼ b; x0 is solution of Ux0 ¼ y, where y is solution of Ly ¼ b. Hence from Lemmas 2.2 and 2.1 we get
kx0 k 6
kyk kbk 6 : ðju1 j ju2 j 1Þ ð1 jl1 j jl2 jÞðju1 j ju2 j 1Þ
3. Correction procedure Note that A ¼ A0 þ C implies that
Ax0 ¼ b þ Cx0 ; 0
where Cx ¼ ½ðl2 þ 0
l1 u2 Þx01
1
ð3:1Þ þ
l1 x02 e1
þ
½ðl2 u2 Þx01
þ
l2 x02 e2 .
We get the solution to Ax ¼ b as follows:
1
x ¼ x aA e1 bA e2 ; where a ¼ ðl2 þ l1 u2 Þx01 þ l1 x02 and b ¼ ðl2 u2 Þx01 þ l2 x02 . If we can find vectors wðiÞ such that AwðiÞ ¼ ei for i ¼ 1; 2 then we can solve Ax ¼ b approximately by adding a linear combination of vectors wð1Þ and wð2Þ to the vector x0 . P Let uðiÞ ¼ tj¼1 r ji ej for i ¼ 1; 2 where t is to be chosen based on a predefined tolerance g. Using the fact that r1 and r 2 are 2 roots of the polynomial x4 þ ax3 þ dx þ ax þ 1, we get for i ¼ 1; 2 tþ1 AuðiÞ ¼ ðr1 et1 ðar tþ1 þ rtþ2 Þet þ ðr t1 þ ar ti Þetþ1 þ r ti etþ2 : i þ aÞe1 e2 r i i i i r2 r2 ð2Þ ð1Þ Let wð1Þ ¼ rr11r ðuð1Þ uð2Þ Þ and wð2Þ ¼ rr11r ððr1 ðr1 1 þ aÞu 2 þ aÞu Þ. 2 2 Then
Awð1Þ ¼ e1 þ
r1 r2 X p ej r 1 r 2 j2C j
ð3:2Þ
Awð2Þ ¼ e2 þ
r1 r2 X q ej ; r 1 r 2 j2C j
ð3:3Þ
and
where C ¼ ft 1; t; t þ 1; t þ 2g,
pt1 ¼ r tþ1 rtþ1 2 1 ; pt ¼ ða þ r 2 Þr tþ1 ða þ r 1 Þr 1tþ1 ; 2 t 1 t ptþ1 ¼ ðr 1 1 þ aÞr 1 ðr 2 þ aÞr 2 ;
ptþ2 ¼ r t1 r t2 and tþ1 tþ1 qt1 ¼ ðr 1 ðr1 2 þ aÞr 1 1 þ aÞr 2 ; tþ1 tþ1 qt ¼ ðr 1 ðr1 2 þ aÞða þ r 1 Þr 1 1 þ aÞða þ r 2 Þr 2 ; 1 t 1 1 t qtþ1 ¼ ðr 1 1 þ aÞðr 2 þ aÞr 2 ðr 2 þ aÞðr 1 þ aÞr 1 ; t 1 t qtþ2 ¼ ðr 1 1 þ aÞr 2 ðr 2 þ aÞr 1 :
Let
~ ¼ x0 awð1Þ bwð2Þ : x
ð3:4Þ
S.S. Nemani / Applied Mathematics and Computation 215 (2010) 3830–3838
3835
~ ¼ Ax0 aAwð1Þ bAwð2Þ . From Eqs. (3.1)–(3.3) we get Then Ax
~b¼ Ax
r1 r2 X ðapj þ bqj Þej : r 1 r 2 j2C
After rearranging the terms, we get
~b¼ Ax
r 1 r2 X ð1Þ t ð2Þ ðd r þ dj r t2 Þej ; r 1 r 2 j2C j 1
ð3:5Þ
where ð1Þ
dt1 ¼ fbðr 1 2 þ aÞ agr 1 ; ð1Þ
dt ¼ fbðr 1 2 þ aÞ agða þ r 1 Þr 1 ; ð1Þ
1 dtþ1 ¼ fa bðr 1 2 þ aÞgðr 1 þ aÞ; ð1Þ
dtþ2 ¼ a bðr 1 2 þ aÞ and ð2Þ
dt1 ¼ fa bðr 1 1 þ aÞgr 2 ; ð2Þ
dt ¼ fa bðr 1 1 þ aÞgða þ r 2 Þr 2 ; ð2Þ
1 dtþ1 ¼ fbðr 1 1 þ aÞ agðr 2 þ aÞ; ð2Þ
dtþ2 ¼ bðr1 1 þ aÞ a: 1 From Theorem 2.1 we know that r 1 1 ¼ r 4 ; r 2 ¼ r 3 and jr 1 j < jr 2 j < 1. Therefore, for all j 2 C,
ð1Þ
jdj j 6 ðjaj þ jbjjr 3 þ ajÞ maxf1; ja þ r 1 j; jr4 þ ajg and ð2Þ
jdj j 6 ðjaj þ jbjjr 4 þ ajÞ maxf1; ja þ r 2 j; jr3 þ ajg: Recall that a ¼ ðl2 þ l1 u2 Þx01 þ l1 x02 and b ¼ ðl2 u2 Þx01 þ l2 x02 . Using triangle inequalty and the definition of the max norm we obtain
jaj 6 ðjl1 j þ jl2 j þ jl1 jju2 jÞkx0 k and
jbj 6 jl2 jðju2 j þ 1Þkx0 k: Then using Eq. (3.5) we get
~ bk 6 kAx
jr 1 r 2 j fM1 þ M 2 gjr 2 jt kx0 k; jr1 r2 j
where M1 ¼ fðjl1 j þ jl2 j þ jl1 jju2 jÞ þ jl2 jðju2 j þ 1Þjr3 þ ajg maxf1; ja þ r 1 j; jr 4 þ ajg and M 2 ¼ fðjl1 j þ jl2 j þ jl1 jju2 jÞ þ jl2 jðju2 j þ 1Þ jr 4 þ ajg maxf1; ja þ r2 j; jr3 þ ajg. Using the bound for kx0 k from Theorem 2.2, we have the following theorem. ~ be as given in Eq. (3.4) . Then Theorem 3.1. Let x
~ bk 6 Mjr 2 jtþ2 kbk; kAx where
M¼
M1 þ M2 : jr 1 r 2 jð1 jl1 j jl2 jÞðju1 j ju2 j 1Þ
Note that jr2 j < 1 so the error decreases exponentially with increasing value of t. The number of components t that needs to be corrected is determined in relation to a specific tolerance g for the desired accuracy. For example, if we required ~ bk 6 gkbk, then from Theorem 3.1, t must be chosen to satisfy kAx
tðgÞ P
log g log M 2: log jr2 j
Substituting the vectors wð1Þ and wð2Þ in Eq. (3.4), we have
~ ¼ x0 x
r1 r2 r1 r2 ða bðr3 þ aÞÞuð1Þ þ ða bðr4 þ aÞÞuð2Þ : r1 r2 r1 r2
ð3:6Þ
3836
S.S. Nemani / Applied Mathematics and Computation 215 (2010) 3830–3838
Table 1 Algorithm for the method. Input
a and d so that A is as given in Eq. (2.1) The size of matrix A ¼ n, the desired tolerance g The right hand side vector b (1) r1 ; r2 ; r3 and r 4 as given in Theorem 2.1 (2) Solve Ly ¼ b, using forward substitution (3) Solve Ux0 ¼ y, using backward substitution l m glog M t ¼ loglog jr2 j
Calculate
(4) Do correction as given in Eq. (3.6) ~ to Ax ¼ b An approximate solution x
Output
Since jr1 j and jr2 j are less than 1, computing the value of ~ xi for increasing value of i will ensure the numerical stability of the ~ from the vector x0 , we need only 4t operations. Hence the total computational efmethod. Note that to obtain the vector x forts for the method given in this paper is 8n þ 4t þ Oð1Þ operations. For sufficiently large value of n, we get t n, implying this method takes 8n þ Oð1Þ operations. The algorithm for this method is given in Table 1. 4. Numerical result Table 2 provides feedback on the number of terms of t to be corrected for a given desired accuracy in the result. The values ~ bk, error, of t were computed with d ¼ ð2jaj þ 2 þ DÞ for different values of a and D. Table 3 presents the values of kAx with the right hand side vector b ¼ ð1; 1; . . . ; 1ÞT and n ¼ 2000. Note that the tolerance g cannot be achieved when tðgÞ > n, and hence the algorithm will break down when the required relative tolerance is sufficiently small and/or jdj is sufficiently close to 2jaj þ 2 (that is, D is sufficiently close to 0). Table 2 indicate that the value of t for large values of D is less ~ bk for large values of D is than the value of t for small values of D. Also note that from Table 3 it is clear that the value of kAx less than the error for small values of D. 5. Remark The method described in this paper is very efficient and stable one. It compares favorably with the regular LU method. Also this method is competitive with other methods [2] in term of arithmetic operations. Following the parallel approach described in [6], parallel implementation can be done in a similar way. Table 2 t-Values for the method.
D
g 1.0000E01
1.0000E03
1.0000E05
1.0000E07
For a ¼ 1 0.001 0.01 0.05 0.1 0.5 2 4 6
831 211 78 50 17 6 4 3
1156 314 124 83 32 14 9 7
1482 417 170 115 47 22 15 12
1808 520 216 148 61 29 21 17
For a ¼ 0:5 0.001 0.01 0.05 0.1 0.5 2 4 6
747 187 68 43 14 5 3 2
1056 285 112 74 28 12 8 7
1365 382 155 105 42 20 14 11
1674 480 199 136 56 27 19 16
For a ¼ 0:1 0.001 0.01 0.05 0.1 0.5 2 4 6
664 163 58 36 12 4 2 1
959 256 100 66 25 11 7 6
1254 349 141 95 38 18 13 10
1549 443 183 125 52 25 18 15
S.S. Nemani / Applied Mathematics and Computation 215 (2010) 3830–3838
3837
Table 3 ~ bk. The values of kAx
g
D
1.0000E01
1.0000E03
1.0000E05
1.0000E07
For a ¼ 1 0.001 0.01 0.05 0.1 0.5 2 4 6
7.8134E03 7.7988E03 7.8132E03 7.9590E03 8.1115E03 8.5965E03 4.7283E03 4.7218E03
7.8849E05 7.8002E05 7.9049E05 7.5814E05 7.5373E05 7.0900E05 8.6140E05 9.2079E05
7.8453E07 7.8016E07 7.9976E07 8.3154E07 7.0038E07 5.7866E07 6.4396E07 7.2304E07
7.8052E09 7.8028E09 8.0914E09 7.9209E09 8.8902E09 8.6133E09 4.8447E09 5.9141E09
For a ¼ 0:5 0.001 0.01 0.05 0.1 0.5 2 4 6
1.0859E02 1.0896E02 1.1051E02 1.1589E02 1.2336E02 1.0927E02 9.0311E03 5.2017E03
1.0849E04 1.0755E04 1.0778E04 1.1584E04 1.2520E04 1.2457E04 1.0229E04 5.7487E05
1.0838E06 1.1128E06 1.1678E06 1.1578E06 1.2701E06 8.2621E07 6.6942E07 1.0113E06
1.0827E08 1.0984E08 1.1390E08 1.1572E08 1.2885E08 1.0194E08 9.9379E09 6.4563E09
For a ¼ 0:1 0.001 0.01 0.05 0.1 0.5 2 4 6
1.7136E02 1.7205E02 1.7422E02 1.8725E02 1.3863E02 9.3424E03 7.3716E03 2.0686E02
1.7106E04 1.7449E04 1.7013E04 1.7609E04 1.7138E04 1.5216E04 1.9746E04 5.2503E05
1.7076E06 1.7696E06 1.8549E06 1.9345E06 1.9696E06 1.2285E06 9.8031E07 9.5138E07
1.7045E08 1.7082E08 1.8114E08 1.8190E08 1.6297E08 1.5108E08 8.3295E09 1.0580E08
Appendix Case (i) d > 0 and a > 0. From Theorem 2.1, we have y2 < y1 < 2. Therefore,
qffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffi y21 4 < y22 4:
Hence x2; < x1; < 0 means jx1; j < jx2; j. Since xi;þ xi; ¼ 1 for i ¼ 1; 2, we get that the signs of xi;þ and xi; are same. Thus, all the roots are negative. Therefore, in this case we will have r 1 ¼ x1; ; r2 ¼ x2; ; r3 ¼ x2;þ and r 4 ¼ x1;þ . Case (ii) d > 0 and a < 0. In this case, 2 < y2 < y1 , therefore, x2;þ < x1;þ . The sign of all the roots in this case is positive. In this case we have r1 ¼ x1; ; r 2 ¼ x2; ; r 3 ¼ x2;þ and r 4 ¼ x1;þ . Case (iii) d < 0. In this case, y1 > 2 and y2 < 2. Thus x1; > 0 and x2; < 0. Subcase (i) jx1; j < jx2;þ j. In this case, we have r1 ¼ x1; ; r 2 ¼ x2;þ ; r 3 ¼ x2; and r4 ¼ x1;þ . Then r1 < r 2 . Therefore, ðr 1 þ r 2 Þ > 0 and r 1 r2 < 0. Consider 1 jr 1 þ r 2 j jr1 r 2 j ¼ 1 þ r1 þ r2 þ r 1 r 2 ¼ ð1 þ r1 Þð1 þ r 2 Þ. Since r 1 > 0; r1 þ 1 > 0 and since 1 < r 2 < 0 we get 1 þ r 2 > 0. Thus, 1 jr1 þ r2 j jr 1 r2 j > 0. Subcase (ii) jx1; j > jx2;þ j. In this case, we have r 1 ¼ x2;þ ; r 2 ¼ x1; ; r3 ¼ x1;þ and r 4 ¼ x2; . Then r 1 < r 2 . Therefore, ðr1 þ r 2 Þ > 0 and r 1 r2 < 0. Consider 1 jr1 þ r 2 j jr 1 r 2 j ¼ 1 r 1 r2 þ r1 r2 ¼ ð1 r1 Þð1 r2 Þ. Since r1 > 0; 1 r 1 > 0 and since 0 < r2 < 1 we get 1 r2 > 0. Thus, 1 jr 1 þ r 2 j jr1 r 2 j > 0.Combine both subcases, we get jr 1 þ r 2 j þ jr 1 r 2 j < 1. Thus take r ¼ r 1 and s ¼ r 2 . References [1] J.H. Ahlberg, E.N. Nilson, J.L. Walsh, The Theory of Splines and their Applications, Academic Press, New York, 1967. [2] F. Diele, L. Lopez, The use of the factorization of five-diagonal matrices by tridiagonal Toeplitz matrices, Appl. Math. Lett. 11 (1998) 51–59. [3] P.G. Patil, Y.S. Swamy, An efficient model for vibration control by piezoelectric smart structure using finite element method, Eur. J. Comput. Sci. Netw. Secu. 8 (11) (2008) 258–264. [4] O. Rojo, A new method for solving symmetric circulant tri-diagonal system of linear equations, Comput. Math. Appl. 20 (1990) 61–67.
3838
S.S. Nemani / Applied Mathematics and Computation 215 (2010) 3830–3838
[5] S.S. Nemani, L.E. Garey, An efficient method for solving second order boundary value problems with two point boundary conditions, Int. J. Comput. Math. 79 (9) (2002) 1001–1008. [6] S.S. Nemani, L.E. Garey, Parallel algorithms for solving tridiagonal and near-circulant systems, Appl. Math. Comput. 130 (2002) 285–294. [7] E.W. Sachs, A.K. Strauss, Efficient solution of a partial integro-differential equation in finance, Appl. Numer. Math. 58 (2008) 1687–1703. [8] P. Whittle, Prediction and Regulation by Linear Least Squares Methods, Van Nostrand, Princeton, NJ, 1963. [9] W.M. Yan, K.L. Chung, A fast algorithm for solving special tri-diagonal systems, Computing 52 (1994) 203–211.
Contents lists available at ScienceDirect
Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc
A fast algorithm for solving Toeplitz penta-diagonal systems q S.S. Nemani Department of Mathematics and Computer Science, Georgian Court University, Lakewood, NJ 08701, United States
a r t i c l e
i n f o
a b s t r a c t In this paper a fast algorithm for solving a large system with a symmetric Toeplitz pentadiagonal coefficient matrix is presented. This efficient method is based on the idea of a system perturbation followed by corrections and is competitive with standard methods. The error analysis is also given. Ó 2009 Elsevier Inc. All rights reserved.
Keywords: Factorization Toeplitz matrix Penta-diagonal matrix Perturbed system
1. Introduction Consider an n n system of linear equations given by
Ax ¼ b; where
2
d
a
1
3
7 6 a 1 7 6a d 7 6 7 61 a d a 1 7 6 7 6 . . . . . .. .. .. .. .. 7 6 7 6 A¼6 7 . . . . . 7 6 . . . . . . . . . . 7 6 7 6 7 6 1 a d a 1 7 6 7 6 1 a d a5 4 1 a d and jdj > 2jaj þ 2. This means the matrix A is strictly diagonally dominant. This special kind of system appears in many applications: spline approximation [1], solution of certain differential equations [3,5,7], prediction [8], etc. The authors in [4,6,9] consider the solution of symmetric circulant tridiagonal and strictly diagonally dominant matrices. They perturb the matrix so that it consists of a near-Toeplitz matrix and a correction matrix. In this paper, a fast algorithm for solving a large system with a symmetric Toeplitz penta-diagonal coefficient matrix is presented. The operation counts for the method considered take into account the basic arithmetic operations and weights them equally. The error analysis is given using the max norm. Some numerical examples are presented to illustrate the efficiency and stability of the method described in this paper. 2. Computation of perturbation matrix Consider an n n penta-diagonal Toeplitz system of linear equations, given by q
Supported by GCU Faculty Summer Research grant. E-mail address: [email protected]
0096-3003/$ - see front matter Ó 2009 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2009.11.026
S.S. Nemani / Applied Mathematics and Computation 215 (2010) 3830–3838
Ax ¼ b;
3831
ð2:1Þ
where
3 d a 1 7 6 a 1 7 6a d 7 6 7 61 a d a 1 7 6 7 6 . . . . . .. .. .. .. .. 7 6 7 6 A¼6 7: . . . . . 7 6 . . . . . . . . . . 7 6 7 6 7 6 1 a d a 1 7 6 7 6 1 a d a5 4 1 a d 2
Conditions on the scalars a and d are given in Theorem 2.1. In order to solve Eq. (2.1), we need to perturb the system as follows:
A ¼ A0 þ C so that A0 has the Toeplitz LU factorization and the matrix C has very few non-zero elements. Let
2
3
1
6 6 l1 6 6 L ¼ 6 l2 6 6 4
7 7 7 7 7 7 7 5
1 l1 1 .. .. .. . . . l2
l1
ð2:2Þ
1
and
2 6 6 6 6 U¼6 6 6 4
u1
u2 1 .. .. .. . . . u1
u2 u1
3 7 7 7 7 ; 17 7 7 u2 5 u1
ð2:3Þ
where l1 ; l2 ; u1 and u2 form the solution of the following non-linear equations:
l2 u1 ¼ 1;
ð2:4Þ
l2 u2 þ l1 u1 ¼ a;
ð2:5Þ
l2 þ l1 u2 þ u1 ¼ d;
ð2:6Þ
l1 þ u2 ¼ a:
ð2:7Þ 0
Then it is easy to show that LU ¼ A where
2
u1 6 6 l1 u1 6 6 1 6 6 6 6 0 A ¼6 6 6 6 6 6 6 4
u2 l1 u2 þ u1 a ..
.
3
1
7 7 7 7 d a 1 7 7 .. .. .. .. 7 . . . . 7 7: .. .. .. .. .. 7 . . . . . 7 7 1 a d a 17 7 7 1 a d a5 1 a d a
1
ð2:8Þ
Thus,
C ¼ ½ðl2 þ l1 u2 Þe1 þ ðl2 u2 Þe2 eT1 þ ½l1 e1 þ l2 e2 eT2 ; where ei ¼ ð0; . . . ; 1; 0; . . . ; 0ÞT , 1 at the ith co-ordinate. The superscript T corresponds to the transpose operation. Note that in order to know the matrix A0 , we need to know l1 ; l2 ; u1 and u2 . Instead of solving the non-linear equations 2 (2.4)–(2.7) directly, we introduced an auxiliary polynomial defined as f ðxÞ ¼ x4 þ ax3 þ dx þ ax þ 1.
3832
S.S. Nemani / Applied Mathematics and Computation 215 (2010) 3830–3838
Theorem 2.1. Let the scalar a and d be such that
jdj > 2jaj þ 2
ð2:9Þ
and for d > 0,
d<
a2 þ 8 : 4
ð2:10Þ 2
If r1 ; r 2 ; r3 and r 4 are the roots of f ðxÞ ¼ x4 þ ax3 þ dx þ ax þ 1 with jr 1 j 6 jr 2 j 6 jr3 j 6 jr 4 j, then jr 1 j < jr 2 j < 1 and 1 < jr 3 j < jr 4 j. Proof. Note that f ð0Þ ¼ 1 means 0 is not a solution to the equation f ðxÞ ¼ 0. Divide the equation by x2 to obtain
x2 þ
1 1 þ a x þ þ d ¼ 0: x2 x
After completing the square and setting y ¼ x þ 1x, the above equation reduces to
y2 þ ay þ ðd 2Þ ¼ 0:
ð2:11Þ
Let y1 ; y2 be solutions of Eq. (2.11). Then
y1 ¼
a þ
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi a2 4d þ 8 2
a
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi a2 4d þ 8 : 2
and
y2 ¼
The condition (2.10) guarantees that y1 and y2 are both real and distinct. Next, we show that jyi j > 2 for i ¼ 1; 2. Case (i) d < 0. pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi From Eq. (2.9), we get a2 4d þ 8 > jaj2 þ 8jaj þ 16 ¼ ðjaj þ 4Þ2 . Therefore, a2 4d þ 8 > jaj þ 4. (This follows from the fact that for any two positive real numbers x and y with x2 > y2 implies x > y.) Subcase (i) a > 0. Since a þ jaj ¼ 0, we get
y1 ¼
a þ
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi a2 4d þ 8 a þ ðjaj þ 4Þ ¼ 2: > 2 2
Since a jaj ¼ 2a < 0, we get
y2 ¼
a
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi a2 4d þ 8 a ðjaj þ 4Þ < 2: < 2 2
Subcase (ii) a < 0. Since a þ jaj ¼ 2a > 0, we get
y1 ¼
a þ
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi a2 4d þ 8 a þ ðjaj þ 4Þ > >2 2 2
and since a jaj ¼ 0, we get
y2 ¼
a
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi a2 4d þ 8 a ðjaj þ 4Þ < ¼ 2: 2 2
Case (ii) d > 0. 2 Using conditions (2.9) and (2.10), we get 2jaj þ 2 < d < a 4þ8. Therefore, 8 < jaj. Thus a2 4d þ 8 < jaj2 8jaj pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi < ðjaj 4Þ2 . This means that a2 4d þ 8 < jðjaj 4Þj ¼ jaj 4. Subcase (i) a > 0.
y1 ¼
a þ
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi a2 4d þ 8 a þ ðjaj 4Þ < ¼ 2: 2 2
Since y2 < y1 , we get y2 < 2. Subcase (ii) a < 0.
y2 ¼
a
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi a2 4d þ 8 a ðjaj 4Þ ¼ 2: > 2 2
Since y2 < y1 , we get y1 > 2.
S.S. Nemani / Applied Mathematics and Computation 215 (2010) 3830–3838
3833
Thus combine all cases, we have jyi j > 2 for i ¼ 1; 2. Since x2 yi x þ 1 ¼ 0 for i ¼ 1; 2 applying quadratic formula gives us
xi; ¼
yi
qffiffiffiffiffiffiffiffiffiffiffiffiffiffi y2i 4 2
:
Since jyi j > 2 for i ¼ 1; 2, the roots xi; are real and distinct. For yi > 2; xi;þ > 1 and jxi; j < 1. For yi < 2; xi; < 1 and jxi;þ j < 1. h From Theorem 2.1 we know that xi; xi;þ ¼ 1. This gives us that r 4 ¼ r 1 and r 3 ¼ r 1 1 2 . Let us assume that r; s 2 fri : 1 6 i 6 4g be such that jr þ sj þ jrsj < 1. (Existence of such r and s for some cases and more information about the roots are described in Appendix.) This means that rs – 1. Therefore, fr i : 1 6 i 6 4g ¼ fr; s; r 1 ; s1 g. Hence we have
f ðxÞ ¼ ðx rÞðx sÞðx s1 Þðx r 1 Þ; which is equal to
x4 ðr þ s þ s1 þ r 1 Þx3 þ ðrs þ rs1 þ 2 þ sr1 þ s1 r 1 Þx2 ðr þ s þ s1 þ r 1 Þx þ 1: Hence comparing the coefficients we get
r þ s þ s1 þ r 1 ¼ a
ð2:12Þ
rs þ rs1 þ 2 þ sr 1 þ s1 r1 ¼ d:
ð2:13Þ
and
Thus, if we set
l1 ¼ ðr þ sÞ;
l2 ¼ rs;
u2 ¼ ðr 1 þ s1 Þ and u1 ¼ r1 s1 :
Then using Eqs. (2.12) and (2.13), it can be shown that l1 ; l2 ; u1 and u2 satisfy Eqs. (2.4)–(2.7). Note that jl1 j þ jl2 j ¼ jr þ sj þ jrsj < 1. Dividing by jrsj both sides in this inequality, we get ju1 j > ju2 j þ 1 and hence both matrices L and U are strictly diagonally dominant. The computation of li and ui for i ¼ 1; 2 gives rise to a Toeplitz LU factorization of the matrix A0 in Oð1Þ time. To solve Eq. (2.1), we first solve A0 x0 ¼ b. To find x0 , we first solve Ly ¼ b and then solve Ux0 ¼ y by forward and backward substitution, respectively. To obtain the vector x0 , it takes 8n operations. The diagonally dominant property of matrices L and U guarantees the numerical stability of the process of solving equation A0 x0 ¼ b. To estimate kx0 k we need the followings. Lemma 2.1. Let L be the matrix as given in Eq. (2.2) and jl1 j þ jl2 j < 1. If Ly ¼ b, then
kyk 6
kbk : ð1 jl1 j jl2 jÞ
Proof Case (i) kyk ¼ jyi j for some i; 3 6 i 6 n. Since Ly ¼ b we get l2 yi2 þ l1 yi1 þ yi ¼ bi . Therefore,
kyk ¼ jyi j 6 jbi j þ jl2 jjyi2 j þ jl1 jjyi1 j 6 kbk þ jl2 jkyk þ jl1 jkyk: Thus
kyk 6
kbk : ð1 jl1 j jl2 jÞ
Case (ii) kyk ¼ jy1 j. Therefore,
kyk ¼ jy1 j ¼ jb1 j 6 kbk: Case (iii) kyk ¼ jy2 j, we have
kyk ¼ jy2 j ¼ jb2 jl1 y1 j 6 kbk þ jl1 jkyk: Thus
kyk 6
kbk : ð1 jl1 jÞ
The result follows from the fact that 1 jl1 j jl2 j 6 1 jl1 j 6 1.
h
Similarly, using the fact that ju1 j ju2 j 1 6 ju1 j ju2 j 6 ju1 j, we get the following lemma.
3834
S.S. Nemani / Applied Mathematics and Computation 215 (2010) 3830–3838
Lemma 2.2. Let U be the matrix as given in Eq. (2.3) and ju1 j > ju2 j þ 1. If Ux0 ¼ y, then
kx0 k 6
kyk : ðju1 j ju2 j 1Þ
Now we are ready to calculate a bound for kx0 k. It is given in the following theorem. Theorem 2.2. Let A0 be the matrix as given in Eq. (2.8) . If A0 x0 ¼ b, then
kx0 k 6
kbk : ð1 jl1 j jl2 jÞðju1 j ju2 j 1Þ
Proof. Since x0 is solution of A0 x0 ¼ b; x0 is solution of Ux0 ¼ y, where y is solution of Ly ¼ b. Hence from Lemmas 2.2 and 2.1 we get
kx0 k 6
kyk kbk 6 : ðju1 j ju2 j 1Þ ð1 jl1 j jl2 jÞðju1 j ju2 j 1Þ
3. Correction procedure Note that A ¼ A0 þ C implies that
Ax0 ¼ b þ Cx0 ; 0
where Cx ¼ ½ðl2 þ 0
l1 u2 Þx01
1
ð3:1Þ þ
l1 x02 e1
þ
½ðl2 u2 Þx01
þ
l2 x02 e2 .
We get the solution to Ax ¼ b as follows:
1
x ¼ x aA e1 bA e2 ; where a ¼ ðl2 þ l1 u2 Þx01 þ l1 x02 and b ¼ ðl2 u2 Þx01 þ l2 x02 . If we can find vectors wðiÞ such that AwðiÞ ¼ ei for i ¼ 1; 2 then we can solve Ax ¼ b approximately by adding a linear combination of vectors wð1Þ and wð2Þ to the vector x0 . P Let uðiÞ ¼ tj¼1 r ji ej for i ¼ 1; 2 where t is to be chosen based on a predefined tolerance g. Using the fact that r1 and r 2 are 2 roots of the polynomial x4 þ ax3 þ dx þ ax þ 1, we get for i ¼ 1; 2 tþ1 AuðiÞ ¼ ðr1 et1 ðar tþ1 þ rtþ2 Þet þ ðr t1 þ ar ti Þetþ1 þ r ti etþ2 : i þ aÞe1 e2 r i i i i r2 r2 ð2Þ ð1Þ Let wð1Þ ¼ rr11r ðuð1Þ uð2Þ Þ and wð2Þ ¼ rr11r ððr1 ðr1 1 þ aÞu 2 þ aÞu Þ. 2 2 Then
Awð1Þ ¼ e1 þ
r1 r2 X p ej r 1 r 2 j2C j
ð3:2Þ
Awð2Þ ¼ e2 þ
r1 r2 X q ej ; r 1 r 2 j2C j
ð3:3Þ
and
where C ¼ ft 1; t; t þ 1; t þ 2g,
pt1 ¼ r tþ1 rtþ1 2 1 ; pt ¼ ða þ r 2 Þr tþ1 ða þ r 1 Þr 1tþ1 ; 2 t 1 t ptþ1 ¼ ðr 1 1 þ aÞr 1 ðr 2 þ aÞr 2 ;
ptþ2 ¼ r t1 r t2 and tþ1 tþ1 qt1 ¼ ðr 1 ðr1 2 þ aÞr 1 1 þ aÞr 2 ; tþ1 tþ1 qt ¼ ðr 1 ðr1 2 þ aÞða þ r 1 Þr 1 1 þ aÞða þ r 2 Þr 2 ; 1 t 1 1 t qtþ1 ¼ ðr 1 1 þ aÞðr 2 þ aÞr 2 ðr 2 þ aÞðr 1 þ aÞr 1 ; t 1 t qtþ2 ¼ ðr 1 1 þ aÞr 2 ðr 2 þ aÞr 1 :
Let
~ ¼ x0 awð1Þ bwð2Þ : x
ð3:4Þ
S.S. Nemani / Applied Mathematics and Computation 215 (2010) 3830–3838
3835
~ ¼ Ax0 aAwð1Þ bAwð2Þ . From Eqs. (3.1)–(3.3) we get Then Ax
~b¼ Ax
r1 r2 X ðapj þ bqj Þej : r 1 r 2 j2C
After rearranging the terms, we get
~b¼ Ax
r 1 r2 X ð1Þ t ð2Þ ðd r þ dj r t2 Þej ; r 1 r 2 j2C j 1
ð3:5Þ
where ð1Þ
dt1 ¼ fbðr 1 2 þ aÞ agr 1 ; ð1Þ
dt ¼ fbðr 1 2 þ aÞ agða þ r 1 Þr 1 ; ð1Þ
1 dtþ1 ¼ fa bðr 1 2 þ aÞgðr 1 þ aÞ; ð1Þ
dtþ2 ¼ a bðr 1 2 þ aÞ and ð2Þ
dt1 ¼ fa bðr 1 1 þ aÞgr 2 ; ð2Þ
dt ¼ fa bðr 1 1 þ aÞgða þ r 2 Þr 2 ; ð2Þ
1 dtþ1 ¼ fbðr 1 1 þ aÞ agðr 2 þ aÞ; ð2Þ
dtþ2 ¼ bðr1 1 þ aÞ a: 1 From Theorem 2.1 we know that r 1 1 ¼ r 4 ; r 2 ¼ r 3 and jr 1 j < jr 2 j < 1. Therefore, for all j 2 C,
ð1Þ
jdj j 6 ðjaj þ jbjjr 3 þ ajÞ maxf1; ja þ r 1 j; jr4 þ ajg and ð2Þ
jdj j 6 ðjaj þ jbjjr 4 þ ajÞ maxf1; ja þ r 2 j; jr3 þ ajg: Recall that a ¼ ðl2 þ l1 u2 Þx01 þ l1 x02 and b ¼ ðl2 u2 Þx01 þ l2 x02 . Using triangle inequalty and the definition of the max norm we obtain
jaj 6 ðjl1 j þ jl2 j þ jl1 jju2 jÞkx0 k and
jbj 6 jl2 jðju2 j þ 1Þkx0 k: Then using Eq. (3.5) we get
~ bk 6 kAx
jr 1 r 2 j fM1 þ M 2 gjr 2 jt kx0 k; jr1 r2 j
where M1 ¼ fðjl1 j þ jl2 j þ jl1 jju2 jÞ þ jl2 jðju2 j þ 1Þjr3 þ ajg maxf1; ja þ r 1 j; jr 4 þ ajg and M 2 ¼ fðjl1 j þ jl2 j þ jl1 jju2 jÞ þ jl2 jðju2 j þ 1Þ jr 4 þ ajg maxf1; ja þ r2 j; jr3 þ ajg. Using the bound for kx0 k from Theorem 2.2, we have the following theorem. ~ be as given in Eq. (3.4) . Then Theorem 3.1. Let x
~ bk 6 Mjr 2 jtþ2 kbk; kAx where
M¼
M1 þ M2 : jr 1 r 2 jð1 jl1 j jl2 jÞðju1 j ju2 j 1Þ
Note that jr2 j < 1 so the error decreases exponentially with increasing value of t. The number of components t that needs to be corrected is determined in relation to a specific tolerance g for the desired accuracy. For example, if we required ~ bk 6 gkbk, then from Theorem 3.1, t must be chosen to satisfy kAx
tðgÞ P
log g log M 2: log jr2 j
Substituting the vectors wð1Þ and wð2Þ in Eq. (3.4), we have
~ ¼ x0 x
r1 r2 r1 r2 ða bðr3 þ aÞÞuð1Þ þ ða bðr4 þ aÞÞuð2Þ : r1 r2 r1 r2
ð3:6Þ
3836
S.S. Nemani / Applied Mathematics and Computation 215 (2010) 3830–3838
Table 1 Algorithm for the method. Input
a and d so that A is as given in Eq. (2.1) The size of matrix A ¼ n, the desired tolerance g The right hand side vector b (1) r1 ; r2 ; r3 and r 4 as given in Theorem 2.1 (2) Solve Ly ¼ b, using forward substitution (3) Solve Ux0 ¼ y, using backward substitution l m glog M t ¼ loglog jr2 j
Calculate
(4) Do correction as given in Eq. (3.6) ~ to Ax ¼ b An approximate solution x
Output
Since jr1 j and jr2 j are less than 1, computing the value of ~ xi for increasing value of i will ensure the numerical stability of the ~ from the vector x0 , we need only 4t operations. Hence the total computational efmethod. Note that to obtain the vector x forts for the method given in this paper is 8n þ 4t þ Oð1Þ operations. For sufficiently large value of n, we get t n, implying this method takes 8n þ Oð1Þ operations. The algorithm for this method is given in Table 1. 4. Numerical result Table 2 provides feedback on the number of terms of t to be corrected for a given desired accuracy in the result. The values ~ bk, error, of t were computed with d ¼ ð2jaj þ 2 þ DÞ for different values of a and D. Table 3 presents the values of kAx with the right hand side vector b ¼ ð1; 1; . . . ; 1ÞT and n ¼ 2000. Note that the tolerance g cannot be achieved when tðgÞ > n, and hence the algorithm will break down when the required relative tolerance is sufficiently small and/or jdj is sufficiently close to 2jaj þ 2 (that is, D is sufficiently close to 0). Table 2 indicate that the value of t for large values of D is less ~ bk for large values of D is than the value of t for small values of D. Also note that from Table 3 it is clear that the value of kAx less than the error for small values of D. 5. Remark The method described in this paper is very efficient and stable one. It compares favorably with the regular LU method. Also this method is competitive with other methods [2] in term of arithmetic operations. Following the parallel approach described in [6], parallel implementation can be done in a similar way. Table 2 t-Values for the method.
D
g 1.0000E01
1.0000E03
1.0000E05
1.0000E07
For a ¼ 1 0.001 0.01 0.05 0.1 0.5 2 4 6
831 211 78 50 17 6 4 3
1156 314 124 83 32 14 9 7
1482 417 170 115 47 22 15 12
1808 520 216 148 61 29 21 17
For a ¼ 0:5 0.001 0.01 0.05 0.1 0.5 2 4 6
747 187 68 43 14 5 3 2
1056 285 112 74 28 12 8 7
1365 382 155 105 42 20 14 11
1674 480 199 136 56 27 19 16
For a ¼ 0:1 0.001 0.01 0.05 0.1 0.5 2 4 6
664 163 58 36 12 4 2 1
959 256 100 66 25 11 7 6
1254 349 141 95 38 18 13 10
1549 443 183 125 52 25 18 15
S.S. Nemani / Applied Mathematics and Computation 215 (2010) 3830–3838
3837
Table 3 ~ bk. The values of kAx
g
D
1.0000E01
1.0000E03
1.0000E05
1.0000E07
For a ¼ 1 0.001 0.01 0.05 0.1 0.5 2 4 6
7.8134E03 7.7988E03 7.8132E03 7.9590E03 8.1115E03 8.5965E03 4.7283E03 4.7218E03
7.8849E05 7.8002E05 7.9049E05 7.5814E05 7.5373E05 7.0900E05 8.6140E05 9.2079E05
7.8453E07 7.8016E07 7.9976E07 8.3154E07 7.0038E07 5.7866E07 6.4396E07 7.2304E07
7.8052E09 7.8028E09 8.0914E09 7.9209E09 8.8902E09 8.6133E09 4.8447E09 5.9141E09
For a ¼ 0:5 0.001 0.01 0.05 0.1 0.5 2 4 6
1.0859E02 1.0896E02 1.1051E02 1.1589E02 1.2336E02 1.0927E02 9.0311E03 5.2017E03
1.0849E04 1.0755E04 1.0778E04 1.1584E04 1.2520E04 1.2457E04 1.0229E04 5.7487E05
1.0838E06 1.1128E06 1.1678E06 1.1578E06 1.2701E06 8.2621E07 6.6942E07 1.0113E06
1.0827E08 1.0984E08 1.1390E08 1.1572E08 1.2885E08 1.0194E08 9.9379E09 6.4563E09
For a ¼ 0:1 0.001 0.01 0.05 0.1 0.5 2 4 6
1.7136E02 1.7205E02 1.7422E02 1.8725E02 1.3863E02 9.3424E03 7.3716E03 2.0686E02
1.7106E04 1.7449E04 1.7013E04 1.7609E04 1.7138E04 1.5216E04 1.9746E04 5.2503E05
1.7076E06 1.7696E06 1.8549E06 1.9345E06 1.9696E06 1.2285E06 9.8031E07 9.5138E07
1.7045E08 1.7082E08 1.8114E08 1.8190E08 1.6297E08 1.5108E08 8.3295E09 1.0580E08
Appendix Case (i) d > 0 and a > 0. From Theorem 2.1, we have y2 < y1 < 2. Therefore,
qffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffi y21 4 < y22 4:
Hence x2; < x1; < 0 means jx1; j < jx2; j. Since xi;þ xi; ¼ 1 for i ¼ 1; 2, we get that the signs of xi;þ and xi; are same. Thus, all the roots are negative. Therefore, in this case we will have r 1 ¼ x1; ; r2 ¼ x2; ; r3 ¼ x2;þ and r 4 ¼ x1;þ . Case (ii) d > 0 and a < 0. In this case, 2 < y2 < y1 , therefore, x2;þ < x1;þ . The sign of all the roots in this case is positive. In this case we have r1 ¼ x1; ; r 2 ¼ x2; ; r 3 ¼ x2;þ and r 4 ¼ x1;þ . Case (iii) d < 0. In this case, y1 > 2 and y2 < 2. Thus x1; > 0 and x2; < 0. Subcase (i) jx1; j < jx2;þ j. In this case, we have r1 ¼ x1; ; r 2 ¼ x2;þ ; r 3 ¼ x2; and r4 ¼ x1;þ . Then r1 < r 2 . Therefore, ðr 1 þ r 2 Þ > 0 and r 1 r2 < 0. Consider 1 jr 1 þ r 2 j jr1 r 2 j ¼ 1 þ r1 þ r2 þ r 1 r 2 ¼ ð1 þ r1 Þð1 þ r 2 Þ. Since r 1 > 0; r1 þ 1 > 0 and since 1 < r 2 < 0 we get 1 þ r 2 > 0. Thus, 1 jr1 þ r2 j jr 1 r2 j > 0. Subcase (ii) jx1; j > jx2;þ j. In this case, we have r 1 ¼ x2;þ ; r 2 ¼ x1; ; r3 ¼ x1;þ and r 4 ¼ x2; . Then r 1 < r 2 . Therefore, ðr1 þ r 2 Þ > 0 and r 1 r2 < 0. Consider 1 jr1 þ r 2 j jr 1 r 2 j ¼ 1 r 1 r2 þ r1 r2 ¼ ð1 r1 Þð1 r2 Þ. Since r1 > 0; 1 r 1 > 0 and since 0 < r2 < 1 we get 1 r2 > 0. Thus, 1 jr 1 þ r 2 j jr1 r 2 j > 0.Combine both subcases, we get jr 1 þ r 2 j þ jr 1 r 2 j < 1. Thus take r ¼ r 1 and s ¼ r 2 . References [1] J.H. Ahlberg, E.N. Nilson, J.L. Walsh, The Theory of Splines and their Applications, Academic Press, New York, 1967. [2] F. Diele, L. Lopez, The use of the factorization of five-diagonal matrices by tridiagonal Toeplitz matrices, Appl. Math. Lett. 11 (1998) 51–59. [3] P.G. Patil, Y.S. Swamy, An efficient model for vibration control by piezoelectric smart structure using finite element method, Eur. J. Comput. Sci. Netw. Secu. 8 (11) (2008) 258–264. [4] O. Rojo, A new method for solving symmetric circulant tri-diagonal system of linear equations, Comput. Math. Appl. 20 (1990) 61–67.
3838
S.S. Nemani / Applied Mathematics and Computation 215 (2010) 3830–3838
[5] S.S. Nemani, L.E. Garey, An efficient method for solving second order boundary value problems with two point boundary conditions, Int. J. Comput. Math. 79 (9) (2002) 1001–1008. [6] S.S. Nemani, L.E. Garey, Parallel algorithms for solving tridiagonal and near-circulant systems, Appl. Math. Comput. 130 (2002) 285–294. [7] E.W. Sachs, A.K. Strauss, Efficient solution of a partial integro-differential equation in finance, Appl. Numer. Math. 58 (2008) 1687–1703. [8] P. Whittle, Prediction and Regulation by Linear Least Squares Methods, Van Nostrand, Princeton, NJ, 1963. [9] W.M. Yan, K.L. Chung, A fast algorithm for solving special tri-diagonal systems, Computing 52 (1994) 203–211.