A modified simple iterative method for nonsymmetric algebraic Riccati equations arising in transport theory

A modified simple iterative method for nonsymmetric algebraic Riccati equations arising in transport theory

Applied Mathematics and Computation 181 (2006) 1499–1504 www.elsevier.com/locate/amc A modified simple iterative method for nonsymmetric algebraic Ric...

225KB Sizes 0 Downloads 40 Views

Applied Mathematics and Computation 181 (2006) 1499–1504 www.elsevier.com/locate/amc

A modified simple iterative method for nonsymmetric algebraic Riccati equations arising in transport theory q Liang Bao

a,b

, Yiqin Lin

c,* ,

Yimin Wei

a,b

a

b

Institute of Mathematics, School of Mathematical Sciences, Fudan University, Shanghai 200433, PR China Key Laboratory of Mathematics for Nonlinear Sciences (Ministry of Education), Fudan University, Shanghai 200433, PR China c School of Mathematics and Computational Science, Sun Yat-Sen University, Guangzhou 510275, PR China

Dedicated to Prof. Zhi-hao Cao on the occasion of his 70th birthday

Abstract In this paper we propose a modified simple iterative method for nonsymmetric algebraic Riccati equations arising in transport theory. It is similar to the method proposed by Lu in [L.Z. Lu, Solution form and simple iteration of a nonsymmetric algebraic Riccati equation arising in transport theory, SIAM J. Matrix Anal. Appl. 26 (2005) 679–685]. The theoretical results and numerical experiments show that the modified version is more efficient than its counterpart.  2006 Elsevier Inc. All rights reserved. Keywords: Nonsymmetric algebraic Riccati equations; Simple iterative method; Minimal positive solution

1. Introduction In this paper, we are interested in iteratively solving the following algebraic Riccati equation arising in transport theory (see [2,3] and the references therein): XCX  XE  AX þ B ¼ 0;

ð1Þ

where A; B; C; E 2 Rnn are given by A ¼ D  eqT ;

B ¼ eeT ;

C ¼ qqT ;

E ¼ D  qeT .

ci , Here e = (1, 1, . . . , 1)T, q = (q1, q2, . . . , qn)T with qi ¼ 2x i

q Y. Lin is supported by Scientific Research Startup Foundation of Sun Yat-Sen University for Young Teacher under Grant 2005-340001131042. Y. Wei is supported by the National Natural Science Foundation of China under Grant 10471027 and Shanghai Education Committee. * Corresponding author. E-mail addresses: [email protected] (L. Bao), [email protected] (Y. Lin), [email protected] (Y. Wei).

0096-3003/$ - see front matter  2006 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2006.02.042

1500

L. Bao et al. / Applied Mathematics and Computation 181 (2006) 1499–1504

(

1 D ¼ diagðd1 ; d2 ; . . . ; dn Þ with di ¼ cxi ð1þaÞ ; 1 D ¼ diagðd 1 ; d 2 ; . . . ; d n Þ with d i ¼ cxi ð1aÞ

and 0 < c 6 1, 0 6 a < 1, 0 < xn <    < x2 < x1 < 1, n X ci ¼ 1; ci > 0; i ¼ 1; 2; . . . ; n. i¼1

It has been shown in [2] and [3] that (1) has positive solutions (in the componentwise sense). Since only the minimal positive solution is physically meaningful, some iterative methods have been developed for computing the minimal positive solution of (1), see [2] and [3]. Recently, Lu [4] has shown that the matrix equation (1) is in essence equivalent to a ‘‘vector equation’’ and developed a simple but efficient iterative procedure to compute the minimal positive solution of (1). In this paper we present a modified version of this simple iterative procedure and prove that the modified simple iterative procedure is more efficient than the original one. This paper is organized as follows. In Section 2, we review the simple iterative procedure proposed by Lu [4] and present the modified simple iterative procedure. Some theoretical results are given in this section. Sections 3 and 4 give some numerical examples and conclusions, respectively. Let A = (ai,j), B = (bi,j) be n · n matrices with real entries. The Hadamard product of A and B is defined by A  B = (ai,j Æ bi,j). 2. Simple iteration and modified simple iteration We rewrite (1) as DX þ XD ¼ ðXq þ eÞðqT X þ eT Þ and let v T ¼ qT X þ e T .

u ¼ Xq þ e;

ð2Þ

It is shown in [2] and [3] that the positive solutions of (1) are of the following form: X ¼ T  ðuvT Þ ¼ ðuvT Þ  T ; where

 T ¼ ðti;j Þ ¼ X ¼ ðxi;j Þ;

ð3Þ

 1 ; di þ d j

u ¼ ðu1 ; u2 ; . . . ; un Þ

T

T

v ¼ ðv1 ; v2 ; . . . ; vn Þ .

and

To find the minimal positive solution of (1), we need to find proper positive vectors u and v in (3). Substituting (3) into (2), we obtain  u ¼ u  ðPvÞ þ e; ð4Þ v ¼ v  ðP~ uÞ þ e; where P ¼ ðpi;j Þ ¼



 qj ; di þ d j

P~ ¼ ð~ pi;j Þ ¼



 qj . dj þ d i

To get the positive vectors u and v from (4), Lu [4] has defined a simple iterative scheme 8 ðkþ1Þ > ¼ uðkÞ  ðPvðkÞ Þ þ e; : ð0Þ u ¼ vð0Þ ¼ 0 and the following results are given.

ð5Þ

L. Bao et al. / Applied Mathematics and Computation 181 (2006) 1499–1504

1501

Lemma 1 [4]. The minimal positive solution X* of (1) can be computed by X  ¼ T  ðu ðv ÞT Þ; where (u*, v*) is the limit of {(u(k), v(k))} defined by the simple iteration (5). Remark 2. Lu has showed that the sequence {(u(k), v(k))} is strictly monotonically increasing and bounded above. Before computing v(k+1), we have obtained u(k+1). It should be a better approximation to u* than u(k). We replace u(k) with u(k+1) in the second equation of (5), and obtain a modified simple iterative scheme 8 uðkþ1Þ ¼ ^ uðkÞ  ðP^vðkÞ Þ þ e; > <^ ð6Þ ^vðkþ1Þ ¼ ^vðkÞ  ðP~ ^ uðkþ1Þ Þ þ e; k ¼ 0; 1; . . . ; > : ð0Þ ð0Þ ^ u ¼ ^v ¼ 0: Then, we have the following result. Lemma 3. For all 0 < c 6 1 and 0 6 a < 1, the sequence fð^uðkÞ ; ^vðkÞ Þg defined by (6) is strictly monotonically increasing, bounded above and thus converging. Proof. It is easy to show by the induction that fð^ uðkÞ ; ^vðkÞ Þg is strictly monotonically increasing. To show that ðkÞ ðkÞ fð^ u ; ^v Þg is bounded above, we observe that ^uð0Þ < uð1Þ ;

^vð0Þ < vð1Þ

and

^ uð1Þ < uð2Þ ;

^vð1Þ < vð2Þ ;

^ð2Þ < uð3Þ ; ^vð2Þ < vð4Þ . Therefore, the statement it follows that u ^ uðiÞ < uð2i1Þ ;

^vðiÞ < vð2iÞ

ð7Þ

is true for i = 2. We now assume that (7) is true for i = k P 2. By (5) and (6) we have ^ uðkÞ  ðP^vðkÞ Þ þ e < uð2kÞ  ðPvð2kÞ Þ þ e ¼ uð2kþ1Þ uðkþ1Þ ¼ ^ and ^vðkþ1Þ ¼ ^vðkÞ  ðP~ ^ uðkþ1Þ Þ þ e < vð2kþ1Þ  ðP~ uð2kþ1Þ Þ þ e ¼ vð2kþ2Þ . We have thus proved that (7) is true for i = k + 1. Hence, by the principle of mathematical induction, (7) is true for all k P 2. Since {(u(k), v(k))} is bounded above (see [4]), so is fð^uðkÞ ; ^vðkÞ Þg. The sequence fð^uðkÞ ; ^vðkÞ Þg is now well defined, monotonically increasing, and bounded above. h Now we present our main result. Theorem 4. The minimal positive solution X* of (1) can be computed by X  ¼ T  ð^ u ð^v ÞT Þ; uðkÞ ; ^vðkÞ Þg defined by the modified simple iteration (6). where ð^ u ; ^v Þ is the limit of fð^ Proof. It is easy to show by the induction that uðkÞ < ^uðkÞ and vðkÞ < ^vðkÞ hold for all k P 3. By (7), we obtain uðkÞ < ^ uðkÞ < uð2k1Þ ;

vðkÞ < ^vðkÞ < vð2kÞ

ðk P 3Þ.

So the sequence fð^ uðkÞ ; ^vðkÞ Þg has the same limit (u*, v*) as the sequence {(u(k), v(k))}. By Lemma 1, we complete the proof. h Remark 5. From the proof of Theorem 4, we can see that the modified simple iteration (6) is more efficient than the simple iteration (5) since uðkÞ < ^ uðkÞ and vðkÞ < ^vðkÞ hold for all k P 3.

1502

L. Bao et al. / Applied Mathematics and Computation 181 (2006) 1499–1504

3. Numerical experiments In this section, we present some numerical experiments to illustrate the performance of the modified simple iterative method (6). Let SI and MSI denote the simple iterative method (5) and the modified simple iterative method (6), respectively. In the following examples, we compare MSI with SI. All numerical experiments are performed on an AMD 1.4 GHz PC with main memory 512 MB. The stopping criterion for two methods is kRk k1 6 1:0e  12; kR0 k1 where Rk denotes the kth residual matrix given by XkCXk  XkE  AXk + B and R0 = B = eeT.

Table 3.1 Example 1, MSI vs. SI (n = 64) c

a 8

6

Method

Iterations

CPU

Residual

10

1  10

SI MSI

17251 13572

27.94 20.38

9.9960e13 9.9969e13

0.001

0.999

SI MSI

746 588

1.08 0.80

9.8010e13 9.9428e13

0.01

0.99

SI MSI

246 195

0.34 0.27

8.9998e13 9.4635e13

0.5

0.5

SI MSI

19 17

0.03 0.02

8.3363e13 6.2768e13

0.85

0.1

SI MSI

8 7

0.01 0.01

2.3443e14 3.4609e13

0 SI MSI _2

log(||Rk||/||R0||)

_4

_6

_8

_10

_12

_14

0

2000

4000

6000

8000

10000

12000

14000

Number of iterations

Fig. 1. Comparison for Example 1 (a = 108, c = 1  106).

16000

18000

L. Bao et al. / Applied Mathematics and Computation 181 (2006) 1499–1504

1503

0 SI MSI _2

log(||Rk||/||R0||)

_4

_6

_8

_10

_12

_14

0

2000

4000

6000

8000

10000

12000

14000

16000

18000

Number of iterations

Fig. 2. Comparison for Example 2 (a = 108, c = 1  106).

Table 3.2 Example 2, MSI vs. SI (n = 128) a

c

Method

Iterations

CPU

Residual

108

1  106

SI MSI

17252 13574

149.34 115.65

9.9711e13 9.9985e13

0.001

0.999

SI MSI

746 588

6.42 4.97

9.8318e13 9.9843e13

0.01

0.99

SI MSI

246 195

2.11 1.68

9.0427e13 9.4950e13

0.5

0.5

SI MSI

19 17

0.16 0.13

8.3379e13 6.2264e13

0.85

0.1

SI MSI

8 7

0.07 0.06

2.3471e14 3.4595e13

Example 1. We consider (1) for n = 64. As in Example 5.2 in [1], the constants ci and xi are given by a numerical quadrature formula on the interval [0, 1], which is obtained by dividing [0, 1] into n/4 subinterval of equal length and applying a Gauss–Legendre quadrature with four nodes to each subinterval. We test five values of (a, c) taken to be (108, 1  106), (0.001, 0.999), (0.01, 0.99), (0.5, 0.5) and (0.85, 0.1). In Table 3.1, we give the results obtained by the two methods. We list the number of iterations, the CPUtime(seconds) required for convergence and the relative residual kRkk1/kR0k1. We can see that the number of iterations for MSI is about 0.75 time of that for SI in case (a, c) = (108, 1  106). See Fig. 1 for a graph of the convergence. Example 2. We now do the same experiment as in the previous example except that n is 128 instead of 64. See Fig. 2 for the computational results. From Table 3.2, we can see that the number of iterations is almost the same for n = 128 as that for n = 64.

1504

L. Bao et al. / Applied Mathematics and Computation 181 (2006) 1499–1504

4. Conclusion In this paper we propose a modified simple iterative method for the minimal positive solution of nonsymmetric Riccati equations arising in transport theory. Numerical experiments show the effectiveness of the proposed method. References [1] C.H. Guo, A.J. Laub, On the iterative solution of a class of nonsymmetric algebraic Riccati equations, SIAM J. Matrix Anal. Appl. 22 (2000) 376–391. [2] J. Juang, Existence of algebraic matrix Riccati equations arising in transport theory, Linear Algebra Appl. 230 (1995) 89–100. [3] J. Juang, W.W. Lin, Nonsymmetric algebraic Riccati equations and Hamiltonian-like matrices, SIAM J. Matrix Anal. Appl. 20 (1999) 228–243. [4] L.Z. Lu, Solution form and simple iteration of a nonsymmetric algebraic Riccati equation arising in transport theory, SIAM J. Matrix Anal. Appl. 26 (2005) 679–685.