# Two-Step Relaxation Newton Method for  Nonsymmetric Algebraic Riccati Equations   Arising from Transport Theory

Authors: Shulin Wu; Chengming Huang
Journal: Mathematical Problems in Engineering
 (2009)

---

## Abstract

We propose a new idea to construct an effective algorithm to compute the minimal positive solution of the nonsymmetric algebraic Riccati equations arising from transport theory. For a class of these equations, an important feature is that the minimal positive solution can be obtained by computing the minimal positive solution of a couple of fixed-point equations with vector form. Based on the fixed-point vector equations, we introduce a new algorithm, namely,two-step relaxation Newton, derived by combining two different relaxation Newton methods to compute the minimal positive solution. The monotone convergence of the solution sequence generated by this new algorithm is established. Numerical results are given to show the advantages of the new algorithm for the nonsymmetric algebraic Riccati equations in vector form.


---

## Body

## 1. Introduction

In this paper, we are interested in iteratively solving the algebraic Riccati equation arising from transport theory (see, e.g., [1–5] and references therein):(1.1)XCX-XE-AX+B=0,
where A,B,C,E∈ℝn×n are given by(1.2)A=Δ-eqT,B=eeT,C=qqT,E=D-qeT.
Here e=(1,1,…,1)T and q=(q1,q2,…,qn)T with qi=ci/2ωi and(1.3)Δ=diag(δ1,δ2,…,δn),δi=1cωi(1+α),D=diag(d1,d2,…,dn),di=1cωi(1-α).
The parameters c and α satisfy 0<c≤1, 0≤α<1, and {ci}i=1n and {ω}i=1n are sets of the composite Gauss-Legendre weights and nodes, respectively, on the interval [0,1] satisfying(1.4)0<ωn<ωn-1<⋯<ω1,ci>0(i=1,2,…,n),∑i=1nci=1;
see, for example, [4] for details. Clearly, it holds that(1.5)0<δ1<δ2<⋯<δn,0<d1<d2<⋯<dn,di=δiforα=0,di>δiforα≠0,i=1,2,…,n.It has been shown that problem (1.1) has positive solution in the sense of component-wise; see [2, 4] for details. Since only the minimal positive solution is physically meaningful, the research in the field of iteration methods centers around the computation of the minimal positive solution of (1.1); see, for example, [3, 5–12]. For the discussion on more general nonsymmetric algebraic Riccati equations arising in real world, we refer the interested reader to [13–17].In the seminal paper by Lu [9], it was shown that the solution X of (1.1) can be written as(1.6)X=T∘(uvT)=(uvT)∘T,
where ∘ denotes the Hadamard product, T is the matrix with elements Ti,j=1/(δi+dj), and u,v are two vectors satisfying the following vector equations:(1.7)f(u,v)=(u-[u∘(Pv)+e]v-[v∘(Qu)+e])=0,
where P=[Pi,j]=qj/(δi+dj) and Q=[Qi,j]=qj/(δj+di).Letw=(uT,vT)T∈ℝ2n with wi=ui and wn+i=vi, i=1,2,…,n, and let g(w)=(g1(w),g2(w),…,g2n(w))T with(1.8)gi(w)={∑j=1nPi,jwn+jfor1≤i≤n,∑j=1nQi-n,jwjforn+1≤i≤2n.
Then the vector equations (1.7) can be uniformly rewritten as(1.9)w=w∘g(w)+e.Based on the vector equations (1.7), Lu [9] investigated the simple iteration (SI) method to compute the minimal positive solution as follows:(1.10)uk+1=uk∘(Pvk)+e,vk+1=vk∘(Quk)+e,v0=0,u0=0.
It was shown that the solution sequence {uk,vk} generated by (1.10) converges monotonically to the minimal positive solution of (1.7). We note that at each step the SI method costs about 4n2 flops (for the definition of flops, see, e.g., [18]), while the Gauss-Jacobi (GJ) scheme defined by Juang [2](1.11)Xk+1=T∘(Xkq+e)(XkTq+e)T
costs about 6n2 flops. Hence, the SI method is more efficient than the GJ method.In (1.10), we note that before computing vk+1, we have obtained uk+1, which should be a better approximation to u than uk. Based on this consideration, Bao et al. [7] constructed the modified simple iteration method as(1.12)uk+1=uk∘(Pvk)+e,vk+1=vk∘(Quk+1)+e,v0=0,u0=0.
It was shown theoretically and numerically that the modified simple iteration procedure (1.12) is more efficient than the original one (1.10).Recently, the so-called nonlinear splitting iteration methods defined as(1.13a)uk+1=uk+1∘[Pvk]+e,vk+1=vk+1∘[Quk]+e,(1.13b)uk+1=uk+1∘[Pvk]+e,vk+1=vk+1∘[Quk+1]+e
were investigated by Bai et al. [6] and independently by Lin [19]. In [6], methods (1.13a) and (1.13b) were called nonlinear block Jacobi (NBJ) iteration method and nonlinear block Gauss-Seidel(NBGS) iteration method, respectively, and it was shown that the solution sequence {uk,vk} generated by the NBJ method and the NBGS method converges monotonically to the minimal positive solution of (1.7). Moreover, numerical experiments given in [6, 19] show that the convergence speed of these two methods is higher than that of the SI method (1.10).In this paper, based on the fixed-point equations (1.7), we present an accelerated version of the NBJ scheme, namely, two-step relaxation Newton method, to compute the minimal positive solution. The construction of the new method is based on the relaxation Newton method introduced by Wu et al. [20, 21]. At each iteration, the new algorithm is composed of two-steps: we use the NBJ method as a simple relaxation Newton method to obtain a coarse approximation of the solution of (1.7), and then with this coarse approximation at hand we use a relative complicated relaxation Newton method to get a finer approximation. It is shown that the solution sequence generated by this method converges monotonically to the minimal positive solution of (1.7). We also prove that the new method is more efficient than the NBJ method and its two-step version.The remainder of this paper is organized as follows. In Section2, we introduce the relaxation Newton method and the its two-step version. Section 3 is concerned with the monotone convergence of the new method. In Section 4, we test some numerical experiments to compare the performance of the new method with the SI method and the NBJ method in the sense of iteration number and CPU time. In the end of the work of this paper, we have constructed another two-step relaxation Newton algorithm which performs much better than the SI, NBJ methods; unfortunately, at the moment we cannot theoretically prove the convergence of this method to the minimal positive solution of the vector equations (1.7). Therefore, we will just report in Section 4 the numerical results of this method.


## 2. The Two-Step Relaxation Newton Method

In this section, we focus on introducing the basic idea of the relaxation Newton method investigated by [20, 21] and its two-step version for the following general nonlinear equations:(2.1)ℱ(x)=0.


### 2.1. The Relaxation Newton Algorithm

The key idea of the relaxation Newton algorithm is to choose some splitting functionF:ℝn×ℝn→ℝn which is minimally assumed to satisfy a consistency condition(2.2)F(x,x)=ℱ(x)
for any x∈ℝn. Then with an initial guess x0 of the unknown solution x*, we start with the previous approximation xk to compute the next approximation xk+1 by solving the following problem:(2.3)F(xk,xk+1)=0
with some conventional method, such as the classical Newton's method, quasi-Newton methods, and Conjugate-Gradient method. Obviously, the generated sequence {xk}0∞ will upon convergence approach to some value x* which satisfies F(x*,x*)=0, that is, ℱ(x)=0. Wu et. al. [20, 21] used the classical Newton's method to solve (2.3), which explains the name: relaxation Newton; the deduced algorithm written compactly is shown in Algorithm 1.Algorithm 1:The relaxation Newton method.
for k=0,1,2,….with a given initial approximation x̃0 of xk+1;for m=0,1,…,MF2(xk,x̃m)Δx̃m=-F(xk,x̃m),x̃m+1=x̃m+Δx̃m,endxk+1=Δx̃MendIn Algorithm1 and what follows(2.4)F2(x,y)=∂F(x,z)∂z|z=y.
If we set x̃0=xk and M=1 in the relaxation algorithm shown in Algorithm 1, by the consistency condition (2.2), the method can be written as(2.5)F2(xk,xk)Δxk=-ℱ(xk),xk+1=xk+Δxk,k=0,1,….
With a special choice of F, the Jacobi matrix F2(x,x) will be a diagonal or block diagonal matrix and invertible in ℝn, and thus iteration method (2.5) can be processed stably and simultaneously with less storage compared with the classical Newton's method. We will see in the next section how to select the splitting function F such that the Jacobi matrix F2(x,x) is a diagonal or block diagonal matrix.


### 2.2. The Two-Step Relaxation Newton Method

Now, suppose that we have two splitting functionsF(x,y), G(x,y,z) and an initial approximation x0 of the solution of (2.1) at hand. We start with the previous approximation xk to compute a mid-approximation xk+1/2 by solving the equations(2.6a)F(xk,xk+1/2)=0,
and then with xk,xk+1/2 at hand we compute xk+1 by solving the equations(2.6b)G(xk,xk+1/2,xk+1)=0.
The splitting functions F and G are assumed to satisfy the consistency conditions(2.7)F(x,x)=ℱ(x),G(x,x,x)=ℱ(x),∀x∈ℝn.
Similar to the relaxation Newton algorithm described above, we use Newton's method with a single iteration to solve equations (2.6a), (2.6b), and the deduced iteration scheme is(2.8a)F2(xk,xk)(xk+1/2-xk)=-ℱ(xk),(2.8b)G3(xk,xk+1/2,xk+1/2)(xk+1-xk+1/2)=-G(xk,xk+1/2,xk+1/2),
where F2 is defined by (2.4), and the Jacobi matrix G3(x,y,z) of the function G is defined by(2.8c)G3(x,y,z)=∂G(x,y,s)∂s|s=z.
Throughout this paper, the iteration scheme (2.8a), (2.8b), (2.8c) is called two-step relaxation Newton (denoted by “TSRN”) algorithm.Specially, for the vector equations (1.7) we consider in this paper the following splitting functions:(2.9a)F(wk,wk+1/2)=(uk+1/2-uk+1/2∘(Pvk)-evk+1/2-vk+1/2∘(Quk)-e),(2.9b)G(wk,wk+1/2,wk+1)=(uk+1-uk+1∘(Pvk+1/2)-uk+1/2∘[Φ(vk+1-vk+1/2)]-evk+1-vk+1∘(Quk+1/2)-vk+1/2∘[Ψ(uk+1-uk+1/2)]-e),
here and hereafter wk+1/2=(uk+1/2T,vk+1/2T)T and wk=(ukT,vkT)T for all k=0,1,2…; Φ and Ψ are diagonal matrices and their diagonal elements are determined as(2.10)Φi,i={Pi,i,ifiisodd,0,ifiiseven,Ψi,i={Qi,i,ifiiseven,0,ifiisodd.
Generally speaking, the splitting function G is a three variables function, but here we consider a special case—only the second and third variables are involved; see (2.9b). In the end of Section 4, we will see another TSRN algorithm where the function G is defined with 3 arguments (see (4.5b)).Define(2.11)Dk=diag(e-Pvk),Λk=diag(e-Quk),D̃k=diag(e-Pvk+1/2),Λ̃k=diag(e-Quk+1/2),Ak=diag(uk+1/2∘[Φe]),Bk=diag(vk+1/2∘[Ψe]).
Then it is clear that(2.12)F2(wk,wk)=[DkΛk],G3(wk,wk+1/2,wk+1/2)=[D̃k-Ak-BkΛ̃k].
Therefore, it follows by substituting the function f defined in (1.7) and the splitting functions F, G defined in (2.9a), (2.9b) into (2.8a), (2.8b), (2.8c) that(2.13a)F2(wk,wk)(uk+1/2vk+1/2)=(ee),(2.13b)G3(wk,wk+1/2,wk+1/2)(uk+1-uk+1/2vk+1-vk+1/2)=-f(uk+1/2,vk+1/2).
We note that with the special splitting functions F and G defined in (2.9a), (2.9b), equations  (2.13a), (2.13b) implies(2.14)F(wk,wk+1/2)=0,G(wk,wk+1/2,wk+1)=0,k=0,1,2,…,
while for general splitting functions this is not always true. For diagonal matrices D̃k, Λ̃k, Ak and Bk, routine calculation yields(2.15)[G3(wk,wk+1/2,wk+1/2)]-1=[[D̃kΛ̃k-AkBk]-1[D̃kΛ̃k-AkBk]-1][Λ̃kAkBkD̃k].
Furthermore, by (2.10) we know that AkBk must be zero matrix for all k=0,1,2…; hence(2.16)[G3(wk,wk+1/2,wk+1/2)]-1=[D̃k-1[D̃kΛ̃k]-1Ak[D̃kΛ̃k]-1BkΛ̃k-1].From (2.9a), (2.9b), (2.10), and (2.16) we obtain the algorithm for implementing the TSRN method as follows.Algorithm 1 (Two-step relaxation Newton iteration).
Starting with the initial value(u0,v0), the solution {uk,vk} is defined, for k=0,1,2,…, by
(1) computing explicitly in the following elementwise fashion:
(2.17)uk+1/2,i=11-[Pvk]i,vk+1/2,i=11-[Quk]i,i=1,2,…,n;
(2) performing fori=1 to n:
ifi is odd
(2.18)uk+1,i=11-[Pvk+1/2]i(1-Φiiuk+1/2,ivk+1/2,i+Φiiuk+1/2,i1-[Quk+1/2]i),vk+1,i=11-[Quk+1/2]i,
else
(2.19)uk+1,i=11-[Pvk+1/2]i,vk+1,i=11-[Quk+1/2]i(1-Ψiiuk+1/2,ivk+1/2,i+Ψiivk+1/2,i1-[Pvk+1/2]i);;
(3) if∥f(uk+1,vk+1)∥∞≤Tol, stop iteration; else set uk=uk+1,vk=vk+1 and go to the first step in Algorithm 1.Clearly, compared with the Newton-type methods investigated by Lu [10] and Lin et al. [8], no LU-factorization is needed at each iteration of the TSRN method.


## 2.1. The Relaxation Newton Algorithm

The key idea of the relaxation Newton algorithm is to choose some splitting functionF:ℝn×ℝn→ℝn which is minimally assumed to satisfy a consistency condition(2.2)F(x,x)=ℱ(x)
for any x∈ℝn. Then with an initial guess x0 of the unknown solution x*, we start with the previous approximation xk to compute the next approximation xk+1 by solving the following problem:(2.3)F(xk,xk+1)=0
with some conventional method, such as the classical Newton's method, quasi-Newton methods, and Conjugate-Gradient method. Obviously, the generated sequence {xk}0∞ will upon convergence approach to some value x* which satisfies F(x*,x*)=0, that is, ℱ(x)=0. Wu et. al. [20, 21] used the classical Newton's method to solve (2.3), which explains the name: relaxation Newton; the deduced algorithm written compactly is shown in Algorithm 1.Algorithm 1:The relaxation Newton method.
for k=0,1,2,….with a given initial approximation x̃0 of xk+1;for m=0,1,…,MF2(xk,x̃m)Δx̃m=-F(xk,x̃m),x̃m+1=x̃m+Δx̃m,endxk+1=Δx̃MendIn Algorithm1 and what follows(2.4)F2(x,y)=∂F(x,z)∂z|z=y.
If we set x̃0=xk and M=1 in the relaxation algorithm shown in Algorithm 1, by the consistency condition (2.2), the method can be written as(2.5)F2(xk,xk)Δxk=-ℱ(xk),xk+1=xk+Δxk,k=0,1,….
With a special choice of F, the Jacobi matrix F2(x,x) will be a diagonal or block diagonal matrix and invertible in ℝn, and thus iteration method (2.5) can be processed stably and simultaneously with less storage compared with the classical Newton's method. We will see in the next section how to select the splitting function F such that the Jacobi matrix F2(x,x) is a diagonal or block diagonal matrix.


## 2.2. The Two-Step Relaxation Newton Method

Now, suppose that we have two splitting functionsF(x,y), G(x,y,z) and an initial approximation x0 of the solution of (2.1) at hand. We start with the previous approximation xk to compute a mid-approximation xk+1/2 by solving the equations(2.6a)F(xk,xk+1/2)=0,
and then with xk,xk+1/2 at hand we compute xk+1 by solving the equations(2.6b)G(xk,xk+1/2,xk+1)=0.
The splitting functions F and G are assumed to satisfy the consistency conditions(2.7)F(x,x)=ℱ(x),G(x,x,x)=ℱ(x),∀x∈ℝn.
Similar to the relaxation Newton algorithm described above, we use Newton's method with a single iteration to solve equations (2.6a), (2.6b), and the deduced iteration scheme is(2.8a)F2(xk,xk)(xk+1/2-xk)=-ℱ(xk),(2.8b)G3(xk,xk+1/2,xk+1/2)(xk+1-xk+1/2)=-G(xk,xk+1/2,xk+1/2),
where F2 is defined by (2.4), and the Jacobi matrix G3(x,y,z) of the function G is defined by(2.8c)G3(x,y,z)=∂G(x,y,s)∂s|s=z.
Throughout this paper, the iteration scheme (2.8a), (2.8b), (2.8c) is called two-step relaxation Newton (denoted by “TSRN”) algorithm.Specially, for the vector equations (1.7) we consider in this paper the following splitting functions:(2.9a)F(wk,wk+1/2)=(uk+1/2-uk+1/2∘(Pvk)-evk+1/2-vk+1/2∘(Quk)-e),(2.9b)G(wk,wk+1/2,wk+1)=(uk+1-uk+1∘(Pvk+1/2)-uk+1/2∘[Φ(vk+1-vk+1/2)]-evk+1-vk+1∘(Quk+1/2)-vk+1/2∘[Ψ(uk+1-uk+1/2)]-e),
here and hereafter wk+1/2=(uk+1/2T,vk+1/2T)T and wk=(ukT,vkT)T for all k=0,1,2…; Φ and Ψ are diagonal matrices and their diagonal elements are determined as(2.10)Φi,i={Pi,i,ifiisodd,0,ifiiseven,Ψi,i={Qi,i,ifiiseven,0,ifiisodd.
Generally speaking, the splitting function G is a three variables function, but here we consider a special case—only the second and third variables are involved; see (2.9b). In the end of Section 4, we will see another TSRN algorithm where the function G is defined with 3 arguments (see (4.5b)).Define(2.11)Dk=diag(e-Pvk),Λk=diag(e-Quk),D̃k=diag(e-Pvk+1/2),Λ̃k=diag(e-Quk+1/2),Ak=diag(uk+1/2∘[Φe]),Bk=diag(vk+1/2∘[Ψe]).
Then it is clear that(2.12)F2(wk,wk)=[DkΛk],G3(wk,wk+1/2,wk+1/2)=[D̃k-Ak-BkΛ̃k].
Therefore, it follows by substituting the function f defined in (1.7) and the splitting functions F, G defined in (2.9a), (2.9b) into (2.8a), (2.8b), (2.8c) that(2.13a)F2(wk,wk)(uk+1/2vk+1/2)=(ee),(2.13b)G3(wk,wk+1/2,wk+1/2)(uk+1-uk+1/2vk+1-vk+1/2)=-f(uk+1/2,vk+1/2).
We note that with the special splitting functions F and G defined in (2.9a), (2.9b), equations  (2.13a), (2.13b) implies(2.14)F(wk,wk+1/2)=0,G(wk,wk+1/2,wk+1)=0,k=0,1,2,…,
while for general splitting functions this is not always true. For diagonal matrices D̃k, Λ̃k, Ak and Bk, routine calculation yields(2.15)[G3(wk,wk+1/2,wk+1/2)]-1=[[D̃kΛ̃k-AkBk]-1[D̃kΛ̃k-AkBk]-1][Λ̃kAkBkD̃k].
Furthermore, by (2.10) we know that AkBk must be zero matrix for all k=0,1,2…; hence(2.16)[G3(wk,wk+1/2,wk+1/2)]-1=[D̃k-1[D̃kΛ̃k]-1Ak[D̃kΛ̃k]-1BkΛ̃k-1].From (2.9a), (2.9b), (2.10), and (2.16) we obtain the algorithm for implementing the TSRN method as follows.Algorithm 1 (Two-step relaxation Newton iteration).
Starting with the initial value(u0,v0), the solution {uk,vk} is defined, for k=0,1,2,…, by
(1) computing explicitly in the following elementwise fashion:
(2.17)uk+1/2,i=11-[Pvk]i,vk+1/2,i=11-[Quk]i,i=1,2,…,n;
(2) performing fori=1 to n:
ifi is odd
(2.18)uk+1,i=11-[Pvk+1/2]i(1-Φiiuk+1/2,ivk+1/2,i+Φiiuk+1/2,i1-[Quk+1/2]i),vk+1,i=11-[Quk+1/2]i,
else
(2.19)uk+1,i=11-[Pvk+1/2]i,vk+1,i=11-[Quk+1/2]i(1-Ψiiuk+1/2,ivk+1/2,i+Ψiivk+1/2,i1-[Pvk+1/2]i);;
(3) if∥f(uk+1,vk+1)∥∞≤Tol, stop iteration; else set uk=uk+1,vk=vk+1 and go to the first step in Algorithm 1.Clearly, compared with the Newton-type methods investigated by Lu [10] and Lin et al. [8], no LU-factorization is needed at each iteration of the TSRN method.


## 3. Convergence Analysis

For our proof of the monotone convergence of the TSRN method, we need the following results proved by Lu [9].Lemma 3.1 (see [9]).
With the functiong defined in (1.8), one has(1)
u*>e,v*>e,(2)
max1≤i≤2n{gi(e)}<1.The first conclusion can be obtained straightforwardly from (1.7).Theorem 3.2.
Let(u0,v0)=(0,0) be the initial point of the TSRN method. The solution sequence {uk,vk}k=1∞ generated by the TSRN method is strictly and monotonically increasing and converges to the minimal positive solution (u*,v*) of the vector equations (1.7); particularly, it holds that  (1)
0≤uk<uk+1/2<uk+1<u*  and  0≤vk<vk+1/2<vk+1<v*, k=0,1,2,…;(2)
limk→+∞uk=u*  and limk→+∞vk=v*.Proof.
We first prove the first conclusion of Theorem3.2 by induction rule and this is completed in two-steps: first, we prove the correctness of this conclusion for k=0, and then under the induction assumption with index k=j, we prove the case for k=j+1.
(1) It is easy to get(u1/2,v1/2)=(e,e). Hence, from the first conclusion of Lemma 3.1, it holds u1/2<u* and v1/2<v*. By the second conclusion of Lemma 3.1, we get e-Pv1/2>0 and e-Qu1/2>0. Therefore, from (2.16) we know [G3(w0,w1/2,w1/2)]-1>0. Since
(3.1)u1/2∘(e-Pv0)-e=0,v1/2∘(e-Qu0)-e=0,u1/2>u0,v1/2>v0,
we have u1/2∘(e-Pv1/2)-e<0 and v1/2∘(e-Qu1/2)-e<0, that is, f(u1/2,v1/2)<0. Hence, it follows from (2.13b) that w1=w1/2-[G3(w0,w1/2,w1/2)]-1f(u1/2,v1/2)>w1/2. We note that for any w̃,w,w-,w+∈ℝ2n it holds
(3.2)G(w̃,w,w-)=G(w̃,w,w+)+G3(w̃,w,w+)(w--w+).
Hence, G(w0,w1/2,w*)+G3(w0,w1/2,w*)(w1-w*)=G(w0,w1/2,w1)=0, that is,
(3.3)w1-w*=-[G3(w0,w1/2,w*)]-1G(w0,w1/2,w*),
where w*=(u*v*). It is easy to verify G3(w0,w1/2,w*)=G3(w0,w1/2,w1/2) and this implies [G3(w0,w1/2,w*)]-1>0. Consider
(3.4)u*∘(e-Pv*)=e,v*∘(e-Qu*)=e,(3.5)u*∘(e-Pv1/2)-u1/2∘Φ(v*-v1/2)-u*∘(e-Pv*)=u*∘[P(v*-v1/2)]-u1/2∘Φ(v*-v1/2),v*∘(e-Qu1/2)-v1/2∘Ψ(u*-u1/2)-v*∘(e-Qu*)=v*∘[Q(u*-u1/2)]-v1/2∘Ψ(u*-u1/2),
and u1/2<u*,v1/2<v*, and we know that the right hand of (3.5) is positive, and this implies G(w0,w1/2,w*)>0. Therefore, from (3.3) we get w1<w*. We have proved that the first conclusion of Theorem 3.2  is correct for k=0.
(2) Assume that the first conclusion of Theorem3.2 is correct for k=j, j≥0 and we will prove that it is also correct for k=j+1. To this end, we note that under the induction assumption it holds
(3.6)e>e-Quj>e-Quj+1/2>e-Quj+1>e-Qu*>0,e>e-Pvj>e-Pvj+1/2>e-Pvj+1>e-Pv*>0,
where the inequalities e-Qu*>0 and e-Pv*>0 follow directly from equality (3.4).
Consider
(3.7)uj+1∘(e-Pvj+1/2)-uj+1/2∘[Φ(vj+1-vj+1/2)]-uj+1∘(e-Pvj+1)=uj+1∘[P(vj+1-vj+1/2)]-uj+1/2∘[Φ(vj+1-vj+1/2)]>0,vj+1∘(e-Quj+1/2)-vj+1/2∘[Ψ(uj+1-uj+1/2)]-vj+1∘(e-Quj+1)=vj+1∘[Q(uj+1-uj+1/2)]-vj+1/2∘[Ψ(uj+1-uj+1/2)]>0,
and since wj+1/2<wj+1, we get
(3.8)uj+1+1/2∘(e-Pvj+1)=e=uj+1∘(e-Pvj+1/2)-uj+1/2∘[Φ(vj+1-vj+1/2)]>uj+1∘(e-Pvj+1),vj+1+1/2∘(e-Quj+1)=e=vj+1∘(e-Quj+1/2)-vj+1/2∘[Ψ(uj+1-uj+1/2)]>uj+1∘(e-Quj+1),
and this coupled with (3.6) gives
(3.9)wj+1<wj+1+1/2.
From (3.6) and wj+1<w*, we have
(3.10)uj+1+1/2∘(e-Pvj+1)=e=u*∘(e-Pv*)<u*∘(e-Pvj+1)vj+1+1/2∘(e-Quj+1)=e=v*∘(e-Qu*)<v*∘(e-Quj+1)⇒uj+1+1/2<u*,vj+1+1/2<v*,
that is, wj+1+1/2<w*; this coupled with (3.6) implies e-Pvj+1+1/2>0 and e-Quj+1+1/2>0. Hence, we arrive at
(3.11)[G3(wj+1,wj+1+1/2,wj+1+1/2)]-1>0.
From (3.9) we have
(3.12)uj+1+1/2∘(e-Pvj+1+1/2)-e<uj+1+1/2∘(e-Pvj+1)-e=0vj+1+1/2∘(e-Quj+1+1/2)-e<vj+1+1/2∘(e-Quj+1)-e=0⇒f(uj+1+1/2,vj+1+1/2)<0;
this coupled with (3.11) gives
(3.13)wj+2=wj+1+1/2-[G3(wj+1,wj+1+1/2,wj+1+1/2)]-1f(uj+1+1/2,vj+1+1/2)>wj+1+1/2.
Next, we provewj+2<w*. To this end, from (3.2) and (2.14) we have
(3.14)G(wj+1,wj+1+1/2,w*)+G3(wj+1,wj+1+1/2,w*)(wj+2-w*)=G(wj+1,wj+1+1/2,wj+2)=0,
and this implies
(3.15)wj+2-w*=-[G3(wj+1,wj+1+1/2,w*)]-1G(wj+1,wj+1+1/2,w*).
From (3.11), we get
(3.16)[G3(wj+1,wj+1+1/2,w*)]-1=[G3(wj+1,wj+1+1/2,wj+1+1/2)]-1>0,
where the first equality can be easily verified from (2.9b). Consider (3.4) and the following relations:
(3.17)u*∘(e-Pvj+1+1/2)-uj+1+1/2∘Φ(v*-vj+1+1/2)-u*∘(e-Pv*)=u*∘[P(v*-vj+1+1/2)]-uj+1+1/2∘Φ(v*-vj+1+1/2)>0,v*∘(e-Quj+1+1/2)-vj+1+1/2∘Ψ(u*-uj+1+1/2)-v*∘(e-Qu*)=v*∘[Q(u*-uj+1+1/2)]-vj+1+1/2∘Ψ(u*-uj+1+1/2)>0,
and since uj+1+1/2<u*,vj+1+1/2<v*, we have G(wj+1,wj+1+1/2,w*)>0; this coupled with (3.15) and (3.16) gives
(3.18)wj+2<w*.
By (3.9), (3.13) and (3.18) we have proved the validity of the first conclusion of Theorem 3.2  for k=j+1. Therefore, by induction rule we have completed the proof of the first conclusion of Theorem 3.2.
From the first conclusion of Theorem3.2, it is obvious that there exist positive vectors û* and v̂* such that
(3.19)limk→+∞uk=limk→+∞uk+1/2=û*,limk→+∞vk=limk→+∞vk+1/2=v̂*.
Therefore, it follows from (2.9a), (2.9b), and (2.14) that
(3.20)û*=û*∘(Pv̂*)+e,v̂*=v̂*∘(Qû*)+e.
This implies that  (û*,v̂*) is a positive solution of the vector equations (1.7). According to the minimal property of (u*,v*), it must hold û*=u* and v̂*=v*.To finish this section, we compare the efficiency of the TSRN method with the NBJ method. To this end, we rewrite the original NBJ iteration scheme (1.13a) into the following form:(3.21)𝒰k+1/2=𝒰k+1/2∘[P𝒱k]+e,𝒱k+1/2=𝒱k+1/2∘[Q𝒰k]+e,𝒰k+1=𝒰k+1∘[P𝒱k+1/2]+e,𝒱k+1=𝒱k+1∘[Q𝒰k+1/2]+e.
Clearly, the TSRN method will reduce into the two-step version NBJ method (3.21) if Ψ=Φ=0 in (2.9b).The following theorem indicates that the TSRN method is more efficient than the NBJ method.Theorem 3.3.
Let both the NBJ method and the TSRN method start with the initial value(0,0), and {uk,vk}, {𝒰k,𝒱k} be the sequences generated by the TSRN method and the NBJ method, respectively. Then it holds that
(3.22)𝒰k+1/2<uk+1/2,𝒱k+1/2<vk+1/2,𝒰k+1<uk+1,𝒱k+1<vk+1
for k≥2.Proof.
It is easy to get(u1/2,v1/2)=(𝒰1/2,𝒱1/2)=(e,e). Hence, we have
(3.23)(u1-𝒰1)∘[e-Pe]=u1/2∘[Φ(v1-v1/2)],(v1-𝒱1)∘[e-Qe]=v1/2∘[Ψ(u1-u1/2)].
This coupled with (2.10), the second conclusion of Lemma 3.1, and the first conclusion of Theorem 3.2 gives
(3.24)u1,i>𝒰1,i,ifiisodd,u1,i=𝒰1,i,ifiiseven,v1,i>𝒱1,i,ifiiseven,v1,i=𝒱1,i,ifiisodd.
Since P,Q>0, it follows from (3.24) that
(3.25)u2+1/2∘[e-Pv1]=e=𝒰2+1/2∘[e-P𝒱1]>𝒰2+1/2∘[e-Pv1],v2+1/2∘[e-Qu1]=e=𝒱2+1/2∘[e-Q𝒰1]>𝒱2+1/2∘[e-Qu1],
where in the last inequality we used the relation (3.6) and the fact established by [6]|the sequence {𝒰k,𝒱k} is strictly and monotonically increasing and converges to the minimal positive solution (u*,v*). From (3.25), we have u2+1/2>𝒰2+1/2 and v2+1/2>𝒱2+1/2.
Consider
(3.26)u3∘[e-Pv2+1/2]-u2+1/2∘[Φ(v3-v2+1/2)]=e=𝒰3∘[e-P𝒱2+1/2]>𝒰3∘[e-Pv2+1/2],v3∘[e-Qu2+1/2]-v2+1/2∘[Ψ(u3-u2+1/2)]=e=𝒱3∘[e-Q𝒰2+1/2]>𝒱3∘[e-Qu2+1/2],
the we arrive at u3>𝒰3 and v3>𝒱3. Therefore, the proof of (3.22) can be completed by using relations (3.25) and (3.26) repeatedly.Remark 3.4.
Since the NBJ method is more feasible than the NBGS method in parallel environment, in this paper we only focus on comparing the efficiency between the TSRN and the NBJ methods.


## 4. Numerical Results

In this section, we compute the minimal positive solution of the vector equations (1.7) by the TSRN method, in order to compare numerically the feasibility and effectiveness of this algorithm with the SI method (1.10) and the NBJ method (1.13a) in the sense of number of iteration (denoted as “IT”) and elapsed CPU time in seconds (denoted as “CPU”). Clearly, these methods are very suitable to be performed in parallel environment. We remark that the Newton type methods coupled with LU-factorization [6, 10], NBGS method (1.13b), and the modified simple iteration (1.12) are sequential methods and not suitable to be implemented in parallel environment.In all experiments, the constantsci and ωi, i=1,2,…,n, are given by the composite Gauss-Legendre quadrature formula on the interval [0,1]. More precisely, the interval [0,1] is first divided into subintervals of equal length, and the composite Gauss-Legendre quadrature formula with 4 nodes is then applied to each subinterval (see [14]). In our implementations, all iterations start with initial value (0,0) and are terminated once the current residual error defined as(4.1)ERRk=max{∥uk-uk∘(Pvk)-e∥∞,∥vk-vk∘(Quk)-e∥∞}
satisfies ERRk≤10-13. All codes are performed in MATLAB (version 7.0) on an Intel (R) Pentium (R) Dual E2110 @ 1.4 GHz PC with memory 1 GB.To make a fair comparison, we rewrite equivalently the SI method (1.10) into two-step fashion as follows:(4.2)uk+1/2=uk∘[Pvk]+e,vk+1/2=vk∘[Quk]+e,uk+1=uk+1/2∘[Pvk+1/2]+e,vk+1=vk+1/2∘[Quk+1/2]+e,
and the vectors uk+1/2,vk+1/2 and uk+1,vk+1 are explicitly computed in elementwise as(4.3)uk+1/2,i=uk,i[Pvk]i+1,vk+1/2,i=vk,i[Quk]i+1,uk+1,i=uk+1/2,i[Pvk+1/2]i+1,vk+1,i=vk+1/2,i[Quk+1/2]i+1,i=1,2,….
The NBJ method is also performed in two-step fashion as defined in (3.21) and the vectors are computed in elementwise as:(4.4)uk+1/2,i=11-[Pvk]i,vk+1/2,i=11-[Quk]i,uk+1,i=11-[Pvk+1/2]i,vk+1,i=11-[Quk+1/2]i,i=1,2,….Example 4.1.
In Table1, for the fixed problem size n=32 but different pairs of (α,c), and in Table 2 for the fixed (α,c)=(10-7,1) but different problem size n, we list ITs and CPUs for the SI, NBJ and TSRN methods, respectively.
We see clearly in Table1 that with fixed problem size, the TSRN method converges faster than the SI and NBJ methods, and as α converges to 0 and c converges to 1, the advantages of the TSRN method are more significant. Moreover, for each pair (α,c), the TSRN method needs less CPU time than the SI and NBJ methods, even through the former needs a little more flops at each iteration. When the parameters α and c are fixed, as the problem size n becomes large, the performance of the TSRN method and the NBJ method becomes close, as shown in Table 2.Table 1
Numerical results forn=32 and different  (α,c).
(α,c)Method(0.1,0.9)(0.001,0.995)(10-5,1-2×10-5)(10-7,1-10-7)(10-12,1-10-12)SI3718123772440571486ITNBJ208410401060931092TSRN208110291049930662SI0.0180.0350.292.77.443CPUNBJ0.0010.0260.2101.5513.525TSRN0.0010.0150.1101.2703.063Table 2
Numerical results for(α,c)=(10-7,1) and different n.
nMethod81632641282565121024SI7031771076714947169771796718897193172014ITNBJ3209031401311233093230724307433074230771TSRN3094530713306643064330637306393067830715SI5.5035.9717.45012.54026.587108.811466.511799.747CPUNBJ2.6872.8653.4416.43711.50147.355203.578751.886TSRN2.3922.5253.2375.60810.56846.458201.795735.711Example 4.2.
In the end of the work of this paper, we have constructed another two-step relaxation Newton algorithm (denoted as TSRN* for the moment) with the splitting functions F and G defined as
(4.5a)F(wk+1/2,wk)=(uk+1/2-uk+1/2∘[Pvk]-evk+1/2-vk+1/2∘[Quk]-e),(4.5b)G(wk,wk+1/2,wk+1)=(uk+1-uk+1∘[Pvk]-(uk+1/2∘[P(vk+1/2-vk)]+Q(uk+1/2-uk))-evk+1-vk+1∘[Quk]-(vk+1/2∘[Q(uk+1/2-uk)]+P(vk+1/2-vk))-e).
With these two splitting functions, we list in Table 3, for fixed problem size n=256 but different (α,c), and in Table 4 for fixed (α,c)=(5×10-9,1-3×10-9) but different problem size, ITs and CPUs for the NBJ, TSRN and TSRN* methods (from the numerical results given in the previous example we believe that the NBJ and TSRN, methods converge faster than the SI method, and thus at the moment we omit the comparison of the TSRN* method with the SI method).
For each problem sizen and (α,c) tested in Tables 3 and 4, it holds that
(4.6)max{∥UTSRN-UTSRN*∥∞,∥VTSRN-VTSRN*∥∞}≤10-14,
where UTSRN,VTSRN and UTSRN*,VTSRN* are the converged vectors generated by the TSRN and TSRN* methods, respectively. Therefore, the TSRN* method really converges to the minimal positive solution of the vector equations (1.7). Moreover, from the numerical results listed in Tables 3 and 4, we see clearly that the TSRN* method performs much better than the TSRN and NBJ methods. Unfortunately, at the moment we cannot prove theoretically the convergence of the TSRN* method to the minimal positive solution of the vector equations (1.7).Table 3
Numerical results forn = 256 and different (a,c).
(α,c)Method(0.01,0.99)(10-6,1-3×10-6)(10-7,1-10-7)(10-8,1-10-9)(10-9,1-10-9)NBJ612417105142946529536ITTSRN602411104892933929310TSRN*46197686122410624114NBJ0.1093.78316.05845.40147.972CPUTSRN0.0983.58115.12143.94844.770TSRN*0.0652.98113.24236.72536.845Table 4
Numerical results for(α,c)=(5×10-9,1-3×10-9) and different n.
nMethod8163264128256384960NBJ2855227834274762732927276272752723427211ITTSRN2740927205271462713027133271432715827196TSRN*2307122659224662236322333223252232222342NBJ2.2782.8273.1075.08112.48142.59595.803480.821CPUTSRN2.1152.4472.9974.75810.92140.98492.919472.829TSRN*1.7781.9712.5033.9588.75033.97076.889441.720



---

*Source: 783920-2009-11-12.xml*