A computational algorithm for solving periodic penta-diagonal linear systems

A computational algorithm for solving periodic penta-diagonal linear systems

Applied Mathematics and Computation 174 (2006) 613–618 www.elsevier.com/locate/amc A computational algorithm for solving periodic penta-diagonal line...

112KB Sizes 0 Downloads 85 Views

Applied Mathematics and Computation 174 (2006) 613–618 www.elsevier.com/locate/amc

A computational algorithm for solving periodic penta-diagonal linear systems A.A. Karawia Department of Mathematics, Faculty of Science, Mansoura University, Mansoura 35516, Egypt

Abstract In this paper we present a new efficient computational algorithm for solving periodic penta-digonal linear systems. The implementation of the algorithm using Computer Algebra Systems (CAS) such as MAPLE, MACSYMA, MATHEMATICA, and MATLAB is straightforward. An example is given in order to illustrate the algorithm.  2005 Elsevier Inc. All rights reserved. Keywords: Periodic penta-diagonal matrices; Linear systems; Determinants; Computer algebra systems (CAS)

1. Introduction Many problems in mathematics and applied science require the solution of linear systems having penta-diagonal coefficient matrices, for example the solution of certain partial differential equations, spline approximation, etc. [1,2]. This kind of linear system is a special class of narrow banded linear systems.

E-mail address: [email protected] 0096-3003/$ - see front matter  2005 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2005.04.098

614

A.A. Karawia / Appl. Math. Comput. 174 (2006) 613–618

It often happens that a matrix has some favorable characteristics. These characteristics should be taken into account when we consider the development of new algorithms for solving problems that involve these kinds of matrices. In the current article, we are going to deal with one of these cases. In more detail, the main objective of the paper is to develop an efficient algorithm for solving general periodic penta-diagonal linear systems of the form: AX ¼ Y ;

ð1:1Þ

where 2



0

0

0

0 ~ a3

 0 .. . .. . .. .

0  .. . .. . .. .

0 0

bn2 ~n1 b 0

d n2

an2

bn1 ~ bn

d n1 bn

Y ¼ ðy 1 ; y 2 ; . . . ; y n ÞT

and

0 ~ a2

d1

a1

~ a1

6b 6 2 6~ 6 b3 6 6 60 6 6 . A¼6 6 .. 6 . 6 . 6 . 6 60 6 6 40 s

d2 b3 .. . .. .

a2 d3 .. . .. .

a3 .. .

..

.

..

.

..

.

 

..

.

..

0



0 0

 

X ¼ ðx1 ; x2 ; . . . ; xn ÞT ;

0

. ~ bn2

 0  

0 .. . .. .

t

3

7 7 7 7 7 7 0 7 7 .. 7 7 . 7; .. 7 7 . 7 7 ~an2 7 7 7 an1 5 0 0

ð1:2Þ

dn n P 4.

This linear system of equations frequently appears in many scientific and engineering areas. A general n · n periodic penta-diagonal matrix A of the form (1.2) can be stored in 5n memory locations by using five vectors ~ a ¼ ð~ a1 ; ~ a2 ; . . . ; ~ an2 ; t; sÞ, a ¼ ða1 ; a2 ; . . . ; an2 ; an1 Þ, ~b ¼ ð~b3 ; ~b4 ; . . . ; ~bn Þ, b ¼ ðb2 ; b3 ; . . . ; bn Þ, and d = (d1, d2, . . . , dn). When considering the system (1.1) it is advantageous to introduce three additional vectors c = (c1, c2, . . . , cn), e = (e1, e2, . . . , en2), and f = (f2, f3, . . . , fn1). These vectors are related to the vectors a, ~ a, b, ~ b, and d. The current paper is organized as follows. In Section 2, a new computational algorithm for the solution of the system (1.1) is given. An illustrative example is presented in Section 3. In Section 4, a Conclusion is given.

2. A new computational algorithm In this section we are going to formulate a new computational algorithm for solving periodic penta-diagonal linear systems of the form (1.1). To do this we begin by considering the LU decomposition [3] of the matrix A in the form:

A.A. Karawia / Appl. Math. Comput. 174 (2006) 613–618

2

~a1 0  0 a2 ~a2 0  ~a3 d 3 a3 0 .. .. .. .. . . . . .. .. .. ..  . . . .    0 ~bn2 bn2 d n2       0 ~bn1 bn1 ~bn 0   0 1 0 0  6 f2 1 0  6 6 . 6 ~b =c f 1 .. 6 3 1 3 6 . . . ¼6 .. .. .. 6 0 6 . .. .. .. 6 . 6 . . . . 6 4 0    0 ~bn1 =cn3 h2    hn3 h1

d1 6 b2 6 6 ~b3 6 6 6 60 6 . 6 . 6 . 6 60 6 40 s 2

a1 d2 b3 .. .

0 0  0 ..

t 0 0

615

3

7 7 7 7 7 7 0 7 .. 7 7 . 7 7 ~an2 7 7 an1 5

. an2 d n1 bn dn 32   0 c1 6   07 76 0 76 6 07 76 ... 76 .. 76 . 0 76 0 6 .7 .. .. 76 . .. 76 0 . 76 fn1 1 0 54 0 0 hn2 hn1 1

a1 e1 ~

0

 .. . .. .

0

v1

3

7 a2 c2 e2 ~ 0 v2 7 7 .. 7 .. .. .. .. . 7 . . . . 7 7. .. ~n3 vn3 7  . cn3 en3 a 7 7       0 cn2 en2 vn2 7 7   0 0 cn1 vn1 5   0 0 0 cn

ð2:1Þ From (2.1) we obtain 8 if i ¼ 1 d > < 1 if i ¼ 2 ci ¼ d 2  f2 e1 > ~ : bi d i  fi ei1  ci2 ~ ai2 if i ¼ 3ð1Þn  1; 8 t if i ¼ 1 > > > > > f v if i ¼ 2 > 2 1 > > < b~i if i ¼ 3ð1Þn  3 vi ¼ ci2 vi2  fi vi1 > ~ > b n2 >~ an2  cn4 vn4  fn2 vn3 if i ¼ n  2 > > > > > : a  ~bn1 v  f v if i ¼ n  1; n1 n1 n2 cn3 n3 8s if i ¼ 1 > c1 > > e1 > >  c2 h1 if i ¼ 2 > > >  < ~a ei1 i2 if i ¼ 3ð1Þn  3 hi ¼  ci hi2 þ ci hi1 > > ~ > bn n4 >  ~acn2 hn4  ecn3 hn3 if i ¼ n  2 > > cn2 n2 > > : bn ~ en2 an3  cn1 hn3  cn1 hn2 if i ¼ n  1 cn1

ð2:2Þ

ð2:3Þ

ð2:4Þ

and cn ¼ d n 

n1 X i¼1

v i hi .

ð2:5Þ

616

A.A. Karawia / Appl. Math. Comput. 174 (2006) 613–618 ~

bi For the sake of programming it is convenient to set ki ¼ ci2 , li ¼ ~ai2 , ci ei1 bi i ¼ 3ð1Þn  1, and ai ¼ ci , bi ¼ ci1 , i ¼ 2ð1Þn  1. It is not difficult to prove that the LU decomposition (2.1) exists only if ci 5 0, i = 1(1)n  1. Moreover the periodic penta-diagonal linear system (1.1) possesses a unique solution if, in addition, cn 5 0. On the other hand, the Determinant of the matrix A is given by n Y detðAÞ ¼ ci ð2:6Þ i¼1

and this shown the importance of the vector c. We may now formulate the following result. Algorithm 2.1. To solve the general periodic penta-diagonal linear system (1.1), we may proceed as follows: step 1: Set c1 ¼ d 1 , f2 ¼ bc12 , e1 ¼ a1 , v1 ¼ t, h1 ¼ cs1 , c2 ¼ d 2  f2 e1 , v2 ¼ f2 v1 , h2 ¼  ec12 h1 , e2 ¼ a2  f2 ~ a1 . step 2: For i = 3, 4, . . . , n  2 Compute ~ bi bi ei2 fi ¼ ci1  ci2 ; ci1 e i ¼ ai  f i ~ ai1 , ~ bi ~ ci ¼ d i  fi ei1  ci2 ai2 . ~

n1 n1  bcn3 step 3: Set fn1 ¼ bcn2

en3 , cn2 ~

n1 ~ cn1 ¼ d n1  fn1 en2  bcn3 an3 .

step 4: For i = 3, 4, . . . , n  3 Compute ~ bi vi ¼ ci2 vi2  fi vi1 ;   ei1 h þ h hi ¼  ~ai2 . i2 i1 ci ci ~

n2 an2  bcn4 vn4  fn2 vn3 , step 5: Set vn2 ¼ ~ ~ bn1 vn1 ¼ an1  cn3 vn3  fn1 vn2 , ~ bn n4 hn2 ¼ cn2  ~acn2 hn4  ecn3 hn3 , n2

step 6: step 7: step 8: step 9:

bn n3  ~acn1 hn3  ecn2 hn2 , hn1 ¼ cn1 n1 Pn1 cn ¼ d n  i¼1 hi vi . Set z1 = y1, z2 = y2  f2z1. For i = 3, 4, . . . , n  1 Compute ~ bi zi ¼ y i  fi zi1  ci2 zi2 . P n1 Compute zn ¼ y n  i¼1 hi zi . Compute the solution vector x using xn1 vn2 xn n1 xn xn ¼ cznn , xn1 ¼ zn1cv , xn2 ¼ zn2 en2cn2 , n1 for i = n  3, n  4, . . . , 1 compute xi ¼ zi ei xiþ1 ~cai i xiþ2 vi xn .

A.A. Karawia / Appl. Math. Comput. 174 (2006) 613–618

617

The new algorithm will be referred to as KPENTA algorithm. For the choice ~ ai ¼ 0, i = 1(1)n  2 and ~ bi ¼ 0, i = 3(1)n, we obtained the PERTRI algorithm [4]. See also the best known algorithm [5]. In [6], Claerbout showed that the two-dimensional Laplacian operator, which appears in 3-D finite-difference migration, has the form of penta-diagonal matrix. If we choice s = t = 0, di = 4, i = 1(1)n and ai ¼ bi ¼ ~ ai ¼ ~ bi ¼ 1 for all possible values of i, we can obtain it. KPENTA algorithm for solving the general periodic pentadiagonal system (1.1) is generally preferable because the conditions ci 5 0, i = 1(1)n are sufficient for the validity of it. The advantage of the vector c is now clear.

3. An illustrative example In this section we are going to give an illustrative example Example 3.1. Solve by 2 1 2 3 6 2 3 4 6 6 6 2 2 1 6 6 0 5 1 6 6 6 0 0 1 6 6 0 0 0 6 6 6 0 0 0 6 6 0 0 0 6 6 4 0 0 0 2 0 0

the periodic penta-diagonal linear system of size 10 given 3

2

10

3

0

0

0

0

0

0

4

5 3

0 2

0 0

0 0

0 0

0 0

2 7

1 5

3 0 2 2

0 0

0 0

2

3

4 1

2

0

0 0

3 0

0 1

7 2

0 0

0 0

0 0

3 0

6 14 7 07 6 7 7 6 7 7 6 6 7 07 6 7 7 6 0 7 7 07 6 7 6 7 7 6 3 7 07 7X ¼ 6 7 6 4 7 7 07 6 7 6 7 7 6 3 7 7 07 6 7 6 4 7 7 27 6 7 6 7 7 4 1 5 15 1 2

8 1 3 4 2 4

5 3

ð3:1Þ

by using the KPENTA algorithm. Solution The application of the KPENTA algorithm gives: • c1 = 1, e1 = 2, v1 = 4, h1 = 2, f2 = 2, c2 = 1, v2 = 8, h2 = 4, e2 = 2 (Step 1).   438 4314 • ½f3 ; f4 ; f5 ; f6 ; f7 ; f8  ¼ 6;  115 ; 17 ; 389 ; 2203 ; 12973 ,  17 73 6246 686 18500  ½e3 ; e4 ; e5 ; e6 ; e7 ; e8  ¼ 33; 5 ; 62 ; 343 ;  2203 ;  56206 , 12973   248 343 2203 12973 67067 ½c3 ; c4 ; c5 ; c6 ; c7 ; c8  ¼ 5; 5 ;  62 ; 686 ; 2203 ; 12973 (Step 2). ; c9 ¼  16634 (Step 3). • f9 ¼ 81446 67067 67067

618

A.A. Karawia / Appl. Math. Comput. 174 (2006) 613–618

  • ½v3 ; v4 ; v5 ; v6 ; v7  ¼ 40; 48; 160 ;  340 ; 6600 , 31 343 2203  17 157  467 1071 ; 686 ;  2203 ; 12973 (Step 4). ½h3 ; h4 ; h5 ; h6 ; h7  ¼ 25 ; 124 , v9 ¼  142185 , h8 ¼ 66386 , h9 ¼  40441 , c10 ¼  105341 (Step 5). • v8 ¼ 17026 12973 67067 67067 8317 8317 • z1 = 10, z2 = 6 (Step 6).   • ½z3 ; z4 ; z5 ; z6 ; z7 ; z8 ; z9  ¼ 10; 8;  37 ; 2987 ; 3276 ; 27887 ;  158819 (Step 7). 31 686 2203 12973 67067 ; • z10 ¼  105341 8317

x10 ¼ 1;

x9 ¼ 1;

x8 ¼ 1 (Step 8).

• [x1, x2, x3, x4, x5, x6, x7, x8, x9, x10] = [1, 1, 1, 1, 1, 1, 1, 1, 1, 1] (Step 9). Also we obtained the determinant of the current matrix det(A) = 421364 by using (2.6).

4. Conclusion The method described here is a very effective, provided optimal LU factorization is used. Our method is competitive the other methods for solving circulant linear systems which appear in many applications.

Acknowledgement I should like to thank Prof. Dr. M. E. A. El-Mikkawy for several comments and suggestions.

References [1] M. Chen, On the solution of circulant linear systems, SIAM J. Numer. Anal. 24 (1987) 668–683. [2] H.L. Stone, Iterative solution of implicit approximations of multidimensional partial differential equations, SIAM J. Numer. Anal. 5 (3) (1968) 530–558. [3] M.B. Allen III, E.L. Isaacson, Numerical Analysis for Applied Science, Wiley-Interscience, John Wiley & Sons, 1997. [4] M.E.A. El-Mikkawy, A new computational algorithm for solving periodic tri-diagonal linear systems, Appl. Math. Comput. 161 (2) (2005) 691–696. [5] M. Chawla, R.R. Khazal, A parallel elimination method for periodic tri-diagonal systems, Int. J. Comput. Math. 79 (4) (2002) 473–484. [6] J.F. Claerbout, Imaging the Earths Interior, Blackwell Scientific Publications, 1985.