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

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.

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):

𝑋𝐶𝑋−𝑋𝐸−𝐴𝑋+𝐵=0, (1.1) where 𝐴,𝐵,𝐶,𝐸∈ℝ^𝑛×𝑛 are given by

𝐴=Δ−𝑒𝑞^𝑇,𝐵=𝑒𝑒^𝑇,𝐶=𝑞𝑞^𝑇,𝐸=𝐷−𝑞𝑒^𝑇. (1.2) Here 𝑒=(1,1,…,1)^𝑇 and 𝑞=(𝑞₁,𝑞₂,…,𝑞_𝑛)^𝑇 with 𝑞_𝑖=𝑐_𝑖/2𝜔_𝑖 and

Δ=diag(𝛿₁,𝛿₂,…,𝛿_𝑛),𝛿_𝑖=1𝑐𝜔_𝑖,(1+𝛼)𝐷=diag(𝑑₁,𝑑₂,…,𝑑_𝑛),𝑑_𝑖=1𝑐𝜔_𝑖.(1−𝛼) (1.3) The parameters 𝑐 and 𝛼 satisfy 0<𝑐≤1, 0≤𝛼<1, and {𝑐_𝑖}^𝑛_𝑖=1 and {𝜔}^𝑛_𝑖=1 are sets of the composite Gauss-Legendre weights and nodes, respectively, on the interval [0,1] satisfying

0<𝜔_𝑛<𝜔_𝑛−1<⋯<𝜔₁,𝑐_𝑖>0(𝑖=1,2,…,𝑛),^𝑛_𝑖=1𝑐_𝑖=1; (1.4) see, for example, [4] for details. Clearly, it holds that

0<𝛿₁<𝛿₂<⋯<𝛿_𝑛,0<𝑑₁<𝑑₂<⋯<𝑑_𝑛,𝑑_𝑖=𝛿_𝑖for𝛼=0,𝑑_𝑖>𝛿_𝑖for𝛼≠0,𝑖=1,2,…,𝑛. (1.5)

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 𝑋 of (1.1) can be written as

𝑋=𝑇∘(𝑢𝑣^𝑇)=(𝑢𝑣^𝑇)∘𝑇, (1.6) where ∘ denotes the Hadamard product, 𝑇 is the matrix with elements 𝑇_𝑖,𝑗=1/(𝛿_𝑖+𝑑_𝑗), and 𝑢,𝑣 are two vectors satisfying the following vector equations:

𝑓(𝑢,𝑣)=𝑢−[𝑢∘(𝑃𝑣)+𝑒]𝑣−[𝑣∘(𝑄𝑢)+𝑒]=0, (1.7) where 𝑃=[𝑃_𝑖,𝑗]=𝑞_𝑗/(𝛿_𝑖+𝑑_𝑗) and 𝑄=[𝑄_𝑖,𝑗]=𝑞_𝑗/(𝛿_𝑗+𝑑_𝑖).

Let 𝑤=(𝑢^𝑇,𝑣^𝑇)^𝑇∈ℝ^2𝑛 with 𝑤_𝑖=𝑢_𝑖 and 𝑤_𝑛+𝑖=𝑣_𝑖, 𝑖=1,2,…,𝑛, and let 𝑔(𝑤)=(𝑔₁(𝑤),𝑔₂(𝑤),…,𝑔_2𝑛(𝑤))^𝑇 with

𝑔_𝑖(𝑤)=^𝑛𝑗=1𝑃_𝑖,𝑗𝑤_𝑛+𝑗for1≤𝑖≤𝑛,_𝑛𝑗=1𝑄_{𝑖−𝑛,𝑗}𝑤_𝑗for𝑛+1≤𝑖≤2𝑛. (1.8) Then the vector equations (1.7) can be uniformly rewritten as

𝑤=𝑤∘𝑔(𝑤)+𝑒. (1.9)

Based on the vector equations (1.7), Lu [9] investigated the simple iteration (SI) method to compute the minimal positive solution as follows:

𝑢_𝑘+1=𝑢_𝑘∘(𝑃𝑣_𝑘𝑣)+𝑒,_𝑘+1=𝑣_𝑘∘(𝑄𝑢_𝑘𝑣)+𝑒,₀=0,𝑢₀=0. (1.10) It was shown that the solution sequence {𝑢_𝑘,𝑣_𝑘} 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 4𝑛² flops (for the definition of flops, see, e.g., [18]), while the Gauss-Jacobi (GJ) scheme defined by Juang [2]

𝑋_𝑘+1=𝑇∘(𝑋_𝑘𝑞+𝑒)(𝑋^𝑇_𝑘𝑞+𝑒)^𝑇 (1.11) costs about 6𝑛² flops. Hence, the SI method is more efficient than the GJ method.

In (1.10), we note that before computing 𝑣_𝑘+1, we have obtained 𝑢_𝑘+1, which should be a better approximation to 𝑢 than 𝑢_𝑘. Based on this consideration, Bao et al. [7] constructed the modified simple iteration method as

𝑢_𝑘+1=𝑢_𝑘∘(𝑃𝑣_𝑘𝑣)+𝑒,_𝑘+1=𝑣_𝑘∘(𝑄𝑢_𝑘+1𝑣)+𝑒,₀=0,𝑢₀=0. (1.12) 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=𝑢_𝑘+1∘[𝑃𝑣_𝑘𝑣]+𝑒,_𝑘+1=𝑣_𝑘+1∘[𝑄𝑢_𝑘𝑢]+𝑒,_𝑘+1=𝑢_𝑘+1∘[𝑃𝑣_𝑘𝑣]+𝑒,_𝑘+1=𝑣_𝑘+1∘[𝑄𝑢_𝑘+1]+𝑒 (1.13a)(1.13b) 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 {𝑢_𝑘,𝑣_𝑘} 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 Section 2, 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:

ℱ(𝑥)=0. (2.1)

2.1. The Relaxation Newton Algorithm

The key idea of the relaxation Newton algorithm is to choose some splitting function 𝐹:ℝ^𝑛×ℝ^𝑛→ℝ^𝑛 which is minimally assumed to satisfy a consistency condition

𝐹(𝑥,𝑥)=ℱ(𝑥) (2.2) for any 𝑥∈ℝ^𝑛. Then with an initial guess 𝑥₀ of the unknown solution 𝑥^∗, we start with the previous approximation 𝑥_𝑘 to compute the next approximation 𝑥_𝑘+1 by solving the following problem:

𝐹(𝑥_𝑘,𝑥_𝑘+1)=0 (2.3) with some conventional method, such as the classical Newton's method, quasi-Newton methods, and Conjugate-Gradient method. Obviously, the generated sequence {𝑥_𝑘}^∞₀ will upon convergence approach to some value 𝑥^∗ which satisfies 𝐹(𝑥^∗,𝑥^∗)=0, that is, ℱ(𝑥)=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.

for 𝑘=0,1,2,….         with a given initial approximation 𝑥0 of 𝑥𝑘+1;         for 𝑚=0,1,…,𝑀       𝐹2(𝑥𝑘,𝑥𝑚)Δ𝑥𝑚=−𝐹(𝑥𝑘,𝑥𝑚),       𝑥𝑚+1=𝑥𝑚+Δ𝑥𝑚,         end         𝑥𝑘+1=Δ𝑥𝑀   end

Algorithm 1 

The relaxation Newton method.

In Algorithm 1 and what follows

𝐹₂(𝑥,𝑦)=𝜕𝐹(𝑥,𝑧)|𝜕𝑧_𝑧=𝑦. (2.4) If we set 𝑥₀=𝑥_𝑘 and 𝑀=1 in the relaxation algorithm shown in Algorithm 1, by the consistency condition (2.2), the method can be written as

𝐹₂(𝑥_𝑘,𝑥_𝑘)Δ𝑥_𝑘=−ℱ(𝑥_𝑘),𝑥_𝑘+1=𝑥_𝑘+Δ𝑥_𝑘,𝑘=0,1,…. (2.5) With a special choice of 𝐹, the Jacobi matrix 𝐹₂(𝑥,𝑥) will be a diagonal or block diagonal matrix and invertible in ℝ^𝑛, 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 𝐹 such that the Jacobi matrix 𝐹₂(𝑥,𝑥) is a diagonal or block diagonal matrix.

2.2. The Two-Step Relaxation Newton Method

Now, suppose that we have two splitting functions 𝐹(𝑥,𝑦), 𝐺(𝑥,𝑦,𝑧) and an initial approximation 𝑥₀ of the solution of (2.1) at hand. We start with the previous approximation 𝑥_𝑘 to compute a mid-approximation 𝑥_𝑘+1/2 by solving the equations

𝐹(𝑥_𝑘,𝑥_𝑘+1/2)=0, (2.6a) and then with 𝑥_𝑘,𝑥_𝑘+1/2 at hand we compute 𝑥_𝑘+1 by solving the equations

𝐺(𝑥_𝑘,𝑥_𝑘+1/2,𝑥_𝑘+1)=0. (2.6b) The splitting functions 𝐹 and 𝐺 are assumed to satisfy the consistency conditions

𝐹(𝑥,𝑥)=ℱ(𝑥),𝐺(𝑥,𝑥,𝑥)=ℱ(𝑥),∀𝑥∈ℝ^𝑛. (2.7) 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

𝐹₂(𝑥_𝑘,𝑥_𝑘)(𝑥_𝑘+1/2−𝑥_𝑘)=−ℱ(𝑥_𝑘𝐺),₃(𝑥_𝑘,𝑥_𝑘+1/2,𝑥_𝑘+1/2)(𝑥_𝑘+1−𝑥_𝑘+1/2)=−𝐺(𝑥_𝑘,𝑥_𝑘+1/2,𝑥_𝑘+1/2), (2.8a)(2.8b) where 𝐹₂ is defined by (2.4), and the Jacobi matrix 𝐺₃(𝑥,𝑦,𝑧) of the function 𝐺 is defined by

𝐺₃(𝑥,𝑦,𝑧)=𝜕𝐺(𝑥,𝑦,𝑠)|𝜕𝑠_𝑠=𝑧. (2.8c) 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:

𝐹(𝑤_𝑘,𝑤_𝑘+1/2𝑢)=^𝑘+1/2−𝑢_𝑘+1/2∘(𝑃𝑣_𝑘𝑣)−𝑒_𝑘+1/2−𝑣_𝑘+1/2∘(𝑄𝑢_𝑘)−𝑒,𝐺(𝑤_𝑘,𝑤_𝑘+1/2,𝑤_𝑘+1𝑢)=^𝑘+1−𝑢_𝑘+1∘(𝑃𝑣_𝑘+1/2)−𝑢_𝑘+1/2∘[Φ(𝑣_𝑘+1−𝑣_𝑘+1/2𝑣)]−𝑒_𝑘+1−𝑣_𝑘+1∘(𝑄𝑢_𝑘+1/2)−𝑣_𝑘+1/2∘[Ψ(𝑢_𝑘+1−𝑢_𝑘+1/2)]−𝑒, (2.9a)(2.9b) here and hereafter 𝑤_𝑘+1/2=(𝑢^𝑇_𝑘+1/2,𝑣^𝑇_𝑘+1/2)^𝑇 and 𝑤_𝑘=(𝑢^𝑇_𝑘,𝑣^𝑇_𝑘)^𝑇 for all 𝑘=0,1,2…; Φ and Ψ are diagonal matrices and their diagonal elements are determined as

Φ_𝑖,𝑖=𝑃_𝑖,𝑖Ψ,if𝑖isodd,0,if𝑖iseven,^𝑖,𝑖=𝑄_𝑖,𝑖,if𝑖iseven,0,if𝑖isodd. (2.10) Generally speaking, the splitting function 𝐺 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 𝐺 is defined with 3 arguments (see (4.5b)).

Define

𝐷_𝑘=diag(𝑒−𝑃𝑣_𝑘),Λ_𝑘=diag(𝑒−𝑄𝑢_𝑘),𝐷_𝑘=diag(𝑒−𝑃𝑣_𝑘+1/2Λ),_𝑘=diag(𝑒−𝑄𝑢_𝑘+1/2),𝐴_𝑘=diag(𝑢_𝑘+1/2∘[Φ𝑒]),𝐵_𝑘=diag(𝑣_𝑘+1/2∘[Ψ𝑒]). (2.11) Then it is clear that

𝐹₂(𝑤_𝑘,𝑤_𝑘𝐷)=^𝑘Λ_𝑘,𝐺₃(𝑤_𝑘,𝑤_𝑘+1/2,𝑤_𝑘+1/2𝐷)=^𝑘−𝐴_𝑘−𝐵_𝑘Λ_𝑘. (2.12) Therefore, it follows by substituting the function 𝑓 defined in (1.7) and the splitting functions 𝐹, 𝐺 defined in (2.9a), (2.9b) into (2.8a), (2.8b), (2.8c) that

𝐹₂(𝑤_𝑘,𝑤_𝑘)𝑢_𝑘+1/2𝑣_𝑘+1/2=𝑒𝑒𝐺₃,(𝑤_𝑘,𝑤_𝑘+1/2,𝑤_𝑘+1/2)𝑢_𝑘+1−𝑢_𝑘+1/2𝑣_𝑘+1−𝑣_𝑘+1/2=−𝑓(𝑢_𝑘+1/2,𝑣_𝑘+1/2). (2.13a)(2.13b) We note that with the special splitting functions 𝐹 and 𝐺 defined in (2.9a), (2.9b), equations (2.13a), (2.13b) implies

𝐹(𝑤_𝑘,𝑤_𝑘+1/2)=0,𝐺(𝑤_𝑘,𝑤_𝑘+1/2,𝑤_𝑘+1)=0,𝑘=0,1,2,…, (2.14) while for general splitting functions this is not always true. For diagonal matrices 𝐷_𝑘, Λ_𝑘, 𝐴_𝑘 and 𝐵_𝑘, routine calculation yields

[𝐺₃(𝑤_𝑘,𝑤_𝑘+1/2,𝑤_𝑘+1/2)]⁻¹=[𝐷_𝑘Λ_𝑘−𝐴_𝑘𝐵_𝑘]⁻¹[𝐷_𝑘Λ_𝑘−𝐴_𝑘𝐵_𝑘]⁻¹Λ_𝑘𝐴_𝑘𝐵_𝑘𝐷_𝑘. (2.15) Furthermore, by (2.10) we know that 𝐴_𝑘𝐵_𝑘 must be zero matrix for all 𝑘=0,1,2…; hence

[𝐺₃(𝑤_𝑘,𝑤_𝑘+1/2,𝑤_𝑘+1/2)]⁻¹=𝐷_𝑘⁻¹[𝐷_𝑘Λ_𝑘]⁻¹𝐴_𝑘[𝐷_𝑘Λ_𝑘]⁻¹𝐵_𝑘Λ_𝑘⁻¹. (2.16)

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 (𝑢₀,𝑣₀), the solution {𝑢_𝑘,𝑣_𝑘} is defined, for 𝑘=0,1,2,…, by
(1) computing explicitly in the following elementwise fashion:
𝑢_{𝑘+1/2,𝑖}=11−[𝑃𝑣_𝑘]_𝑖,𝑣_{𝑘+1/2,𝑖}=11−[𝑄𝑢_𝑘]_𝑖,𝑖=1,2,…,𝑛; (2.17) ?
(2) performing for 𝑖=1 to 𝑛:
if 𝑖 is odd
𝑢_𝑘+1,𝑖=11−[𝑃𝑣_𝑘+1/2]_𝑖(1−Φ_𝑖𝑖𝑢_{𝑘+1/2,𝑖}𝑣_{𝑘+1/2,𝑖}+Φ_𝑖𝑖𝑢_{𝑘+1/2,𝑖}1−[𝑄𝑢_𝑘+1/2]_𝑖𝑣),_𝑘+1,𝑖=11−[𝑄𝑢_𝑘+1/2]_𝑖, (2.18) else 𝑢_𝑘+1,𝑖=11−[𝑃𝑣_𝑘+1/2]_𝑖,𝑣_𝑘+1,𝑖=11−[𝑄𝑢_𝑘+1/2]_𝑖(1−Ψ_𝑖𝑖𝑢_{𝑘+1/2,𝑖}𝑣_{𝑘+1/2,𝑖}+Ψ_𝑖𝑖𝑣_{𝑘+1/2,𝑖}1−[𝑃𝑣_𝑘+1/2]_𝑖); (2.19);
(3) if ‖𝑓(𝑢_𝑘+1,𝑣_𝑘+1)‖_∞≤Tol, stop iteration; else set 𝑢_𝑘=𝑢_𝑘+1,𝑣_𝑘=𝑣_𝑘+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 function 𝑔 defined in (1.8), one has(1)𝑢^∗>𝑒,𝑣^∗>𝑒,(2)max_{1≤𝑖≤2𝑛}{𝑔_𝑖(𝑒)}<1.

The first conclusion can be obtained straightforwardly from (1.7).

Theorem 3.2. Let (𝑢₀,𝑣₀)=(0,0) be the initial point of the TSRN method. The solution sequence {𝑢_𝑘,𝑣_𝑘}^∞_𝑘=1 generated by the TSRN method is strictly and monotonically increasing and converges to the minimal positive solution (𝑢^∗,𝑣^∗) of the vector equations (1.7); particularly, it holds that (1)0≤𝑢_𝑘<𝑢_𝑘+1/2<𝑢_𝑘+1<𝑢^∗ and 0≤𝑣_𝑘<𝑣_𝑘+1/2<𝑣_𝑘+1<𝑣^∗, 𝑘=0,1,2,…;(2)lim_𝑘→+∞𝑢_𝑘=𝑢^∗ and lim_𝑘→+∞𝑣_𝑘=𝑣^∗.

Proof. We first prove the first conclusion of Theorem 3.2 by induction rule and this is completed in two-steps: first, we prove the correctness of this conclusion for 𝑘=0, and then under the induction assumption with index 𝑘=𝑗, we prove the case for 𝑘=𝑗+1.
(1) It is easy to get (𝑢_1/2,𝑣_1/2)=(𝑒,𝑒). Hence, from the first conclusion of Lemma 3.1, it holds 𝑢_1/2<𝑢^∗ and 𝑣_1/2<𝑣^∗. By the second conclusion of Lemma 3.1, we get 𝑒−𝑃𝑣_1/2>0 and 𝑒−𝑄𝑢_1/2>0. Therefore, from (2.16) we know [𝐺₃(𝑤₀,𝑤_1/2,𝑤_1/2)]⁻¹>0. Since
𝑢_1/2∘(𝑒−𝑃𝑣₀)−𝑒=0,𝑣_1/2∘(𝑒−𝑄𝑢₀)−𝑒=0,𝑢_1/2>𝑢₀,𝑣_1/2>𝑣₀, (3.1) we have 𝑢_1/2∘(𝑒−𝑃𝑣_1/2)−𝑒<0 and 𝑣_1/2∘(𝑒−𝑄𝑢_1/2)−𝑒<0, that is, 𝑓(𝑢_1/2,𝑣_1/2)<0. Hence, it follows from (2.13b) that 𝑤₁=𝑤_1/2−[𝐺₃(𝑤₀,𝑤_1/2,𝑤_1/2)]⁻¹𝑓(𝑢_1/2,𝑣_1/2)>𝑤_1/2. We note that for any 𝑤,𝑤,𝑤₋,𝑤₊∈ℝ^2𝑛 it holds 𝐺(𝑤,𝑤,𝑤₋)=𝐺(𝑤,𝑤,𝑤₊)+𝐺₃(𝑤,𝑤,𝑤₊)(𝑤₋−𝑤₊). (3.2) Hence, 𝐺(𝑤₀,𝑤_1/2,𝑤^∗)+𝐺₃(𝑤₀,𝑤_1/2,𝑤^∗)(𝑤₁−𝑤^∗)=𝐺(𝑤₀,𝑤_1/2,𝑤₁)=0, that is, 𝑤₁−𝑤^∗=−[𝐺₃(𝑤₀,𝑤_1/2,𝑤^∗)]⁻¹𝐺(𝑤₀,𝑤_1/2,𝑤^∗), (3.3) where 𝑤^∗=(^𝑢∗_𝑣∗). It is easy to verify 𝐺₃(𝑤₀,𝑤_1/2,𝑤^∗)=𝐺₃(𝑤₀,𝑤_1/2,𝑤_1/2) and this implies [𝐺₃(𝑤₀,𝑤_1/2,𝑤^∗)]⁻¹>0. Consider 𝑢^∗∘(𝑒−𝑃𝑣^∗)=𝑒,𝑣^∗∘(𝑒−𝑄𝑢^∗𝑢)=𝑒,_∗∘(𝑒−𝑃𝑣_1/2)−𝑢_1/2∘Φ(𝑣^∗−𝑣_1/2)−𝑢^∗∘(𝑒−𝑃𝑣^∗)=𝑢^∗∘[𝑃(𝑣^∗−𝑣_1/2)]−𝑢_1/2∘Φ(𝑣^∗−𝑣_1/2𝑣),_∗∘(𝑒−𝑄𝑢_1/2)−𝑣_1/2∘Ψ(𝑢^∗−𝑢_1/2)−𝑣^∗∘(𝑒−𝑄𝑢^∗)=𝑣^∗∘[𝑄(𝑢^∗−𝑢_1/2)]−𝑣_1/2∘Ψ(𝑢^∗−𝑢_1/2), (3.4)(3.5) and 𝑢_1/2<𝑢^∗,𝑣_1/2<𝑣^∗, and we know that the right hand of (3.5) is positive, and this implies 𝐺(𝑤₀,𝑤_1/2,𝑤^∗)>0. Therefore, from (3.3) we get 𝑤₁<𝑤^∗. We have proved that the first conclusion of Theorem 3.2 is correct for 𝑘=0.
(2) Assume that the first conclusion of Theorem 3.2 is correct for 𝑘=𝑗, 𝑗≥0 and we will prove that it is also correct for 𝑘=𝑗+1. To this end, we note that under the induction assumption it holds
𝑒>𝑒−𝑄𝑢_𝑗>𝑒−𝑄𝑢_𝑗+1/2>𝑒−𝑄𝑢_𝑗+1>𝑒−𝑄𝑢^∗>0,𝑒>𝑒−𝑃𝑣_𝑗>𝑒−𝑃𝑣_𝑗+1/2>𝑒−𝑃𝑣_𝑗+1>𝑒−𝑃𝑣^∗>0, (3.6) where the inequalities 𝑒−𝑄𝑢^∗>0 and 𝑒−𝑃𝑣^∗>0 follow directly from equality (3.4).
Consider
𝑢_𝑗+1∘(𝑒−𝑃𝑣_𝑗+1/2)−𝑢_𝑗+1/2∘[Φ(𝑣_𝑗+1−𝑣_𝑗+1/2)]−𝑢_𝑗+1∘(𝑒−𝑃𝑣_𝑗+1)=𝑢_𝑗+1∘[𝑃(𝑣_𝑗+1−𝑣_𝑗+1/2)]−𝑢_𝑗+1/2∘[Φ(𝑣_𝑗+1−𝑣_𝑗+1/2𝑣)]>0,_𝑗+1∘(𝑒−𝑄𝑢_𝑗+1/2)−𝑣_𝑗+1/2∘[Ψ(𝑢_𝑗+1−𝑢_𝑗+1/2)]−𝑣_𝑗+1∘(𝑒−𝑄𝑢_𝑗+1)=𝑣_𝑗+1∘[𝑄(𝑢_𝑗+1−𝑢_𝑗+1/2)]−𝑣_𝑗+1/2∘[Ψ(𝑢_𝑗+1−𝑢_𝑗+1/2)]>0, (3.7) and since 𝑤_𝑗+1/2<𝑤_𝑗+1, we get 𝑢_𝑗+1+1/2∘(𝑒−𝑃𝑣_𝑗+1)=𝑒=𝑢_𝑗+1∘(𝑒−𝑃𝑣_𝑗+1/2)−𝑢_𝑗+1/2∘[Φ(𝑣_𝑗+1−𝑣_𝑗+1/2)]>𝑢_𝑗+1∘(𝑒−𝑃𝑣_𝑗+1𝑣),_𝑗+1+1/2∘(𝑒−𝑄𝑢_𝑗+1)=𝑒=𝑣_𝑗+1∘(𝑒−𝑄𝑢_𝑗+1/2)−𝑣_𝑗+1/2∘[Ψ(𝑢_𝑗+1−𝑢_𝑗+1/2)]>𝑢_𝑗+1∘(𝑒−𝑄𝑢_𝑗+1), (3.8) and this coupled with (3.6) gives 𝑤_𝑗+1<𝑤_𝑗+1+1/2. (3.9) From (3.6) and 𝑤_𝑗+1<𝑤^∗, we have 𝑢_𝑗+1+1/2∘(𝑒−𝑃𝑣_𝑗+1)=𝑒=𝑢^∗∘(𝑒−𝑃𝑣^∗)<𝑢^∗∘(𝑒−𝑃𝑣_𝑗+1)𝑣_𝑗+1+1/2∘(𝑒−𝑄𝑢_𝑗+1)=𝑒=𝑣^∗∘(𝑒−𝑄𝑢^∗)<𝑣^∗∘(𝑒−𝑄𝑢_𝑗+1)⇒𝑢_𝑗+1+1/2<𝑢^∗,𝑣_𝑗+1+1/2<𝑣^∗, (3.10) that is, 𝑤_𝑗+1+1/2<𝑤^∗; this coupled with (3.6) implies 𝑒−𝑃𝑣_𝑗+1+1/2>0 and 𝑒−𝑄𝑢_𝑗+1+1/2>0. Hence, we arrive at [𝐺₃(𝑤_𝑗+1,𝑤_𝑗+1+1/2,𝑤_𝑗+1+1/2)]⁻¹>0. (3.11) From (3.9) we have 𝑢_𝑗+1+1/2∘(𝑒−𝑃𝑣_𝑗+1+1/2)−𝑒<𝑢_𝑗+1+1/2∘(𝑒−𝑃𝑣_𝑗+1𝑣)−𝑒=0_𝑗+1+1/2∘(𝑒−𝑄𝑢_𝑗+1+1/2)−𝑒<𝑣_𝑗+1+1/2∘(𝑒−𝑄𝑢_𝑗+1)−𝑒=0⇒𝑓(𝑢_𝑗+1+1/2,𝑣_𝑗+1+1/2)<0; (3.12) this coupled with (3.11) gives 𝑤_𝑗+2=𝑤_𝑗+1+1/2−[𝐺₃(𝑤_𝑗+1,𝑤_𝑗+1+1/2,𝑤_𝑗+1+1/2)]⁻¹𝑓(𝑢_𝑗+1+1/2,𝑣_𝑗+1+1/2)>𝑤_𝑗+1+1/2. (3.13)?
Next, we prove 𝑤_𝑗+2<𝑤^∗. To this end, from (3.2) and (2.14) we have
𝐺(𝑤_𝑗+1,𝑤_𝑗+1+1/2,𝑤^∗)+𝐺₃(𝑤_𝑗+1,𝑤_𝑗+1+1/2,𝑤^∗)(𝑤_𝑗+2−𝑤^∗)=𝐺(𝑤_𝑗+1,𝑤_𝑗+1+1/2,𝑤_𝑗+2)=0, (3.14) and this implies 𝑤_𝑗+2−𝑤^∗=−[𝐺₃(𝑤_𝑗+1,𝑤_𝑗+1+1/2,𝑤^∗)]⁻¹𝐺(𝑤_𝑗+1,𝑤_𝑗+1+1/2,𝑤^∗). (3.15) From (3.11), we get [𝐺₃(𝑤_𝑗+1,𝑤_𝑗+1+1/2,𝑤^∗)]⁻¹=[𝐺₃(𝑤_𝑗+1,𝑤_𝑗+1+1/2,𝑤_𝑗+1+1/2)]⁻¹>0, (3.16) where the first equality can be easily verified from (2.9b). Consider (3.4) and the following relations: 𝑢^∗∘(𝑒−𝑃𝑣_𝑗+1+1/2)−𝑢_𝑗+1+1/2∘Φ(𝑣^∗−𝑣_𝑗+1+1/2)−𝑢^∗∘(𝑒−𝑃𝑣^∗)=𝑢^∗∘[𝑃(𝑣^∗−𝑣_𝑗+1+1/2)]−𝑢_𝑗+1+1/2∘Φ(𝑣^∗−𝑣_𝑗+1+1/2𝑣)>0,_∗∘(𝑒−𝑄𝑢_𝑗+1+1/2)−𝑣_𝑗+1+1/2∘Ψ(𝑢^∗−𝑢_𝑗+1+1/2)−𝑣^∗∘(𝑒−𝑄𝑢^∗)=𝑣^∗∘[𝑄(𝑢^∗−𝑢_𝑗+1+1/2)]−𝑣_𝑗+1+1/2∘Ψ(𝑢^∗−𝑢_𝑗+1+1/2)>0, (3.17) and since 𝑢_𝑗+1+1/2<𝑢^∗,𝑣_𝑗+1+1/2<𝑣^∗, we have 𝐺(𝑤_𝑗+1,𝑤_𝑗+1+1/2,𝑤^∗)>0; this coupled with (3.15) and (3.16) gives 𝑤_𝑗+2<𝑤^∗. (3.18) By (3.9), (3.13) and (3.18) we have proved the validity of the first conclusion of Theorem 3.2 for 𝑘=𝑗+1. Therefore, by induction rule we have completed the proof of the first conclusion of Theorem 3.2.
From the first conclusion of Theorem 3.2, it is obvious that there exist positive vectors 𝑢^∗ and 𝑣^∗ such that
lim_𝑘→+∞𝑢_𝑘=lim_𝑘→+∞𝑢_𝑘+1/2=𝑢^∗,lim_𝑘→+∞𝑣_𝑘=lim_𝑘→+∞𝑣_𝑘+1/2=𝑣^∗. (3.19) Therefore, it follows from (2.9a), (2.9b), and (2.14) that 𝑢^∗=𝑢^∗∘(𝑃𝑣^∗𝑣)+𝑒,_∗=𝑣^∗∘(𝑄𝑢^∗)+𝑒. (3.20) This implies that (𝑢^∗,𝑣^∗) is a positive solution of the vector equations (1.7). According to the minimal property of (𝑢^∗,𝑣^∗), it must hold 𝑢^∗=𝑢^∗ and 𝑣^∗=𝑣^∗.

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:

𝒰_𝑘+1/2=𝒰_𝑘+1/2∘[𝑃𝒱_𝑘𝒱]+𝑒,_𝑘+1/2=𝒱_𝑘+1/2∘[𝑄𝒰_𝑘𝒰]+𝑒,_𝑘+1=𝒰_𝑘+1∘[𝑃𝒱_𝑘+1/2𝒱]+𝑒,_𝑘+1=𝒱_𝑘+1∘[𝑄𝒰_𝑘+1/2]+𝑒. (3.21) 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 {𝑢_𝑘,𝑣_𝑘}, {𝒰_𝑘,𝒱_𝑘} be the sequences generated by the TSRN method and the NBJ method, respectively. Then it holds that 𝒰_𝑘+1/2<𝑢_𝑘+1/2,𝒱_𝑘+1/2<𝑣_𝑘+1/2,𝒰_𝑘+1<𝑢_𝑘+1,𝒱_𝑘+1<𝑣_𝑘+1 (3.22) for 𝑘≥2.

Proof. It is easy to get (𝑢_1/2,𝑣_1/2)=(𝒰_1/2,𝒱_1/2)=(𝑒,𝑒). Hence, we have (𝑢₁−𝒰₁)∘[𝑒−𝑃𝑒]=𝑢_1/2∘[Φ(𝑣₁−𝑣_1/2)],(𝑣₁−𝒱₁)∘[𝑒−𝑄𝑒]=𝑣_1/2∘[Ψ(𝑢₁−𝑢_1/2)]. (3.23) This coupled with (2.10), the second conclusion of Lemma 3.1, and the first conclusion of Theorem 3.2 gives 𝑢_1,𝑖>𝒰_1,𝑖𝑢,if𝑖isodd,_1,𝑖=𝒰_1,𝑖𝑣,if𝑖iseven,_1,𝑖>𝒱_1,𝑖𝑣,if𝑖iseven,_1,𝑖=𝒱_1,𝑖,if𝑖isodd. (3.24) Since 𝑃,𝑄>0, it follows from (3.24) that 𝑢_2+1/2∘[𝑒−𝑃𝑣₁]=𝑒=𝒰_2+1/2∘[𝑒−𝑃𝒱₁]>𝒰_2+1/2∘[𝑒−𝑃𝑣₁𝑣],_2+1/2∘[𝑒−𝑄𝑢₁]=𝑒=𝒱_2+1/2∘[𝑒−𝑄𝒰₁]>𝒱_2+1/2∘[𝑒−𝑄𝑢₁], (3.25) where in the last inequality we used the relation (3.6) and the fact established by [6]|the sequence {𝒰_𝑘,𝒱_𝑘} is strictly and monotonically increasing and converges to the minimal positive solution (𝑢^∗,𝑣^∗). From (3.25), we have 𝑢_2+1/2>𝒰_2+1/2 and 𝑣_2+1/2>𝒱_2+1/2.
Consider
𝑢₃∘[𝑒−𝑃𝑣_2+1/2]−𝑢_2+1/2∘[Φ(𝑣₃−𝑣_2+1/2)]=𝑒=𝒰₃∘[𝑒−𝑃𝒱_2+1/2]>𝒰₃∘[𝑒−𝑃𝑣_2+1/2𝑣],₃∘[𝑒−𝑄𝑢_2+1/2]−𝑣_2+1/2∘[Ψ(𝑢₃−𝑢_2+1/2)]=𝑒=𝒱₃∘[𝑒−𝑄𝒰_2+1/2]>𝒱₃∘[𝑒−𝑄𝑢_2+1/2], (3.26) the we arrive at 𝑢₃>𝒰₃ and 𝑣₃>𝒱₃. 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 constants 𝑐_𝑖 and 𝜔_𝑖, 𝑖=1,2,…,𝑛, 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

ERR_𝑘=max‖𝑢_𝑘−𝑢_𝑘∘(𝑃𝑣_𝑘)−𝑒‖_∞,‖𝑣_𝑘−𝑣_𝑘∘(𝑄𝑢_𝑘)−𝑒‖_∞ (4.1) satisfies ERR_𝑘≤10⁻¹³. 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:

𝑢_𝑘+1/2=𝑢_𝑘∘[𝑃𝑣_𝑘𝑣]+𝑒,_𝑘+1/2=𝑣_𝑘∘[𝑄𝑢_𝑘𝑢]+𝑒,_𝑘+1=𝑢_𝑘+1/2∘[𝑃𝑣_𝑘+1/2𝑣]+𝑒,_𝑘+1=𝑣_𝑘+1/2∘[𝑄𝑢_𝑘+1/2]+𝑒, (4.2) and the vectors 𝑢_𝑘+1/2,𝑣_𝑘+1/2 and 𝑢_𝑘+1,𝑣_𝑘+1 are explicitly computed in elementwise as

𝑢_{𝑘+1/2,𝑖}=𝑢_𝑘,𝑖[𝑃𝑣_𝑘]_𝑖𝑣+1,_{𝑘+1/2,𝑖}=𝑣_𝑘,𝑖[𝑄𝑢_𝑘]_𝑖𝑢+1,_𝑘+1,𝑖=𝑢_{𝑘+1/2,𝑖}[𝑃𝑣_𝑘+1/2]_𝑖𝑣+1,_𝑘+1,𝑖=𝑣_{𝑘+1/2,𝑖}[𝑄𝑢_𝑘+1/2]_𝑖+1,𝑖=1,2,…. (4.3) The NBJ method is also performed in two-step fashion as defined in (3.21) and the vectors are computed in elementwise as:

𝑢_{𝑘+1/2,𝑖}=11−[𝑃𝑣_𝑘]_𝑖,𝑣_{𝑘+1/2,𝑖}=11−[𝑄𝑢_𝑘]_𝑖,𝑢_𝑘+1,𝑖=11−[𝑃𝑣_𝑘+1/2]_𝑖,𝑣_𝑘+1,𝑖=11−[𝑄𝑢_𝑘+1/2]_𝑖,𝑖=1,2,…. (4.4)

Example 4.1. In Table 1, for the fixed problem size 𝑛=32 but different pairs of (𝛼,𝑐), and in Table 2 for the fixed (𝛼,𝑐)=(10⁻⁷,1) but different problem size 𝑛, we list ITs and CPUs for the SI, NBJ and TSRN methods, respectively.
We see clearly in Table 1 that with fixed problem size, the TSRN method converges faster than the SI and NBJ methods, and as 𝛼 converges to 0 and 𝑐 converges to 1, the advantages of the TSRN method are more significant. Moreover, for each pair (𝛼,𝑐), 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 𝑐 are fixed, as the problem size 𝑛 becomes large, the performance of the TSRN method and the NBJ method becomes close, as shown in Table 2.

Table 1 

Numerical results for 𝑛=32 and different (𝛼,𝑐).

Table 2 

Numerical results for (𝛼,𝑐)=(10⁻⁷,1) and different n .

Example 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 𝐹 and 𝐺 defined as 𝐹(𝑤_𝑘+1/2,𝑤_𝑘𝑢)=^𝑘+1/2−𝑢_𝑘+1/2∘[𝑃𝑣_𝑘𝑣]−𝑒_𝑘+1/2−𝑣_𝑘+1/2∘[𝑄𝑢_𝑘]−𝑒,𝐺(𝑤_𝑘,𝑤_𝑘+1/2,𝑤_𝑘+1𝑢)=^𝑘+1−𝑢_𝑘+1∘[𝑃𝑣_𝑘]−(𝑢_𝑘+1/2∘[𝑃(𝑣_𝑘+1/2−𝑣_𝑘)]+𝑄(𝑢_𝑘+1/2−𝑢_𝑘𝑣))−𝑒_𝑘+1−𝑣_𝑘+1∘[𝑄𝑢_𝑘]−(𝑣_𝑘+1/2∘[𝑄(𝑢_𝑘+1/2−𝑢_𝑘)]+𝑃(𝑣_𝑘+1/2−𝑣_𝑘))−𝑒. (4.5a)(4.5b) With these two splitting functions, we list in Table 3, for fixed problem size 𝑛=256 but different (𝛼,𝑐), and in Table 4 for fixed (𝛼,𝑐)=(5×10⁻⁹,1−3×10⁻⁹) 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 size 𝑛 and (𝛼,𝑐) tested in Tables 3 and 4, it holds that max{‖𝑈_TSRN−𝑈_TSRN^∗‖_∞,‖𝑉_TSRN−𝑉_TSRN^∗‖_∞}≤10⁻¹⁴, (4.6) where 𝑈_TSRN,𝑉_TSRN and 𝑈_TSRN^∗,𝑉_TSRN^∗ 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 for n = 256 and different (a,c).

Table 4 

Numerical results for (𝛼,𝑐)=(5×10⁻⁹,1−3×10⁻⁹) and different n.

Acknowledgments

This work was supported by the NSF of China (no. 10971077) and by the program for NCET, the State Education Ministry of China.