Iterative Algorithm for Solving a System of Nonlinear Matrix Equations

Abstract

We discuss the positive definite solutions for the system of nonlinear matrix equations 𝑋−𝐴^∗𝑌^−𝑛𝐴=𝐼 and 𝑌−𝐵^∗𝑋^−𝑚𝐵=𝐼, where 𝑛, 𝑚 are two positive integers. Some properties of solutions are studied, and the necessary and sufficient conditions for the existence of positive definite solutions are given. An iterative algorithm for obtaining positive definite solutions of the system is proposed. Moreover, the error estimations are found. Finally, some numerical examples are given to show the efficiency of the proposed iterative algorithm.

1. Introduction

In this paper, we consider the system of nonlinear matrix equations that can be expressed in the form 𝑋−𝐴^∗𝑌^−𝑛𝐴=𝐼,𝑌−𝐵^∗𝑋^−𝑚𝐵=𝐼, (1.1) where 𝑛,𝑚 are two positive integers, 𝑋,𝑌 are 𝑟×𝑟 unknown matrices, 𝐼 is the 𝑟×𝑟 identity matrix, and 𝐴,𝐵 are given nonsingular matrices. All matrices are defined over the complex field.

System of nonlinear matrix equations of the form of (1.1) is a special case of the system of algebraic discrete-type Riccati equations of the form 𝑋_𝑖=𝑉^∗_𝑖𝑋_𝑖𝑉_𝑖+𝑄_𝑖−𝑉^∗_𝑖𝑋_𝑖𝐵_𝑖+𝐴_𝑖𝑅_𝑖+𝐵^∗_𝑖𝑋_𝑖𝐵_𝑖⁻¹𝐵^∗_𝑖𝑋_𝑖𝑉_𝑖+𝐴^∗_𝑖, (1.2) where 𝑖=1,2,…,𝑘 [1, 2], when 𝑚=𝑛=1. It is well known that the algebraic Riccati equations often arise in control theory, stability theory, communication system, dynamic programming, signal processing, statistics, and so forth [1–3]. In the recent years, some special case of the system (1.2) has been studied in many papers [4–14]. For example, Costa and Marques [5] have studied the maximal and stabilizing hermitian solutions for discrete-time-coupled algebraic Riccati equations. Czornik and Świerniak [7, 8] have studied the lower and the upper bounds on the solution of coupled algebraic Riccati equation. In [13] Mukaidani et al. proposed a numerical algorithm for finding solution of cross-coupled algebraic Riccati equations. In [4] Aldubiban has studied the properties of special case of Sys. (1.2) and obtained the sufficient conditions for existence of a positive definite solution and proposed an iterative algorithm to calculate the solutions. In [10] Ivanov proposed a method to solve the discrete-time-coupled algebraic Riccati equations.

This paper is organized as following: in Section 2, we derive the necessary and sufficient conditions of existence solutions for the Sys. (1.1). In Section 3, we introduce an iterative algorithm to obtain the positive definite solutions of Sys. (1.1). We discuss the convergence of the proposed iterative algorithm and study the convergence of the algorithm. Some numerical examples are given to illustrate the efficiency for suggested algorithm in Section 4.

The following notations are used throughout the rest of the paper. The notation 𝐴≥0(𝐴>0) means that 𝐴 is a positive semidefinite (positive definite), 𝐴^⋆ denotes the complex conjugate transpose of 𝐴, and 𝐼 is the identity matrix. Moreover, 𝐴≥𝐵(𝐴>𝐵) is used as a different notation for 𝐴−𝐵≥0(𝐴−𝐵>0). We denote by 𝜌(𝐴) the spectral radius of 𝐴, 𝜆_𝑟(𝑋), and 𝜇_𝑟(𝑌) mean the eigenvalues of 𝑋 and 𝑌, respectively. The norm used in this paper is the spectral norm of the matrix 𝐴, that is, ‖𝐴‖=𝜌(𝐴𝐴^⋆) unless otherwise noted.

2. Conditions for Existence of the Solutions

In this section, we will discuss some properties of the solutions for the Sys. (1.1), and we obtain the necessary and sufficient conditions for the existence of the solutions of Sys. (1.1).

Theorem 2.1. If 𝜆₋,𝜆₊ are the smallest and the largest eigenvalues of a solution X of Sys. (1.1), respectively, and 𝜇₋,𝜇₊ are the smallest and the largest eigenvalues of a solution Y of Sys. (1.1), respectively, 𝜂,𝜉 are eigenvalues of 𝐴,𝐵, then 𝜆₋−1𝜇^𝑛₋≤𝜂≤𝜆₊−1𝜇^𝑛₊,𝜇₋−1𝜆^𝑚₋≤𝜉≤𝜇₊−1𝜆^𝑚₊. (2.1)

Proof. Let 𝜈 be an eigenvector corresponding to an eigenvalue 𝜂 of the matrix 𝐴 and |𝜈|=1, and let 𝜔 be an eigenvector corresponding to an eigenvalue 𝜉 of the matrix 𝐵 and |𝜔|=1. Since the solution (𝑋,𝑌) of Sys. (1.1) is a positive definite solution then (𝜆₋−1)𝐼≤𝑋−𝐼≤(𝜆₊−1)𝐼 and (𝜇₋−1)𝐼≤𝑌−𝐼≤(𝜇₊−1)𝐼.
From the Sys. (1.1), we have ⟨𝑋𝜈,𝜈⟩−𝐴^∗𝑌^−𝑛𝐴=⟨⟩,⟨⟩−⟨⟩=𝑌𝜈,𝜈𝐼𝜈,𝜈𝑋𝜈,𝜈𝐼𝜈,𝜈^−𝑛,)⟩=𝜂𝐴𝜈,𝐴𝜈⟨(𝑋−𝐼𝜈,𝜈²𝑌^−𝑛𝜈,𝜈, (2.2) that is 𝜆₋−1𝜇^𝑛₋≤𝜂²≤𝜆₊−1𝜇^𝑛₊. (2.3) Hence 𝜆₋−1𝜇^𝑛₋≤𝜂≤𝜆₊−1𝜇^𝑛₊. (2.4) Also, from the Sys. (1.1), we have ⟨𝑌𝜔,𝜔⟩−𝐵^∗𝑋^−𝑚𝐵=⟨⟩,⟨⟩−⟨⟩=𝑋𝜔,𝜔𝐼𝜔,𝜔𝑌𝜔,𝜔𝐼𝜔,𝜔^−𝑚,)⟩=𝜉𝐵𝜔,𝐵𝜔⟨(𝑌−𝐼𝜔,𝜔²𝑋^−𝑚𝜔,𝜔, (2.5) that is 𝜇₋−1𝜆^𝑚₋≤𝜉²≤𝜇₊−1𝜆^𝑚₊. (2.6) Hence 𝜇₋−1𝜆^𝑚₋≤𝜉≤𝜇₊−1𝜆^𝑚₊. (2.7)

Theorem 2.2. If Sys. (1.1) has a positive definite solution (𝑋,𝑌), then 𝐼<𝑋<𝐼+𝐴^∗𝐴,𝐼<𝑌<𝐼+𝐵^∗𝐵. (2.8)

Proof. Since (𝑋,𝑌) is a positive definite solution of Sys. (1.1), then 𝑋>𝐼,𝐴^∗𝑌^−𝑛𝐴>0,𝑌>𝐼,𝐵^∗𝑋^−𝑚𝐵>0. (2.9) From the inequality 𝑌⁻¹<𝐼, we have 𝐴^∗𝑌^−𝑛𝐴<𝐴^∗𝐴, (2.10) that is 𝑋=𝐼+𝐴^∗𝑌^−𝑛𝐴,<𝐼+𝐴^∗𝐴. (2.11) Hence 𝐼<𝑋<𝐼+𝐴^∗𝐴. (2.12) From the inequality 𝑋⁻¹<𝐼, we have 𝐵^∗𝑋^−𝑚𝐵<𝐵^∗𝐵, (2.13) that is 𝑌=𝐼+𝐵^∗𝑋^−𝑚𝐵,<𝐼+𝐵^∗𝐵. (2.14) Hence 𝐼<𝑌<𝐼+𝐵^∗𝐵. (2.15)

Theorem 2.3. Sys. (1.1) has a positive definite solution (𝑋,𝑌) if and only if the matrices 𝐴,𝐵 have the factorization 𝐴=𝑃^∗𝑃^𝑛/2𝑁,𝐵=𝑄^∗𝑄^𝑚/2𝑀, (2.16) where 𝑃,𝑄 are nonsingular matrices satisfying the following system: 𝑄^∗𝑄−𝑁^∗𝑃𝑁=𝐼,_∗𝑃−𝑀^∗𝑀=𝐼. (2.17) In this case the solution is (𝑄^∗𝑄,𝑃^∗𝑃).

Proof. Let Sys. (1.1) have a positive definite solution (𝑋,𝑌); then 𝑋=𝑄^∗𝑄,𝑌=𝑃^∗𝑃, where 𝑄,𝑃 are nonsingular matrices. Then Sys. (1.1) can be rewritten as 𝑄^∗𝑄−𝐴^∗𝑃^∗𝑃^−𝑛𝑃𝐴=𝐼,_∗𝑃−𝐵^∗𝑄^∗𝑄^−𝑚𝑄𝐵=𝐼,_∗𝑄−𝐴^∗𝑃^∗𝑃^−𝑛/2𝑃^∗𝑃^−𝑛/2𝑃𝐴=𝐼,_∗𝑃−𝐵^∗𝑄^∗𝑄^−𝑚/2𝑄^∗𝑄^−𝑚/2𝐵=𝐼. (2.18) Letting 𝑁=(𝑃^∗𝑃)^−𝑛/2𝐴, 𝑀=(𝑄^∗𝑄)^−𝑚/2𝐵, then 𝐴=(𝑃^∗𝑃)^𝑛/2𝑁, 𝐵=(𝑄^∗𝑄)^𝑚/2𝑀, then the Sys. (1.1) is an equivalent to Sys. (2.17).
Conversely, if 𝐴,𝐵 have the factorization (2.16) and satisfy Sys. (2.17), let 𝑋=𝑄^∗𝑄,𝑌=𝑃^∗𝑃, then 𝑋,𝑌 are positive definite matrices, and we have 𝑋−𝐴^∗𝑌^−𝑛𝐴=𝑄^∗𝑄−𝐴^∗𝑃^∗𝑃^−𝑛𝐴=𝑄^∗𝑄−𝐴^∗𝑃^∗𝑃^−𝑛/2𝑃^∗𝑃^−𝑛/2𝐴=𝑄^∗𝑄−𝑁^∗𝑁=𝐼,𝑌−𝐵^∗𝑋^−𝑚𝐵=𝑃^∗𝑃−𝐵^∗𝑄^∗𝑄^−𝑚𝐵=𝑃^∗𝑃−𝐵^∗𝑄^∗𝑄^−𝑚/2𝑄^∗𝑄^−𝑚/2𝐵=𝑃^∗𝑃−𝑀^∗𝑀=𝐼. (2.19)(2.20) Hence Sys. (1.1) has a positive definite solution.

3. Iterative Algorithm for Solving the System

In this section, we will investigate the iterative solution of the Sys. (1.1). From this section to the end of the paper we will consider that the matrices 𝐴,𝐵 are normal satisyfing 𝐴⁻¹𝐵=𝐵𝐴⁻¹ and 𝐴⁻¹𝐵^∗=𝐵^∗𝐴⁻¹.

Let us consider the following iterative algorithm.

Algorithm 3.1. Take 𝑋₀=𝐼,𝑌₀=𝐼.
For𝑠=0,1,2,… compute 𝑋_𝑠+1=𝐼+𝐴^∗𝑌_𝑠^−𝑛𝑌𝐴,_𝑠+1=𝐼+𝐵^∗𝑋_𝑠^−𝑚𝐵. (3.1)

Lemma 3.2. For the Sys. (1.1), we have 𝐴𝑋_𝑠=𝑋_𝑠𝐴,𝐵𝑌_𝑠=𝑌_𝑠𝐵,𝐴𝑌_𝑠⁻¹=𝑌_𝑠⁻¹𝐴,𝐵𝑋_𝑠⁻¹=𝑋_𝑠⁻¹𝐵, (3.2) where {𝑋_𝑠},{𝑌_𝑠},𝑠=0,1,2,…, are determined by Algorithm 3.1.

Proof. Since 𝑋₀=𝑌₀=𝐼, then 𝐴𝑋₀=𝑋₀𝐴,𝐵𝑌₀=𝑌₀𝐵,𝐴⁻¹𝑌₀=𝑌₀𝐴⁻¹,𝐵⁻¹𝑋₀=𝑋₀𝐵⁻¹. (3.3) Using the conditions 𝐴𝐴^∗=𝐴^∗𝐴,𝐵𝐵^∗=𝐵^∗𝐵, we obtain 𝐴𝑋₁=𝐴𝐼+𝐴^∗𝐴=𝐴+𝐴𝐴^∗𝐴=𝐴+𝐴^∗=𝐴𝐴𝐼+𝐴^∗𝐴𝐴=𝑋₁𝐴. (3.4) Also, we have 𝐵𝑌₁=𝑌₁𝐵. (3.5) Using the conditions 𝐴⁻¹𝐵=𝐵𝐴⁻¹,𝐴⁻¹𝐵^∗=𝐵^∗𝐴⁻¹, we obtain 𝐴⁻¹𝑌₁=𝐴⁻¹𝐼+𝐵^∗𝐵=𝐴⁻¹+𝐴⁻¹𝐵^∗𝐵=𝐴⁻¹+𝐵^∗𝐴⁻¹𝐵=𝐴⁻¹+𝐵^∗𝐵𝐴⁻¹=𝐼+𝐵^∗𝐵𝐴⁻¹=𝑌₁𝐴⁻¹. (3.6) By the same manner, we get 𝐵⁻¹𝑋₁=𝑋₁𝐵⁻¹. (3.7) Further, assume that for each 𝑘 it is satisfied that 𝐴𝑋_𝑘−1=𝑋_𝑘−1𝐴,𝐵𝑌_𝑘−1=𝑌_𝑘−1𝐵,𝐴⁻¹𝑌_𝑘−1=𝑌_𝑘−1𝐴⁻¹,𝐵⁻¹𝑋_𝑘−1=𝑋_𝑘−1𝐵⁻¹. (3.8) Now, by induction, we will prove 𝐴𝑋_𝑘=𝑋_𝑘𝐴,𝐵𝑌_𝑘=𝑌_𝑘𝐵,𝐴⁻¹𝑌_𝑘=𝑌_𝑘𝐴⁻¹,𝐵⁻¹𝑋_𝑘=𝑋_𝑘𝐵⁻¹. (3.9) Since the two matrices 𝐴,𝐵 are normal and using the equalities (3.8), therefor 𝐴𝑋_𝑘=𝐴𝐼+𝐴^∗𝑌^−𝑛_𝑘−1𝐴=𝐴+𝐴𝐴^∗𝑌^−𝑛_𝑘−1𝐴=𝐴+𝐴^∗𝐴𝑌^−𝑛_𝑘−1𝐴=𝐴+𝐴^∗𝑌^−𝑛_𝑘−1=𝐴𝐴𝐼+𝐴^∗𝑌^−𝑛_𝑘−1𝐴𝐴=𝑋_𝑘𝐴. (3.10) Similarly 𝐵𝑌_𝑘=𝑌_𝑘𝐵. (3.11) By using the conditions 𝐴⁻¹𝐵=𝐵𝐴⁻¹,𝐴⁻¹𝐵^∗=𝐵^∗𝐴⁻¹ and the equalities (3.8), we have 𝐴⁻¹𝑌_𝑘=𝐴⁻¹𝐼+𝐵^∗𝑋^−𝑚_𝑘−1𝐵=𝐴⁻¹+𝐴⁻¹𝐵^∗𝑋^−𝑚_𝑘−1𝐵=𝐴⁻¹+𝐵^∗𝐴⁻¹𝑋^−𝑚_𝑘−1𝐵=𝐴⁻¹+𝐵^∗𝑋^−𝑚_𝑘−1𝐴⁻¹𝐵=𝐴⁻¹+𝐵^∗𝑋^−𝑚_𝑘−1𝐵𝐴⁻¹=𝐼+𝐵^∗𝑋^−𝑚_𝑘−1𝐵𝐴⁻¹=𝑌_𝑘𝐴⁻¹. (3.12) Also, we can prove 𝐵⁻¹𝑋_𝑘=𝑋_𝑘𝐵⁻¹. (3.13) Therefore, the equalities (3.2) are true for all 𝑠=0,1,2,….

Corollary 3.3. From Lemma 3.2, we have 𝐴𝑌_𝑠^−𝑛=𝑌_𝑠^−𝑛𝐴,𝐵𝑋_𝑠^−𝑚=𝑋_𝑠^−𝑚𝐵, (3.14) where {𝑋_𝑠},{𝑌_𝑠},𝑠=0,1,2,…, are determined by Algorithm 3.1.

Lemma 3.4. For the Sys. (1.1), we have 𝑋_𝑠𝑋_𝑠+1=𝑋_𝑠+1𝑋_𝑠,𝑌_𝑠𝑌_𝑠+1=𝑌_𝑠+1𝑌_𝑠, (3.15) where {𝑋_𝑠},{𝑌_𝑠},𝑠=0,1,2,…, are determined by Algorithm 3.1.

Proof. Since 𝑋₀=𝑌₀=𝐼, then 𝑋₀𝑋₁=𝑋₁𝑋₀,𝑌₀𝑌₁=𝑌₁𝑌₀.
By using the equalities (3.14), we have 𝑋₁𝑋₂=𝐼+𝐴^∗𝐴𝐼+𝐴^∗𝑌₁^−𝑛𝐴=𝐼+𝐴^∗𝐴+𝐴^∗𝑌₁^−𝑛𝐴+𝐴^∗𝐴𝐴^∗𝑌₁^−𝑛𝐴=𝐼+𝐴^∗𝐴+𝐴^∗𝑌₁^−𝑛𝐴+𝐴^∗𝐴𝑌₁^−𝑛𝐴^∗𝐴=𝐼+𝐴^∗𝐴+𝐴^∗𝑌₁^−𝑛𝐴+𝐴^∗𝑌₁^−𝑛𝐴𝐴^∗𝐴=𝐼+𝐴^∗𝑌₁^−𝑛𝐴𝐼+𝐴^∗𝐴=𝑋₂𝑋₁. (3.16) Similarly we get 𝑌₁𝑌₂=𝑌₂𝑌₁. (3.17) Further, assume that for each 𝑘 it is satisfied that 𝑋_𝑘−1𝑋_𝑘=𝑋_𝑘𝑋_𝑘−1,𝑌_𝑘−1𝑌_𝑘=𝑌_𝑘𝑌_𝑘−1. (3.18) Now, we will prove 𝑋_𝑘𝑋_𝑘+1=𝑋_𝑘+1𝑋_𝑘,𝑌_𝑘𝑌_𝑘+1=𝑌_𝑘+1𝑌_𝑘. (3.19) From the equalities (3.18), we have 𝑋^−𝑚_𝑘−1𝑋_𝑘^−𝑚=𝑋_𝑘^−𝑚𝑋^−𝑚_𝑘−1,𝑌^−𝑛_𝑘−1𝑌_𝑘^−𝑛=𝑌_𝑘^−𝑛𝑌^−𝑛_𝑘−1. (3.20) By using the equalities (3.14) and (3.20), we have 𝑋_𝑘𝑋_𝑘+1=𝐼+𝐴^∗𝑌^−𝑛_𝑘−1𝐴𝐼+𝐴^∗𝑌_𝑘^−𝑛𝐴=𝐼+𝐴^∗𝑌^−𝑛_𝑘−1𝐴+𝐴^∗𝑌_𝑘^−𝑛𝐴+𝐴^∗𝑌^−𝑛_𝑘−1𝐴𝐴^∗𝑌_𝑘^−𝑛𝐴=𝐼+𝐴^∗𝑌^−𝑛_𝑘−1𝐴+𝐴^∗𝑌_𝑘^−𝑛𝐴+𝐴^∗𝐴𝑌^−𝑛_𝑘−1𝑌_𝑘^−𝑛𝐴^∗𝐴=𝐼+𝐴^∗𝑌^−𝑛_𝑘−1𝐴+𝐴^∗𝑌_𝑘^−𝑛𝐴+𝐴^∗𝐴𝑌_𝑘^−𝑛𝑌^−𝑛_𝑘−1𝐴^∗𝐴=𝐼+𝐴^∗𝑌^−𝑛_𝑘−1𝐴+𝐴^∗𝑌_𝑘^−𝑛𝐴+𝐴^∗𝑌_𝑘^−𝑛𝐴𝐴^∗𝑌^−𝑛_𝑘−1𝐴=𝐼+𝐴^∗𝑌_𝑘^−𝑛𝐴𝐼+𝐴^∗𝑌^−𝑛_𝑘−1𝐴=𝑋_𝑘+1𝑋_𝑘. (3.21) By the same manner, we can prove 𝑌_𝑘𝑌_𝑘+1=𝑌_𝑘+1𝑌_𝑘. (3.22) Therefore, the equalities (3.15) are true for all 𝑠=0,1,2,….

Theorem 3.5. If 𝐴,𝐵 are satisfying the following conditions: (i)‖𝐴‖²(1+‖𝐴‖²)^𝑚−1<1/𝑚, (ii)‖𝐵‖²(1+‖𝐵‖²)^𝑛−1<1/𝑛, then the Sys. (1.1) has a positive definite solution (𝑋,𝑌), which satisfy 𝑋_2𝑠<𝑋<𝑋_2𝑠+1,𝑌_2𝑠<𝑌<𝑌_2𝑠+1,𝑠=0,1,2,…,max𝑋−𝑋_2𝑠,𝑋_2𝑠+1−𝑋<𝑞^𝑠‖𝐴‖²,max𝑌−𝑌_2𝑠,𝑌_2𝑠+1−𝑌<𝑞^𝑠‖𝐵‖², (3.23) where 𝑞=𝑛𝑚‖𝐴‖²‖𝐵‖²(1+‖𝐴‖²)^𝑚−1(1+‖𝐵‖²)^𝑛−1<1, {𝑋_𝑠},{𝑌_𝑠},𝑠=0,1,2,…, are determined by Algorithm 3.1.

Proof. For 𝑋₁,𝑌₁ we have 𝑋₁=𝐼+𝐴^∗𝐴>𝑋₀ and 𝑌₁=𝐼+𝐵^∗𝐵>𝑌₀.
Since 𝑋₁>𝑋₀,𝑌₁>𝑌₀ then 𝑋₁^−𝑚<𝑋₀^−𝑚,𝑌₁^−𝑛<𝑌₀^−𝑛 and 𝐵^∗𝑋₁^−𝑚𝐵<𝐵^∗𝑋₀^−𝑚𝐵,𝐴^∗𝑌₁^−𝑛𝐴<𝐴^∗𝑌₀^−𝑛𝐴, hence 𝐼=𝑋₀<𝑋₂=𝐼+𝐴^∗𝑌₁^−𝑛𝐴<𝐼+𝐴^∗𝑌₀^−𝑛𝐴=𝑋₁, 𝐼=𝑌₀<𝑌₂=𝐼+𝐵^∗𝑋₁^−𝑚𝐵<𝐼+𝐵^∗𝑋₀^−𝑚𝐵=𝑌₁, that is, 𝑋₀<𝑋₂<𝑋₁,𝑌₀<𝑌₂<𝑌₁. (3.24) We find the relation between 𝑋₂,𝑋₃,𝑋₄, and 𝑋₅ and the relation between 𝑌₂,𝑌₃,𝑌₄, and 𝑌₅.
Since 𝑋₀<𝑋₂,𝑌₀<𝑌₂, then 𝑋₃=𝐼+𝐴^∗𝑌₂^−𝑛𝐴<𝐼+𝐴^∗𝑌₀^−𝑛𝐴=𝑋₁, 𝑋₃=𝐼+𝐴^∗𝑌₂^−𝑛𝐴>𝐼+𝐴^∗𝑌₁^−𝑛𝐴=𝑋₂, 𝑌₃=𝐼+𝐵^∗𝑋₂^−𝑚𝐵<𝐼+𝐵^∗𝑋₀^−𝑚𝐵=𝑌₁, 𝑌₃=𝐼+𝐵^∗𝑋₂^−𝑚𝐵>𝐼+𝐵^∗𝑋₁^−𝑚𝐵=𝑌₂.
Since 𝑋₂<𝑋₃<𝑋₁,𝑌₂<𝑌₃<𝑌₁, then 𝑋₄=𝐼+𝐴^∗𝑌₃^−𝑛𝐴>𝐼+𝐴^∗𝑌₁^−𝑛𝐴=𝑋₂, 𝑋₄=𝐼+𝐴^∗𝑌₃^−𝑛𝐴<𝐼+𝐴^∗𝑌₂^−𝑛𝐴=𝑋₃, 𝑌₄=𝐼+𝐵^∗𝑋₃^−𝑚𝐵>𝐼+𝐵^∗𝑋₁^−𝑚𝐵=𝑌₂, 𝑌₄=𝐼+𝐵^∗𝑋₃^−𝑚𝐵<𝐼+𝐵^∗𝑋₂^−𝑚𝐵=𝑌₃.
Also since 𝑋₂<𝑋₄<𝑋₃,𝑌₂<𝑌₄<𝑌₃, then 𝑋₅=𝐼+𝐴^∗𝑌₄^−𝑛𝐴<𝐼+𝐴^∗𝑌₂^−𝑛𝐴=𝑋₃, 𝑋₅=𝐼+𝐴^∗𝑌₄^−𝑛𝐴>𝐼+𝐴^∗𝑌₃^−𝑛𝐴=𝑋₄, 𝑌₅=𝐼+𝐵^∗𝑋₄^−𝑚𝐵<𝐼+𝐵^∗𝑋₂^−𝑚𝐵=𝑌₃, 𝑌₅=𝐼+𝐵^∗𝑋₄^−𝑚𝐵>𝐼+𝐵^∗𝑋₃^−𝑚𝐵=𝑌₄.
Thus we get 𝑋₀<𝑋₂<𝑋₄<𝑋₅<𝑋₃<𝑋₁,𝑌₀<𝑌₂<𝑌₄<𝑌₅<𝑌₃<𝑌₁. (3.25) So, assume that for each 𝑘 it is satisfied that 𝑋₀<𝑋_2𝑘<𝑋_2𝑘+2<𝑋_2𝑘+3<𝑋_2𝑘+1<𝑋₁,𝑌₀<𝑌_2𝑘<𝑌_2𝑘+2<𝑌_2𝑘+3<𝑌_2𝑘+1<𝑌₁. (3.26) Now, we will prove 𝑋_2𝑘+2<𝑋_2𝑘+4<𝑋_2𝑘+5<𝑋_2𝑘+3 and 𝑌_2𝑘+2<𝑌_2𝑘+4<𝑌_2𝑘+5<𝑌_2𝑘+3.
By using the inequalities (3.26) we have 𝑋_2𝑘+4=𝐼+𝐴^∗𝑌^−𝑛_2𝑘+3𝐴>𝐼+𝐴^∗𝑌^−𝑛_2𝑘+1𝐴=𝑋_2𝑘+2,𝑋_2𝑘+4=𝐼+𝐴^∗𝑌^−𝑛_2𝑘+3𝐴<𝐼+𝐴^∗𝑌^−𝑛_2𝑘+2𝐴=𝑋_2𝑘+3. (3.27) Also we have 𝑌_2𝑘+4=𝐼+𝐵^∗𝑋^−𝑚_2𝑘+3𝐵>𝐼+𝐵^∗𝑋^−𝑚_2𝑘+1𝐵=𝑌_2𝑘+2,𝑌_2𝑘+4=𝐼+𝐵^∗𝑋^−𝑚_2𝑘+3𝐵<𝐼+𝐵^∗𝑋^−𝑚_2𝑘+2𝐵=𝑌_2𝑘+3. (3.28) Similarly 𝑋_2𝑘+5=𝐼+𝐴^∗𝑌^−𝑛_2𝑘+4𝐴<𝐼+𝐴^∗𝑌^−𝑛_2𝑘+2𝐴=𝑋_2𝑘+3,𝑋_2𝑘+5=𝐼+𝐴^∗𝑌^−𝑛_2𝑘+4𝐴>𝐼+𝐴^∗𝑌^−𝑛_2𝑘+3𝐴=𝑋_2𝑘+4. (3.29) Also we have 𝑌_2𝑘+5=𝐼+𝐵^∗𝑋^−𝑚_2𝑘+4𝐵<𝐼+𝐵^∗𝑋^−𝑚_2𝑘+2𝐵=𝑌_2𝑘+3,𝑌_2𝑘+5=𝐼+𝐵^∗𝑋^−𝑚_2𝑘+4𝐵>𝐼+𝐵^∗𝑋^−𝑚_2𝑘+3𝐵=𝑌_2𝑘+4. (3.30) Therefore, the inequalities (3.26) are true for all 𝑠=0,1,2,…; consequently the subsequences {𝑋_2𝑠},{𝑋_2𝑠+1},{𝑌_2𝑠}, and {𝑌_2𝑠+1} are monotonic and bounded. therefore they are convergent to positive definite matrices. To prove that the sequences {𝑋_2𝑠},{𝑋_2𝑠+1} have a common limit, we have 𝑋_2𝑠+1−𝑋_2𝑠=𝐴^∗𝑌^−𝑛_2𝑠𝐴−𝐴^∗𝑌^−𝑛_2𝑠−1𝐴=𝐴^∗𝑌^−𝑛_2𝑠−𝑌^−𝑛_2𝑠−1𝐴≤‖𝐴‖²𝑌^−𝑛_2𝑠𝑌^𝑛_2𝑠−1−𝑌^𝑛_2𝑠𝑌^−𝑛_2𝑠−1≤‖𝐴‖²𝑌^−𝑛_2𝑠𝑌^−𝑛_2𝑠−1𝑌^𝑛_2𝑠−1−𝑌^𝑛_2𝑠=‖𝐴‖²𝑌^−𝑛_2𝑠𝑌^−𝑛_2𝑠−1𝑌_2𝑠−1−𝑌_2𝑠^𝑛_𝑖=1𝑌^𝑛−𝑖_2𝑠−1𝑌^𝑖−1_2𝑠. (3.31) Since 𝐼<𝑌_𝑠<𝐼+𝐵^∗𝐵, then we have 𝑌_𝑠^−𝑛<1,𝑌_𝑠<1+‖𝐵‖². (3.32) Consequently 𝑋_2𝑠+1−𝑋_2𝑠<‖𝐴‖²𝑌_2𝑠−1−𝑌_2𝑠^𝑛_𝑖=1𝑌_2𝑠−1^𝑛−𝑖𝑌_2𝑠^𝑖−1<𝑛‖𝐴‖²1+‖𝐵‖^2𝑛−1𝑌_2𝑠−1−𝑌_2𝑠. (3.33) Also, to prove that the sequences {𝑌_2𝑠},{𝑌_2𝑠+1} have a common limit, we have 𝑌_2𝑠+1−𝑌_2𝑠=𝐵^∗𝑋^−𝑚_2𝑠𝐵−𝐵^∗𝑋^−𝑚_2𝑠−1𝐵=𝐵^∗𝑋^−𝑚_2𝑠−𝑋^−𝑚_2𝑠−1𝐵≤‖𝐵‖²𝑋^−𝑚_2𝑠𝑋^𝑚_2𝑠−1−𝑋^𝑚_2𝑠𝑋^−𝑚_2𝑠−1≤‖𝐵‖²𝑋^−𝑚_2𝑠𝑋^−𝑚_2𝑠−1𝑋^𝑚_2𝑠−1−𝑋^𝑚_2𝑠=‖𝐵‖²𝑋^−𝑚_2𝑠𝑋^−𝑚_2𝑠−1𝑋_2𝑠−1−𝑋_2𝑠^𝑚_𝑖=1𝑋^𝑚−𝑖_2𝑠−1𝑋^𝑖−1_2𝑠. (3.34) Since 𝐼<𝑋_𝑠<𝐼+𝐴^∗𝐴, then we have 𝑋_𝑠^−𝑚<1,𝑋_𝑠<1+‖𝐴‖². (3.35) Consequently 𝑌_2𝑠+1−𝑌_2𝑠<‖𝐵‖²𝑋_2𝑠−1−𝑋_2𝑠^𝑚_𝑖=1𝑋_2𝑠−1^𝑚−𝑖𝑋_2𝑠^𝑖−1<𝑚‖𝐵‖²1+‖𝐴‖^2𝑚−1𝑋_2𝑠−1−𝑋_2𝑠. (3.36) By using (3.36) in (3.33) and (3.33) in (3.36), we have 𝑋_2𝑠+1−𝑋_2𝑠<𝑛𝑚‖𝐴‖²‖𝐵‖²1+‖𝐴‖^2𝑚−11+‖𝐵‖^2𝑛−1𝑋_2𝑠−1−𝑋_2𝑠−2,𝑌_2𝑠+1−𝑌_2𝑠<𝑛𝑚‖𝐴‖²‖𝐵‖²1+‖𝐴‖^2𝑚−11+‖𝐵‖^2𝑛−1𝑌_2𝑠−1−𝑌_2𝑠−2. (3.37) Therefore 𝑋_2𝑠+1−𝑋_2𝑠<𝑞𝑋_2𝑠−1−𝑋_2𝑠−2<⋯<𝑞^𝑠𝑋₁−𝑋₀,𝑌_2𝑠+1−𝑌_2𝑠<𝑞𝑌_2𝑠−1−𝑌_2𝑠−2<⋯<𝑞^𝑠𝑌₁−𝑌₀. (3.38) Consequently the subsequences {𝑋_2𝑠},{𝑋_2𝑠+1} are convergent and have a common positive definite limit 𝑋. Also, the subsequences {𝑌_2𝑠},{𝑌_2𝑠+1} are convergent and have a common positive definite limit 𝑌. Therefore (𝑋,𝑌) is a positive definite solution of Sys. (1.1).

4. Numerical Examples

In this section the numerical examples are given to display the flexibility of the method. The solutions are computed for some different matrices 𝐴,𝐵 with different orders. In the following examples we denote by 𝑋,𝑌 the solutions which are obtained by iterative Algorithm 3.1, and 𝜖₁(𝑋)=‖𝑋−𝑋_𝑠‖, 𝜖₂(𝑋)=‖𝑋_𝑠−𝐴^∗𝑌_𝑠^−𝑛𝐴−𝐼‖, 𝜖₁(𝑌)=‖𝑌−𝑌_𝑠‖, and 𝜖₂(𝑌)=‖𝑌_𝑠−𝐵^∗𝑋_𝑠^−𝑚𝐵−𝐼‖.

Example 4.1. Consider Sys. (1.1) with 𝑛=5, 𝑚=5 and matrices 𝐴=0410−1801−2,𝐵=0−1−1−1−2−683−1. (4.1) By computation, we get 𝑋=,100019−5.999840−5.9998452.0729𝑌=1−2.49178×10⁻⁷−1.30736×10⁻⁶−2.49178×10⁻⁷21−1.30736×10⁻⁶,𝜆12.00002_𝑟(𝑋)={53.1277,17.9452,1},𝜇_𝑟(𝑌)={3.00002,1.00001,1}. (4.2) The results are given in the Table 1.

Table 1 

Example 4.2. Consider Sys. (1.1) with 𝑛=22, 𝑚=12 and matrices 𝐴=−0.10211240010421020,𝐵=−0.11212200010012010. (4.3) By computation, we get ,,𝜆𝑋=1.017170.02510930.0089920.004810470.02510931.11444−0.0619579−0.01276460.008992−0.06195791.124610.05232160.00481047−0.01276460.05232161.02878𝑌=1.0460.007308320.0167610.008766690.007308321.034040.01642780.03227750.0167610.01642781.016050.01284310.008766690.03227750.01284311.0343_𝑟(𝑋)={1.19574,1.07941,1.00868,1.00118},𝜇_𝑟(𝑌)={1.08177,1.04193,1.0064,1.0003}. (4.4) The results are given in the Table 2.

Table 2 

5. Conclusion

In this paper we considered the system of nonlinear matrix equations (1.1) where 𝑛,𝑚 are two positive integers. We achieved the general conditions for the existence of a positive definite solution. Moreover, we discussed an iterative algorithm from which a solution can always be calculated numerically, whenever the system is solvable. The numerical examples included in this paper showed the efficiency of the iterative algorithm which is described.