# On the Solution of a Class of Nonlinear Systems Governed by anM-Matrix

Authors: Woula Themistoclakis; Antonia Vecchio
Journal: Discrete Dynamics in Nature and Society
 (2012)

---

## Abstract

We consider a weakly nonlinear system of the form(I+φ(x)A)x=p, where φ(x) is a real function of the unknown vector x, and (I+φ(x)A) is an M-matrix. We propose to solve it by means of a sequence of linear systems defined by the iteration procedure (I+φ(xr)A)xr+1=p, r=0,1,…. The global convergence is proved by considering a related fixed-point problem.


---

## Body

## 1. Introduction

In recent years [1, 2], we approached the numerical resolution of the following nonstandard integrodifferential problems arising in the study of the kinetic theory of dusty plasmas:(1.1)ν(t)y(t)=∫0∞k(s)y(s)ds(D(t)y′(t))′+p(t),t>0,y′(0)=0,y(∞)=0,
whose discretization, by means of a difference scheme and a quadrature rule [2], leads to a particular kind of nonlinear system of equations. Solving it by means of a fixed-point (FP) iteration process, we noted that such a procedure seems to globally converge, that is, it converges independently of the choice of the starting point. Our aim, here, is to explain the reason of this “nice” behavior. Throughout the paper, the notation  M≥0, with M∈RN×K, means that each element of M is nonnegative, whereas M>0 means M≥0 and M≠0. The nonlinear system under investigation is(1.2)(I+φ(x)A)x=p,
where I is the N-dimensional identity matrix, x∈RN is an unknown vector, and we make the following assumptions which will be sufficient for the global convergence of the iterative process we propose (see Theorem 3.3):(i)
p∈RNwith  p>0;(ii)
A∈RN×N  has nonnegative main diagonal and nonpositive off-diagonal entries (Z-matrix);(iii)
A is rowwise weakly diagonal dominant;(iv)
φ:RN→R is a differentiable function such that φ(x)≥0, for all x≥0 and it is homogeneous of degree γ≤1, that is, for all α>0φ(αx)=αγφ(x);(v)
grad(φ(x))>0for  allx≥0, where grad(φ(x)) denotes the (row vector) gradient of φ(x).Systems of type (1.2) fall into the class of weakly nonlinear systems (see, e.g., [3, 4]) and they also arise from the discretization, by finite-difference methods, of mathematical models less cumbersome than (1.1), as, for instance, the differential problems of the form:(1.3)y(t)+b(y)a(t)y′(t)=p(t),t≥0,y(0)=0.A central rule in our results is played by the matrix(I+φ(x)A) appearing in (1.2) which, as it can be immediately seen from the assumption (ii)‒(iv), is an M-matrix for  allx≥0 (see, e.g., [5, page 137, prop. M35]).The organization of the paper is the following. In Section2, we investigate a particular class of one variable functions f such that the FP iterations tr+1=f(tr) converge independently of the choice of the starting value t0. By exploiting this result jointly with the properties of the matrix I+φ(x)A, we are lead to state that, under the assumptions (i)–(v), the FP iteration procedure (I+φ(xr)A)xr+1=p globally converges to a solution of (1.2). This constitutes our main result which is contained in Section 3. Finally, some numerical experiments showing the sharpness of the required conditions and the performance of the iteration scheme are reported in Section 4.


## 2. The One Variable Case

In order to prove our main theorem, we need the following result on a particular class of numerical sequences.Theorem 2.1.
Letf be a one variable scalar function such that(a)
f:t∈X⊆[0,+∞)→(0,Mf), with Mf>0, X closed and f(X)⊂X;(b)
f continuous;(c)
tf(t) strictly increasing.Then, starting from any t0∈X, the sequence obtained by the functional iteration process
(2.1)tr+1=f(tr),r=0,1,2,…,
converges to t*, a fixed point of f.In this section, we are going to prove this theorem.Denote byg the function g(t)=t-f(t), then we have(2.2)g(tr)=tr-f(tr)=tr-tr+1.
Of course, if there exists a point tr such that g(tr)=0, then tr=tr+n for all n=0,1,… and the statement is trivial. Similarly, if the sequence g(tr) is ultimately nonnegative (nonpositive), we have that the sequence tr is ultimately decreasing (increasing) and it trivially converges by hypothesis (a). Therefore, from now on, we assume that the sequence g(tr) satisfies(2.3)∀r≥0,∃m>r:g(tr)g(tm)<0.
The following lemma holds.Lemma 2.2.
Ifg(tr)<0(resp.>0) and g(tr+m)>0(resp.<0),  m=1,…,n,  then
(2.4)tr<tr+m+2(resp.tr>tr+m+2),m=-1,0,1,…,n.Proof.
Let us prove (2.4) by induction on n≥1. Assume n=1, that is, we are assuming that
(2.5)g(tr)<0,g(tr+1)>0,
and we want to prove
(2.6)tr<tr+m+2,m=-1,0,1.
From the first of (2.5) and (2.2), we immediately get (2.6) with m=-1, that is, we have
(2.7)tr<tr+1,
and so, from the hypothesis (c), we deduce
(2.8)trf(tr)<tr+1f(tr+1),
which, recalling that f(tr)=tr+1, implies
(2.9)trtr+1<tr+1tr+2,
and, therefore, (2.6) holds with m=0, that is,
(2.10)tr<tr+2.
By using (c) once again, (2.10) gives trf(tr)<tr+2f(tr+2) or equivalently
(2.11)trtr+1<tr+2tr+3.
But the second of (2.5) implies tr+2<tr+1, so that from (2.11), we obtain
(2.12)trtr+2<trtr+1<tr+2tr+3,
which assures that (2.6) holds for m=1 too.
Now, we assume that the lemma holds form=n, that is, suppose that
(2.13)g(tr)<0,g(tr+m)>0,m=1,…,n,
implies
(2.14)tr<tr+m+2,m=-1,0,1,…,n.
In order to prove the lemma for m=n+1, assume
(2.15)g(tr)<0,g(tr+m)>0,m=1,…,n+1.
Of course (2.14) holds, and the second of (2.15) implies
(2.16)tr+1>tr+2>⋯>tr+n+2.
But from (2.14) with m=n and (c), we have
(2.17)trtr+1<tr+n+2tr+n+3,
which in view of (2.16) leads to
(2.18)trtr+n+2<tr+n+2tr+n+3.
Hence, tr<tr+n+3, that is, (2.14) holds for m=n+1 too, and the desired result is proved.Now, let us partition the sequence{tr}r into two subsequences {tri}i and {tsj}j defined by(2.19)g(tri)<0,g(tsj)>0,i,j=1,2,…In order to state Theorem2.1, we are going to prove that these subsequences converge to the same limit. First of all, the following result holds.Lemma 2.3.
The sequences{tri}i and {tsj}j are strictly increasing and decreasing respectively.Proof.
Let us prove that the sequence{tri} is strictly increasing, the proof about {tsj} is analogous, and it is omitted. Consider ri such that g(tri)<0. Then, two cases may occur: (I) g(tri+1)<0or (II) g(tri+1)>0.(I)
In this case, it is clear thattri+1 belongs to the same subsequence {tri}i and we immediately have
(2.20)tri<tri+1=tri+1.(II)
From (2.3), there exists n≥2 such that tri+1=tri+n, that is, we have(2.21)g(tri)<0,g(tri+m)>0form=1,…,n-1,g(tri+n)<0.
Then, Lemma 2.2 with r=ri assures that
(2.22)tri<tri+m+2,m=-1,0,1,…,n-1,
which, for m=n-2≥0, gives
(2.23)tri<tri+n=tri+1.
The desired result comes out from (2.20) and (2.23).Lemma 2.4.
The sequences{tri}i and {tsj}j provide two separate sets, that is,
(2.24)tri≤tsj,∀i,j∈N.Proof.
Suppose, for example,ri<sj and precisely let sj=ri+k. If k=1, we get (2.24) directly from the assumption 0>g(tri)=tri-tri+1. If k≥2, observe that the numerical sequence {g(tri),g(tri+1),…,g(tri+k)} presents at least one change of sign, in view of (2.19). Let us examine the case of an unique change of sign, by supposing that there exists an integer l, with 0≤l≤k-1, such that
(2.25)g(tri+m)<0form=0,1,…,l,g(tri+l+m)>0form=1,…,k-l,
the other cases being analogous.
By Lemma2.3, we get
(2.26)tri≤tri+1≤⋯≤tri+l,tri+k≤tri+k-1≤⋯≤tri+l+1,
and by applying Lemma 2.2 with r=ri+l, we deduce
(2.27)tri+l<tri+l+m+2,m=-1,0,1,…,k-l.
Hence, by taking m=k-l-2≥-1, we conclude
(2.28)tri<tri+1<⋯<tri+l<tri+k<tri+k-1<⋯<tri+l+1.By the previous lemmas, we get the next result which completes the proof of the convergence of the FP iteration process.Theorem 2.5.
For any choice oft0∈X,  the sequence {tr}r defined in (2.1) converges.Proof.
It is sufficient to prove that the subsequences{tri} and {tsj} defined in (2.19) converge to the same limit, satisfying
(2.29)lim⁡itri=lim⁡jtsj=lim⁡rtr.
By the previous lemmas, we partitioned the sequence {tr} into two subsequences, {tri} and {tsj}, which are separate, strictly monotone, and bounded sequences. Hence, they converge and
(2.30)0≤sup⁡itri=lim⁡itri=:l1≤l2∶=lim⁡jtsj=inf⁡jtsj≤Mf.
Assume ab absurdo that
(2.31)l1<l2.
Let us consider the subsequence {tρi} of {tri} and the subsequence {tσj} of {tsj}, such that
(2.32)g(tρi-1)>0,g(tρi)<0,g(tσj-1)<0,g(tσj)>0.
In other words, tρi and tσj are those elements of {tri} and {tsj},  respectively, whose previous elements, in the main sequence {tr}, belong to the other subsequence {tsj} and {tri} respectively, that is, we have that
(2.33)∃k∈Nsuch  thattρi-1=tσk,∃h∈Nsuch  thattσj-1=tρh.
Of course, as {tρi} is a subsequence of the convergent sequence {tri}, there results
(2.34)lim⁡itρi=l1,
and the same is true for {tσj}, that is,
(2.35)lim⁡jtσj=l2,
where l1 and l2 are given in (2.30). On the other hand, recalling that tρi=f(tρi-1), from (b) and (2.33), we have
(2.36)l1=lim⁡itρi=lim⁡if(tρi-1)=lim⁡kf(tσk)=f(l2),l2=lim⁡jtσj=lim⁡jf(tσj-1)=lim⁡hf(tρh)=f(l1).
But by the hypothesis (c), the assumption (2.31) implies l1f(l1)<l2f(l2), that is, in view of (2.36), we have l1l2<l2l1, which is absurd. Thus, the result l1=l2=limrtr is achieved.


## 3. The Solution of the Nonlinear System

In this section, we come back to the nonlinear system (1.2). First of all, we recall that, under the hypotheses (ii)–(iv), for any x≥0, the matrix (I+φ(x)A)) appearing in (1.2) is an M-matrix and it satisfies (see, e.g., [5, prop. N38])(3.1)(I+φ(x)A)-1>0,∀x≥0.
Now, we define the function F as(3.2)F:q∈[0,+∞)⟶φ((I+qA)-1p),
and we are going to prove that it satisfies the properties stated in the following theorem, where γ is the homogeneity degree of φ appearing in (iv) and the classical notation e=(1,…,1)T∈RN is used.Theorem 3.1.
Assume that (i)–(v) hold.  Then,F is a continuous, differentiable function, and for any q≥0,  it satisfies
(3.3)F(q)∈(0,MF)withMF=‖p‖∞γφ(e),(3.4)F′(q)=-grad(φ((I+qA)-1p))(I+qA)-1A(I+qA)-1p,(3.5)s(q)∶=qγF(q)strictlyincreasing.Proof.
As the inverse of(I+qA) is defined for all q≥0, then F is clearly a continuous function. Moreover, it can be easily seen that
(3.6)ddq(I+qA)-1=-(I+qA)-1A(I+qA)-1,
hence, (3.4) is immediately true. In view of Theorem A.1 in [6], we observe that (ii) and (iii) assures
(3.7)∥(I+qA)-1∥∞≤1,∀q≥0.
Therefore, for any q≥0, the vector x(q):=(I+qA)-1p satisfies ∥x(q)∥∞≤∥p∥∞ and because of (v), we get F(q)=φ(x(q))≤φ(∥p∥∞e). Consequently, assertion (3.3) holds taking into account (iv). Finally, again in view of (iv), we have
(3.8)s(q)∶=qγF(q)=qγφ(x(q))=φ(qx(q))⟹s′(q)=grad(φ(qx(q)))(x(q)+qx′(q)),
where, recalling (3.6), it results
(3.9)x(q)+qx′(q)=[(I+qA)-1-(I+qA)-1qA(I+qA)-1]p=(I+qA)-2p.
Thus, by (v) and (3.1), we deduce (3.5) from
(3.10)s′(q)=grad(φ(qx(q)))(I+qA)-2p>0.Using this theorem, we can easily claim that the functionF in (3.2) has the same characteristics as the function f introduced in the previous section.Corollary 3.2.
Assume that (i)–(v) hold. Then, the functionF verifies all the hypotheses (a)–(c).Proof.
Assertions (a) and (b) directly follow from Theorem3.1. Moreover, in view of (3.5), (c) is immediately obtained by qF(q)=s(q)q1-γ, recalling that γ<1 implies q1-γ strictly increasing.Now, we are ready to state our main result about the convergence of the sequence{xr} defined as(3.11)xr+1=(I+φ(xr)A)-1p,r=0,1,…,x0≥0.Theorem 3.3.
Assume that the hypotheses (i)–(v) hold. Then, the sequence{xr} converges to a solution x* of (1.2).Proof.
Let us start from an arbitraryx0≥0 and put
(3.12)qr=φ(xr),r=0,1,…
From (3.11), we have that (3.12) can also be written as qr+1=F(qr),r=0,1,… Therefore, Corollary 3.2 assures that the sequence {qr} is convergent. Denoted by q* its limit, we set x*:=(I+q*A)-1p, and we obtain the statement by observing that
(3.13)lim⁡rxr+1=lim⁡r(I+φ(xr)A)-1p=lim⁡r(I+qrA)-1p=(I+q*A)-1p=x*.Remark 3.4.
Of course, in view of (3.7), any solution of (1.2) satisfies ∥x*∥∞≤∥p∥∞, but nothing seems to imply its uniqueness. Anyway, we conjecture it, because in a large variety of experimental tests (including also high-dimensional systems), we never found more than one solution. Of course, a sufficient condition for the uniqueness is ∥grad(φ(x))∥1∥A∥∞∥p∥∞<1,for any x≥0,which assures |F′(q)|<1, for any q≥0, by virtue of (3.4) and (3.7). Nevertheless, as it can be observed in the next section, it seems that the solution is unique even if such a condition is largely not satisfied.


## 4. Numerical Experiments

In order to verify our theoretical results and to check the performances of the FP iteration scheme on problem (1.2), we carried out a large variety of numerical experiments. Here, the most significant are reported. From the previous section, it appears clear that the iterative scheme (3.11), which furnishes a vector at each iteration, requires the same computational effort as (3.12), which instead provides scalar values. Hence, for the sake of simplicity, all the numerical tests reported here, are referred to (3.12), that is, we test the iterative procedure:(4.1)qr+1=φ((I+qrA)-1p)=F(qr).
In particular, in all the section, we assume that the convergence is reached whenever |qr+1-qr|<tol and we put tol=10-6.Let us consider the following two problems:(4.2)N=3,φ(x)=‖x‖2,A=(1-10-.51-.5-.5-12),p=(5110),γ=1;(4.3)N=3,φ(x)=x1exp⁡(x1.1x1+100x2+x3),x=(x1,x2,x3)T,A=(10-100-1015-.50-12),p=(.51002),γ=12
which, as it can be easily checked, satisfy the hypotheses of Theorem 3.3. In Tables 1 and 2, we report the number of iterations we performed in order to compute the fixed point q* of problems (4.2) and (4.3), respectively, for different choices of the starting value q0. Moreover, we also report the number of iterations performed by the Matlab routine fzero. This is given by the sum of the number of iterations to find an interval containing the zero (first addendum in Table 1) and the number of zero-finding iterations (second addendum in Table 1).Table 1FP iterations.15141415q0.2.41100fzero iterations.21+419+416+412+6Table 2FP iterations.1919191720q0.05.112100fzero iterations.ABABAB11ABFP +fzero iterations.1+131+151+123+42+13From these two tests, we observe that, according to Theorem3.3, the FP iteration method converges for any q0 in both the cases. Of course, when we know the interval where a zero lies, fzero does converge and it requires much less iterations, but in many cases this information is unavailable and, see Table 2, fzero fails to provide a value. In fact, at the third line of Table 2, the symbol AB stays for “aborted” and means that the MATLAB routine does not succeed in finding an interval containing a sign change. The last line of Table 2 shows that the best results can be obtained by a combined solution strategy, which starts with a certain number of FP iterations (first addendum) and takes the final steps using fzero iterations (second addendum). At the fifth column, we take the right balance between FP and fzero iterations, while in the remaining columns, we displayed the minimum number of FP iterations which are necessary to get a value from fzero. Thus, we can fruitfully exploit the global convergence of our FP iteration method in order to create more efficient combined strategy with fzero. Of course, this could be done also by means of other well-known methods like Newton’s one.In order to verify whether the convergence of the FP procedure could be ascribed to the Banach contraction theorem, we plotted the absolute value of the derivative of the functionF in (4.1) (not reported) and we noted that the starting value in bold, that is, q0=.2, q0=.4 in Table 1 and q0=.05, q0=.1 in Table 2, belongs to intervals where |F′(q)|>1. Hence we can assert that, at least in these cases, the convergence of the FP method is not “helped” by the contraction principle.As we already wrote in Remark3.4, from these numerical experiments and many others, starting from very different values of q0, the method seems to converge to the same fixed point. For example, in problems (4.2) and (4.3), we obtain q*≈4.7954 and q*≈2.6085, respectively, for a very large number of initial guesses (much larger than the ones reported in Tables 1 and 2). This numerical evidence let us conjecture that the assumptions in Theorem 3.3 are also sufficient to assure the uniqueness of the solution of problem (1.2).In Figures1 and 2, we plotted the first 100 elements of the error sequence {|qr-q*|}, where qr is computed by applying the FP iterations to the following two problems:(4.4)N=3,φ(x)=x1sin⁡2(20x1.1x1+x2+x3),q*=1.1210,A=(10-100-1015-.50-12),p=(.51002),γ=1;(4.5)N=1,φ(x)=x3/2,q*=0.0898,A=100,p=2,γ=32.
It can be easily seen that both the problems do not satisfy the hypotheses of Theorem 3.3. To be more precise, (4.4) does not fulfill (v), whereas (4.5) does not verify (iv), having γ>1. The figures clearly show that the method fails to converge, suggesting us that the assumptions of our main theorem are not too restrictive.Figure 1
The error sequence{|qr-q*|}r=1,…,100 for problem (4.4).Figure 2
The error sequence{|qr-q*|}r=1,…,100 for problem (4.5).Moreover, we point out that, in our experience, some problems arising in the the kinetic theory of dusty plasmas, are more cumbersome than (1.1), and therefore, the application of Newton or other available iterative procedures is not always simple, convenient, or possible. A comparison between FP and other iteration processes on problem of type (1.1) will be the subject of future investigations. Finally, in order to test the performance of our method for larger systems arising from the applications, we discretized a problem of type (1.1) with stepsize h on a sufficiently large interval (0,T), with T=Nh. In this way, we get a nonlinear system of type (1.2), which has dimension N. By applying FP iteration method to such a system for different values of N, with 102<N<105, we observe that the number of iterations does not vary with N and so with h, but only with the choice of the starting point q0. Moreover, we underline that in this case the linear systems arising from the FP iterations are tridiagonal and then very easy to be solved.



---

*Source: 412052-2012-06-08.xml*