A new iterative method for solving linear systems

A new iterative method for solving linear systems

Applied Mathematics and Computation 179 (2006) 725–730 www.elsevier.com/locate/amc A new iterative method for solving linear systems q Nenad Ujevic...

128KB Sizes 5 Downloads 190 Views

Applied Mathematics and Computation 179 (2006) 725–730 www.elsevier.com/locate/amc

A new iterative method for solving linear systems

q

Nenad Ujevic´ Department of Mathematics, University of Split, Teslina 12/III, 21000 Split, Croatia

Abstract A new iterative method for solving linear systems is derived. It can be considered as a modification of the Gauss–Seidel method. The modified method can be two times faster than the original one. Ó 2006 Published by Elsevier Inc. Keywords: Linear system; Iterative method; Gauss–Seidel method; Modification

1. Introduction In recent years a number of authors have considered iterative methods for solving linear systems. In Refs. [1–7] we mentioned only few of published articles. In this paper, we consider the linear system of equations Ax = b, where A is a positive definite matrix of order n and b 2 Rn is a given element. As we know, such a system can be solved by the well-known Gauss– Seidel method. This method can be derived such that we seek an approximate solution of the problem 1 f ðxÞ ¼ ðAx; xÞ  ðb; xÞ ! inf . 2 We form the sequence xkþ1 ¼ xk þ aki ei ;

k ¼ 0; 1; 2; . . . ; i ¼ 1; 2; . . . ; n;

where e1 = (1, 0, 0, . . . , 0), e2 = (0, 1, 0, . . . , 0), . . . , en = (0, 0, . . . , 0, 1), aki 2 R and x0 2 Rn is a given element. The parameters aki are chosen such that the method converges. In such a way, the Gauss–Seidel method examine equations of the system Ax = b one at a time in sequence and previously computed results are used as soon as they are available. Here, we give a new iterative method for solving linear systems. It can be considered as a modification of the Gauss–Seidel method. The modified method updates two components of the approximate solution xk at the same time. In the modification of the Gauss–Seidel method, considered in this paper, a sequence of approximate solutions is formed by the formula q

This paper is in final form and no version of it will be submitted for publication elsewhere. E-mail address: [email protected]

0096-3003/$ - see front matter Ó 2006 Published by Elsevier Inc. doi:10.1016/j.amc.2005.11.128

N. Ujevic´ / Applied Mathematics and Computation 179 (2006) 725–730

726

xkþ1 ¼ xk þ aki ei þ cki ej ; where cki 2 R and j depends on i. As we know, the Gauss–Seidel method requires 2n2 arithmetic operations per one iteration cycle. The modified method requires only 4n arithmetic operations more than the Gauss– Seidel method, per one iteration cycle (it does not require 2n2 + 2n2 operations, as we could expect). At the same time, it can be two times faster than the original method, since it updates two components of the approximate solution. Besides, we consider a general method, i.e. we form the sequence xkþ1 ¼ xk þ aki ei þ cki qj ; where qj 2 Rn is an arbitrary element. We also show that this general method gives a better reduction of the error xk  x* (Ax* = b) than the Gauss–Seidel method. Of course, the above mentioned modification also give a better reduction of the error. 2. Preliminary results We consider the linear system of equations: Ax ¼ b;

ð2:1Þ n

where A is a symmetric positive definite matrix of order n and b 2 R is a given element. We also consider the problem: 1 f ðxÞ ¼ ðAx; xÞ  ðb; xÞ ! inf : 2 0

ð2:2Þ

00

Since f (x) = Ax  b and f (x) = A, the system (2.1) and the problem (2.2) have the same solution x* 2 Rn. We seek an approximate solution of the problem (2.2). For that purpose, we form a sequence (xk) by the procedure: for

k ¼0; 1; 2; . . . y 1 ¼xk

for i ¼1; 2; . . . ; n y iþ1 ¼y i þ ai ei end for i xkþ1 ¼y nþ1 stopping criteria end

for k

where x0 2 Rn is a given element, ai are chosen in a given way and e1 = (1, 0, 0, . . . , 0), e2 = (0, 1, 0, . . . , 0), . . . , en = (0, 0, . . . , 0, 1). If we choose p ai ¼  i ; pi ¼ ðai ; y i Þ  bi ; i ¼ 1; 2; . . . ; n ð2:3Þ aii then we get the Gauss–Seidel method. We have 1 f ðy þ hÞ ¼ f ðyÞ þ ðf 0 ðyÞ; hÞ þ ðAh; hÞ 2 if the function f is given by (2.2). Substituting y = yi, h = aiei in (2.4) we get 1 1 f ðy i þ ai ei Þ ¼ f ðy i Þ þ ai ðAy i  b; ei Þ þ a2i ðAei ; ei Þ ¼ f ðy i Þ þ ai pi þ a2i aii . 2 2

ð2:4Þ

ð2:5Þ

N. Ujevic´ / Applied Mathematics and Computation 179 (2006) 725–730

727

From (2.3) and (2.5) we get f ðy iþ1 Þ  f ðy i Þ ¼ 

1 p2i 6 0; 2 aii

ð2:6Þ

i.e. f(yi+1) 6 f(yi). On the other hand, we can show that the reduction of f is equivalent to the reduction of the error gi = yi  x*. Let us do that. If A is a symmetric positive definite matrix of order n then A defines a new scalar product in Rn by the formula hx; yi ¼ ðAx; yÞ;

x; y 2 Rn .

The corresponding norm is kxk2A ¼ ðAx; xÞ;

x 2 Rn .

We have ky iþ1  x k2A  ky i  x k2A ¼ ðAy iþ1  Ax ; y iþ1  x Þ  ðAy i  Ax ; y i  x Þ; ¼ ðAy iþ1 ; y iþ1 Þ  2ðb; y iþ1 Þ  ½ðAy i ; y i Þ  2ðb; y i Þ ¼ 2f ðy iþ1 Þ  2f ðy i Þ. Hence, the reduction of the function f is equivalent to the reduction of the error in the A-norm. We have 2

2

kgiþ1 kA  kgi kA ¼ 

p2i 6 0; aii

ð2:7Þ

i.e. kgi+1kA 6 kgikA. A consequence is the convergence of Gauss–Seidel method. 3. A general method Let y iþ1 ¼ y i þ hi ;

hi ¼ 

pi ei ; aii

ð3:1Þ

where pi ¼ ðai ; y i Þ  bi .

ð3:2Þ

Let us consider the element ziþ1 ¼ y i þ hi þ ci qi ;

ð3:3Þ

where ci 2 R, qi 2 Rn. From (2.4) we get 1 f ðy i þ hi þ ci qi Þ ¼ f ðy i Þ þ ðf 0 ðy i Þ; hi þ ci qi Þ þ ðAðhi þ ci qi Þ; hi þ ci qi Þ 2 1 1 ¼ f ðy i Þ þ ðAy i  b; hi Þ þ ci ðAy i  b; qi Þ þ ðAhi ; hi Þ þ ci ðAqi ; hi Þ þ c2i ðAqi ; qi Þ 2 2 1 2 ¼ f ðy iþ1 Þ þ ci ½ðAy i  b; qi Þ þ ðAqi ; hi Þ þ ci ðAqi ; qi Þ. 2

ð3:4Þ

We now define the function 1 gðcÞ ¼ c½ðAy i  b; qi Þ þ ðAqi ; hi Þ þ c2 ðAqi ; qi Þ 2

ð3:5Þ

N. Ujevic´ / Applied Mathematics and Computation 179 (2006) 725–730

728

such that g0 ðcÞ ¼ðAy i  b; qi Þ þ ðAqi ; hi Þ þ cðAqi ; qi Þ; g00 ðcÞ ¼ðAqi ; qi Þ > 0. From the equation g 0 (c) = 0 we get ci ¼ 

ðAy i  b; qi Þ þ ðAqi ; hi Þ . ðAqi ; qi Þ

ð3:6Þ

From (3.6) and (3.5) we have gðci Þ ¼  ¼

½ðAy i  b; qi Þ þ ðAqi ; hi Þ2 1 ½ðAy i  b; qi Þ þ ðAqi ; hi Þ2 þ ðAqi ; qi Þ 2 2 ðAqi ; qi Þ ðAqi ; qi Þ 1 2 ½ðAy i  b; qi Þ þ ðAqi ; hi Þ 6 0. 2ðAqi ; qi Þ

ð3:7Þ

From (3.3), (3.4) and (3.7) we get f ðziþ1 Þ  f ðy i Þ ¼ f ðy iþ1 Þ  f ðy i Þ 

1 2 ½ðAy i  b; qi Þ þ ðAqi ; hi Þ . 2ðAqi ; qi Þ

ð3:8Þ

From (3.8) we see that the element zi+1 gives a better reduction of the function f than the element yi+1. As we know, this means that the element zi+1 gives better reduction of the error than the element yi+1. From (3.8) and (2.7) we get kziþ1  x k2A  ky i  x k2A ¼

p2i aii



ð3:9Þ

1 2 ½ðAy i  b; qi Þ þ ðAqi ; hi Þ . ðAqi ; qi Þ

We now describe the modification of the Gauss–Seidel method. We form a sequence (xk) by the procedure for

k ¼0; 1; 2; . . . z1 ¼xk

for i ¼1; 2; . . . ; n ziþ1 ¼zi þ hi þ ci qi end

for i

xkþ1 ¼znþ1 end

for k

where ^ pi ; ^ pi ¼ ðai ; zi Þ  bi ; aii ðAzi  b; qi Þ þ ðAqi ; hi Þ . ci ¼  ðAqi ; qi Þ

hi ¼ 

ð3:10Þ ð3:11Þ

Note that qi can be any element of Rn. However, we choose the elements qi such that the modified method is effective in applications.

N. Ujevic´ / Applied Mathematics and Computation 179 (2006) 725–730

729

4. A particular method We choose: qi = ej, where j depends on i, j 5 i. Then we have ziþ1 ¼ zi þ hi þ ci ej ; ^ pi pi ¼ ðai ; zi Þ  bi ; hi ¼  ei ; ^ aii

ð4:1Þ ð4:2Þ ^ p

ci ¼ 

^ pj  aiii aij ðAzi  b; ej Þ þ ðAej ; hi Þ aii ^pj  ^pi ¼ ¼ ; ðAej ; ej Þ ajj ajj aii

ci ¼ 

^ aij pj  ^pi ; ajj aii ajj

i.e. ð4:3Þ

where ^ pj ¼ ðaj ; zi Þ  bj .

ð4:4Þ

We can choose the indices j in different ways. Here, we choose j = i  1, if i = 2, 3, . . . , n and j = n if i = 1. Then we have zi ¼ zi1 þ hi1 þ ci1 ei2 ;

hi1 ¼ 

^ pi1 ai1;i1

ð4:5Þ

ei1 ;

where c0 = cn, e0 = en, e1 = en1 and h0 = hn, ^ p0 ¼ ^pn , a00 = ann. From (4.4) and (4.5) we get ^ pj ¼ ðaj ; zi Þ  bj ¼ ðai1 ; zi1 þ hi1 þ ci1 ei2 Þ  bi1 ¼ ðai1 ; zi1 Þ  bi1 

^pi1 ai1;i1

ðai1 ; ei1 Þ þ ci1 ðai1 ; ei2 Þ ¼ ^pi1 

¼ ci1 ai1;i2 ;

^pi1 ai1;i1

ai1;i1 þ ci1 ai1;i2 ð4:6Þ

where a0,1 = an,n1, a1,0 = a1,n. From (4.3) and (4.6) we get ci ¼ ci1

^ ai1;i2 pi1 þ ai;i1 . ai1;i1 ai1;i1 aii

The above results provide the next algorithm. Algorithm 1. Choose x0 = (x1, . . . , xn) 2 Rn and set x = x0. for i = 1, 2, . . . , n j=i1 m=i2 if i = 1 then m=n1 j=n endif if i = 2 then m = n ajm ri ¼  aijajj ti ¼ aii ajj end for i pn = (an,x)  bn p1 = (a1,x)  b1

ð4:7Þ

730

N. Ujevic´ / Applied Mathematics and Computation 179 (2006) 725–730

pn a1n þ p1 ann a11 ann for k = 0, 1, . . . for i = 1, 2, . . . , n j=i1 if i = 1 then j = n pi = (ai,x)  bi p xi ¼ xi  i aii ci = cjri + piti xj = xj + ci end for i stopping criteria end for k cn ¼ 

We must note that this modification of the Gauss–Seidel method requires only four arithmetic operations (additions and multiplications) per one iteration step more than the original Gauss–Seidel method. On the other hand, it can be two times faster than the Gauss–Seidel method as we can see if we solve the next example. Example 1. Let the matrix A be given by aii ¼ 4n; ai;iþ1 ¼ aiþ1;i ¼ n; aij ¼ 0:5; j 6¼ i; i þ 1. Pn Let b ¼ k¼1 aik such that x = (1, 1, . . . , 1) is the exact solution of the system Ax = b. If we choose n = 1000, the initial guess x0 = (x1, . . .,xn), xi = 0.001 * i, i = 1, 2, . . . , n, the stopping criteria kxk+1  xkk < 0.000001 and solve the above problem by the Gauss–Seidel method and by the modified method then we shall see that the modified method is two times faster than the Gauss–Seidel method. In fact, the Gauss–Seidel method stops after k = 11 iterations and the modified method stops after k = 5 iterations. The time of execution for the Gauss–Seidel method is 1.42 s and for the modified method it is 0.632 s. References [1] G.H. Golub, C.F. Van Loan, Matrix Computations, The Johns Hopkins University Press, Baltimore and London, 1996. [2] A. Hadjidimos, Successive overrelaxation (SOR) and related methods, J. Comput. Appl. Math. 123 (2000) 177–199. [3] W. Li, The convergence of the modified Gauss–Seidel methods for consistent linear systems, J. Comput. Appl. Math. 154 (2003) 97– 105. [4] W. Li, A note on the preconditioned Gauss–Seidel (GS) method for linear systems, J. Comput. Appl. Math. 182 (2005) 81–90. [5] W. Li, W. Sun, Modified Gauss–Seidel type methods and Jacobi type methods for Z-matrices, Lin. Alg. Appl. 317 (2000) 227–240. [6] H. Niki, K. Harada, M. Morimoto, M. Sakakihara, The survey of preconditioners used for accelerating the rate of convergence in the Gauss–Seidel method, J. Comput. Appl. Math. 164–165 (2004) 587–600. ¨ zban, Improved convergence criteria for Jacobi and Gauss–Seidel iterations, Appl. Math. Comput. 152 (3) (2004) 693–700. [7] A.Y. O