# Design of High-Order Iterative Methods for Nonlinear Systems by Using Weight Function Procedure

Authors: Santiago Artidiello; Alicia Cordero; Juan R. Torregrosa; María P. Vassileva
Journal: Abstract and Applied Analysis
 (2015)

---

## Abstract

We present two classes of iterative methods whose orders of convergence are four and five, respectively, for solving systems of nonlinear equations, by using the technique of weight functions in each step. Moreover, we show an extension to higher order, adding only one functional evaluation of the vectorial nonlinear function. We perform numerical tests to compare the proposed methods with other schemes in the literature and test their effectiveness on specific nonlinear problems. Moreover, some real basins of attraction are analyzed in order to check the relation between the order of convergence and the set of convergent starting points.


---

## Body

## 1. Introduction

To find the solution of systems of nonlinear equationsF(x)=0, F:Rn→Rn, is a common, ancient, and important problem in mathematics and engineering (see, e.g., [1]). Previous experimentation on some of these applied problems shows that high-order methods, associated with floating point arithmetics of multiprecision, are very useful because they carry a clear reduction in the number of required iterations (see, e.g., [2, 3]).Over the last years there have been numerous contributions of different authors that have designed iterative methods trying to solve nonlinear systems, making several modifications to the classical schemes to accelerate the convergence and reduce the number of operations and functional evaluations at each step of the iterative process. The extension of the variants of Newton’s method for scalar problems described by Weerakoon and Fernando in [4], by Özban in [5], and by Gerlach in [6], to functions of several variables has been developed in [7–10]. All these processes have yield generating multipoint methods with Newton’s method as a predictor. However, a general procedure called pseudocomposition is designed in [11, 12] by using a generic predictor and the Gaussian quadrature as the corrector.On the other hand, one of the most used techniques to generate high-order methods for solving nonlinear equations or systems was introduced by Traub in [13], that is, the composition of two iterative methods of orders p1 and p2, respectively, to obtain a method of order p1p2. With the direct application of this technique, it is usually necessary to evaluate the nonlinear function and its associated Jacobian matrix to increase the order of convergence and the new method usually ends up being less efficient than the original ones.So, for designing a more efficient method, it is usual to estimate the new evaluation of the Jacobian matrix, in terms of Jacobian matrices previously used. With this in mind, the new method generally has a lower order of convergence thanp1p2 but the number of functional evaluations per iteration decreases.Now, let us consider the problem of finding a real zero of a functionF:Ω⊆Rn→Rn, that is, a solution ξ∈Ω of the nonlinear system of equations F(x)=0. In general, this solution must be approximated by means of iterative methods. The most known and widely used method is classical Newton’s method:(1)xk+1=xk-F′xk-1Fxk.It is known that this method has second order of convergence and it uses one functional evaluation of the vectorial function F(x) and also of the associated Jacobian matrix F′(x) at the previous iterate, to generate a new one.In the numerical section, we will use some multipoint iterative schemes to compare with ours. These have been developed by different authors: the first one is the fourth-order method designed by Jarratt in [14], whose iterative expression is(2)yk=xk-23F′xk-1Fxk,xk+1=xk-6F′yk-2F′xk-1M·3F′yk+F′xkF′xk-1Fxk.In the following, we will denote it by JM4. Let us also consider the fourth-order method designed by Sharma et al. in [15] whose iterative formula is(3)yk=xk-23F′xk-1Fxk,xk+1=xk-12-I+94F′yk-1MMMMMM+34F′xk-1F′ykM·F′xk-1Fxk,denoted by SHM4. We will also use some fifth-order methods to compare with ours: the scheme obtained by Vassileva [12], generalizing Jarratt’s method, improves this convergence by adding n functional evaluations (one functional evaluation of F in the second step of the process), being its iterative expression:(4)yk=xk-αF′xk-1Fxk,xk+1=xk-a1F′xk+a2F′yk-1M·b1F′xk+b2F′ykF′xk-1Fyk,where α=1, a1=1, a2=5, b1=-3, and b2=-1 and denoted by VM5, and the method designed by Sharma and Gupta in [16] whose iterative formula is(5)yk=xk-12F′xk-1Fxk,zk=xk-F′yk-1Fxk,xk+1=zk-aF′yk+bF′xk-1M·Fzk,where a=2 and b=-1. We denote this method as SHM5.In this work we have designed some iterative methods with orders of convergence four and five and presented an analysis in Section2. We compare them in Section 3 with other known processes on academic examples, with the aim that the results of this comparison will give us an idea of how robust our methods are. The stability is studied under the point of view of the real dynamics on an specific 2-dimensional problem. An applied problem, the partial differential equation of molecular interaction, is used to check the conclusions reached in previous sections.


## 2. Design and Analysis of Convergence

LetF:Ω⊆Rn→Rn, n>1, be a sufficiently differentiable function on a convex set Ω∈Rn and let ξ be a solution of the nonlinear system of equations F(x)=0. We propose the following iterative scheme:(6)yk=xk-αF′xk-1Fxk,xk+1=yk-HηkF′xk-1Fxk,where(7)ηk=m1I+m2F′yk-1F′xk+m3F′xk-1F′yk,with m1, m2, and m3 being real parameters.This scheme has been designed as a modified double-Newton method with a frozen evaluation ofF and F′ and a matrix weight function H(t) at the second step. We denote this class of methods (whose order of convergence is four under certain conditions) by MS4. As it has been stated, H is a matrix function with matrix variable. Specifically, if X=Rn×n denotes the Banach space of real square n×n matrices, then we can define H:X→X such that its Frechet derivatives satisfy(a)
H′(u)(v)=H1uv, where H′:X→L(X) and H1∈R,(b)
H′′(u,v)(w)=H2uvw, where H′′:X×X→L(X) and H2∈R.Let us note thatL(X) denotes the space of linear mappings of X on itself.To prove the local convergence of the proposed iterative method to the solutionξ of F(x)=0, we will use a Taylor series expansion of the involved functions around the solution, by assuming that Jacobian matrix F′(x) is nonsingular in a neighborhood of ξ,(8)Fξ+h=F′ξh+∑q=2p-1Cqhq+Ohp,where Cq=1/q![F′(ξ)]-1F(q)(ξ), q≥2.So,F′ can be expressed (in a neighborhood of ξ) as(9)F′ξ+h=F′ξI+∑q=2p-1qCqhq-1+Ohp,where I is the identity matrix. This technique, including the notation used, has been detailed in [17].Theorem 1.
LetF:Ω⊆Rn→Rn, n≥2, be a sufficiently differentiable function in an open convex set Ω and let ξ∈Ω be a solution of the system of nonlinear equations F(x)=0. One supposes that F′(x) is continuous and nonsingular at ξ. Let H:Rn×n→Rn×n be a sufficiently differentiable matrix weight function such that H(I)=1/3I, H1=3/4(m2-m3), and H2=3(m2-3m3)/4(m2-m3)3, with I being the identity matrix. Let one consider x(0) as an initial guess, sufficiently close to ξ. Then, the sequence {x(k)}k≥0 obtained using expression (6) converges to ξ with order of convergence four if α=2/3 and m1+m2+m3=1. Moreover, the error equation is(10)ek+1=34m2-m3213m22-26m2m3+45m32C23MM34m2-m32-C2C3+19C24ek4+Oek5,where Cq=1/q![F′(ξ)]-1F(q)(ξ), q≥2, and ek=x(k)-ξ.Proof.
From (8) and (9) we obtain, respectively,(11)Fxk=F′ξek+∑j=2pCjekj+Oekp+1,F′xk=F′ξI+∑j=2p-1Cjekj-1+Oekp.Then, forcing F′(xk)[F′(xk)]-1=[F′xk]-1F′(x(k))=I, the Taylor expansion of the inverse of the Jacobian matrix at the iterate x(k) is obtained as follows:(12)F′xk-1=I+∑j=2p-1Xjekj-1F′ξ-1+Oekp,where Xm=-∑j=2m‍iXm-i+1Cj, m=2,3,…. Then,(13)F′xk-1Fxk=ek+∑j=2pAjekj+Oekp+1,where A2=-C2 and As=Cs+∑j=3s‍Xs-j+2Cj-1+Xs, s=3,4,…. These expressions allow us to obtain(14)yk=xk-αF′xk-1Fxk=ξ+1-αek-αA+Oekp+1,where A=∑j=2p‍Ajekj. Using result (14) and the Taylor series expansion around ξ, we obtain(15)F′yk=F′ξI+∑j=1p-1Djekj+Oekp,F′yk-1=I+∑j=2p-1Y2ekjF′ξ-1+Oekp,where D1=2βC2, D2=2αC22+3β2C3, D3=-2αC2A3+6αβC3C2+4β5C4, D4=-2αC2A4+3α2C3C22+12αβC3A3-3αβ2C4C2+5β4C5, β=1-α, Y1=-D1, and Ys=Ds-∑j=2s-1‍YjDs-j, s=2,3,…. Thus, we obtain(16)uk=F′yk-1F′xk=I+∑j=1pEjekj+Oekp+1,where(17)E1=2C2+Y1,Es=s+1Cs+1+∑j=1s-1s-j+1Yj+1Cs-j+1+Ys+1,s=2,3,….Moreover,(18)vk=F′xk-1F′yk=I+∑j=1pFjekj+Oekp+1,where(19)F1=D1+X2,Fk=Dk+∑j+1k-1Xj+1Dk-j+Xk+1,k=2,3,….Then,(20)ηk=mI+∑j=2pNjekp+Oekp+1,where N1=m-1, Ni=m2Ei-1+m3Fi-1, i=2,3,…, and m=m1+m2+m3.
By using Taylor series expansion ofH around I,(21)Hηk=HI+H1ηk-I+12H2ηk-I2,we obtain the expression(22)Hηk=h1I+h2ek+h3ek2+h4ek3+h5ek4+Oek5,where h1=H(I)+H1N1+1/2H2N12, hj=H1Nj+1/2H2Mj-1, j=2,3,…, M1=N1N2+N2N1, M2=N1N3+N22+N3N1, M3=N1N4+N2N3+N3N2+N4N1, M4=N1N5+N2N4+N32+N4N2+N5N1, and M5=N1N6+N2N5+N3N4+N4N3+N5N2+N6N1.
So, we get(23)HηkF′xk-1Fxk=L1ek+L2ek2+L3ek3+L4ek4+Oek5,where L1=h1, Ls=∑j=1s-1‍hjAs+j-1, s=2,3,…. This allows us to obtain the error equation of this iterative scheme:(24)ek+1=W1ek+W2ek2+W3ek3+W4ek4+Oek5,where W1=β-L1 and Ws=-(αAs+Ls), s=2,3,….
To get order of convergence of at least four, we imposeW1=0, W2=0, and W3=0: from W1=1-α-L1=0, we get h1=1-α that gives us second order of convergence. Then, H(I)+H1N1+1/2H2N12=1-α and, as N1=0, H0=1-α is obtained.
Moreover, asW2=αA2+L2=0, by using h1=1-α and M1=0,(25)αA2+L2=αA2+h1A2+h2=αA2+1-αA2+h2=A2+h2=A2+H1N2+12H2M1=-C2+H1m2E1+m3F1=-1+H12αm2-m3C2.By replacing H1=1/2α(m2-m3) in (25), we get third order of convergence.
Forcing alsoW3=αA3+L3=0,(26)αA3+L3M=αA3+h1A3+h2A2+h3M=α+h1A3+h2A2+h3M=A3-h2C2+h3M=2C22-C3-H1N2+12H2M1C2MM+H1N3+12H2M2M=-2C22-C3-H1N2C2-N3+12H2M2M=1-32αC2MM+-2+2αm2m2-m3+2α2H2m2-m32C22.
Then, replacingα=2/3 and H2=3(m2-m3)/4(m2-m3)3 in (26), the order of convergence will be four.
Finally, the error equation is(27)ek+1=34m2-m3213m22-26m2m3+45m32C23MM34m2-m32-C2C3+19C24ek4+Oek5.So, the proof is completed.Starting from the general family of iterative methods (6) and under the hypothesis of Theorem 1, we develop two specific fourth-order iterative schemes.(i)
By requiringm1=m3=0 and m2=1, it is directly obtained from Theorem 1 that H(I)=1/3I and H1=H2=3/4. So, we get the following iterative expression denoted by MS4A:(28)yk=xk-23F′xk-1Fxk,xk+1=yk-HηkF′xk-1Fxk,Hηk=-124I+38F′yk-1M·F′xkF′yk-1F′xk.(ii)
The second iterative scheme is denoted by MS4B and is defined by imposing the following conditions:m1=1/4, m2=0, and m3=3/4 so that Theorem 1 implies that H(I)=1/3I, H1=-1, and H2=4. The resultant iterative expression that we call MS4B is(29)yk=xk-23F′xk-1Fxk,xk+1=yk-HηkF′xk-1Fxk,Hηk=5324I-3F′xk-1F′ykM+98F′xk-1F′ykF′xk-1M·F′yk.At this point we wonder if we can obtain a method of order of convergence higher than four slightly modifying the iterative scheme (6), just adding one functional evaluation of F in its second step:(30)yk=xk-αF′xk-1Fxk,xk+1=yk-HηkF′xk-1Fyk,where(31)ηk=m1I+m2F′yk-1F′xk+m3F′xk-1F′yk.We denote this class of methods by MS5. In the following, we prove that its order of convergence is five under certain conditions.Theorem 2.
LetF:Ω⊆Rn→Rn, n≥2, be a sufficiently differentiable function in an open convex set Ω and let ξ∈Ω be a solution of the system of nonlinear equations F(x)=0. One supposes that F′(x) is continuous and nonsingular at ξ. Let H:Rn×n→Rn×n be a sufficiently differentiable matrix weight function such that H(I)=I, H1=1/(m2-m3), and H2=m2-5m3/2(m2-m3)3, where I is the identity matrix. Let one consider x(0) as an initial guess, sufficiently close to ξ. Then, the sequence {x(k)}k≥0 obtained using the expression (30) converges to ξ with order of convergence five if α=1 and m1+m2+m3=1. Moreover, the error equation is(32)ek+1=2m22+14m32m2-m3C24+12C2C3C2-3C3C2ek5+Oek6,where Cq=1/q![F′(ξ)]-1F(q)(ξ), q≥2, and ek=x(k)-ξ.Proof.
Let us remark that the only difference between methods (30) and (6) is that, in its second step, F(x(k)) is replaced by F(y(k)). So, we need the Taylor series expansion of(33)Fyk=F′ξB1ek+B2ek2+B3ek3+B4ek4+B5ek5+Oek6,where B1=β=1-α, B2=(α+β2)C2, B3=-αA3+2αβC22+β3C3, B4=-αA4+α2C23-2αβC2A3+3αβ2C3C2+β4C4, and B5=-αA5-2αβC2A4+α2(C22A3+C2A3C2)+3α2βC3C22-3αβ2C3A3+4αβ3C4C2+β5C5. Then, we obtain the following expressions:(34)F′xk-1Fyk=R1ek+R2ek2+R3ek3+R4ek4+R5ek5+Oek6,HηkF′xk-1Fyk=S1ek+S2ek2+S3ek3+S4ek4+S5ek5+Oek6,where R1=B1, Rs=Bs+∑j=1s‍Xj+1Bj-s, S1=h1R1, and Sm=∑j=1m‍hjRm-j+1, s,m=2,3,….
So, the expression of the error equation in the last step is(35)ek+1=T1ek+T2ek2+T3ek3+T4ek4+T5ek5+Oek6,where T1=β-S1 and Tj=-αAj-Sj, j=2,3,…. Forcing T1=0, T2=0, T3=0, and T4=0, we obtain order of convergence five. Solving this system of equations we obtain the conditions that guarantee order of convergence five as they appear in the hypothesis of this theorem. Then, the error equation takes the form(36)ek+1=2m22+14m32m2-m3C24+12C2C3C2-3C3C2ek5+Oek6and the proof is completed.From the family of iterative methods (30) under the conditions imposed in Theorem 2, we get two particular fifth-order iterative schemes.(i)
By requiring thatm1=m3=1/2 and m2=0, conditions from Theorem 2 yield H(I)=I, H1=-2, and H2=10. The expression of the resulting iterative method that we call MS5A is(37)yk=xk-F′xk-1Fxk,xk+1=yk-HηkF′xk-1Fyk,Hηk=134I-72F′xk-1F′ykM+54F′xk-1F′ykF′xk-1M·F′yk.(ii)
The second scheme will be denoted by MS5B and results from setting the parametersm1=1/4, m2=0, and m3=3/4 in Theorem 2. Then, H(I)=I, H1=1, and H2=4. Then the iterative expression is(38)yk=xk-23F′xk-1Fxk,xk+1=yk-HηkF′xk-1Fyk,Hηk=14I+12F′xk-1F′ykM+14F′xk-1F′ykF′xk-1M·F′yk.It is interesting that, based on Theorem2, we can generate methods whose order of convergence is higher than 5. This process consists of adding a new step to scheme (30), keeping the weight function H(η(k)) unaltered. In this way, the iterative expression is(39)yk=xk-αF′xk-1Fxk,zk=yk-HηkF′xk-1Fyk,xk+1=zk-aHηkF′xk-1Fzk,where η(k)=m1I+m2F′yk-1F′(xk)+m3F′xk-1F′(y(k)).Assuming that hypotheses of Theorem2 are satisfied, the error in the second step of (30) can be expressed as E=z(k)-ξ=Mek5+Nek6+O(ek7). Then, the error equation of (39) takes the form(40)ek+1=Mek5+Nek6-aHηkF′xk-1Fzk=Mek5+Nek6-aI+h2ek+h3ek2M·I+X2ek+X3ek2Mek5+Nek6=M-aMek5+N-aNek6+Rek7+Oek8and the iterative method has order of convergence five. If a=1, the order of convergence of the iterative scheme (39) is seven.Based on this result, we state the following theorem which generalizes this procedure.Theorem 3.
Under the same hypothesis of Theorem2 and adding a=1, sequence {x(k)}k≥0 obtained using the following expression converges to ξ with order of convergence p+2:(41)yk=xk-αF′xk-1Fxk,zk=yk-HηkF′xk-1Fyk,uk+1=zk-aHηkF′xk-1Fzk,uk+1⋮xk+1=wk-aHηkF′xk-1Fwk,where η(k)=m1I+m2[F′yk]-1F′(xk)+m3[F′xk]-1F′(y(k)), and the order of convergence of the penultimate step is p with error equation E=w(k)-ξ=Mekp+Nekp+1+O(ekp+2).Proof.
In general, the error equation of the last step takes the form(42)ek+1=E-aHηkF′xk-1Fwk=Mekp+Nekp+1-aI+h2ek+h3ek2M·I+X2ek+X3ek2Mekp+Nekp+1=M-aMekp+N-aNekp+1+Rekp+2+Oekp+3.This shows that the order of convergence of the method will be p+2 if a=1.In order to calculate the different efficiency indices, it is necessary to take into account not only the order of convergence, but also the number of functional evaluations (n2: each Jacobian matrix F′ and n: each vectorial function F) per iteration and the amount of product-quotients per step. This is obtained by observing the iterative expression of the method and using the fact that the number of product-quotients needed for solving a set of m linear systems with the same coefficient matrix is 1/3n3+mn2-1/3n and the product between a matrix and a vector implies n2 products-quotients.With respect to the classical efficiency index, it is clear that one of all (both proposed and known) fourth-order methods isI4th=41/2n2+n and, for all fifth-order schemes, I5th=51/2n2+2n. It is also easy to check that I4th<I5th, for any n>2.In Table1, the computational efficiency indices are showed, and also the needed information is showed: the number of functional evaluations of matrices (FEM, n2 evaluations) or nonlinear functions (FEF, n evaluations), amount of linear systems with different Jacobian matrices (NS1, 1/3n3+mn2-1/3n products-quotients), the number of linear systems with the same Jacobian matrix (NS2, n2 products-quotients), and also how many times a matrix-vector product appears (MxV, n2 products).Table 1
Functional evaluations and product-quotients of the methods.
Methods
FEM
FEF
NS1
NS2
MxV
IEC
MS4A
2
1
2
1
2
4
(
1
/
(
(
2
/
3
)
n
3
+
7
n
2
+
(
1
/
3
)
n
)
)
MS4B
2
1
1
2
2
4
(
1
/
(
(
1
/
3
)
n
3
+
7
n
2
-
(
1
/
3
)
n
)
)
SHM4
2
1
2
1
1
4
(
1
/
(
(
2
/
3
)
n
3
+
6
n
2
-
(
1
/
3
)
n
)
)
JM4
2
1
2
0
1
4
(
1
/
(
(
2
/
3
)
n
3
+
5
n
2
-
(
1
/
3
)
n
)
)
MS5A
2
2
1
3
2
5
(
1
/
(
(
1
/
3
)
n
3
+
8
n
2
+
(
5
/
3
)
n
)
)
MS5B
2
2
1
3
2
5
(
1
/
(
(
1
/
3
)
n
3
+
8
n
2
+
(
5
/
3
)
n
)
)
SHM5
2
2
3
0
0
5
(
1
/
(
n
3
+
5
n
2
+
n
)
)
VM5
2
2
2
1
1
5
(
1
/
(
(
2
/
3
)
n
3
+
6
n
2
+
(
4
/
3
)
n
)
)Figure1 shows the computational efficiency index, for different sizes of the system, of fourth- (MS4A and MS4B) and fifth-order (MS5A and MS5B) proposed methods, joint with the known ones (JM4, SHM4, VM5, and SHM5), for comparison purposes. It can be observed that, in both cases, known methods have better computational efficiency indices for small sizes of the system, but in larger cases the new methods show better behavior.Figure 1
Computational efficiency index of fourth- and fifth-order methods.
(a)
(b)


## 3. Numerical Results

In this section we test the developed methods to check their effectiveness compared with other known ones. Here we can see the list of systems of nonlinear equations used in these numerical tests, conducted in Matlab R2010a by using variable precision arithmetics with 2000 digits of mantissa andx(k+1)-x(k)<10-500 or F(x(k+1))<10-500 as a stopping criterion. Consider(i)
F1(x1,x2)=(exp(x1)exp(x2)+x1cos(x2),x1+x2-1), ξ≈(3.47063096,-2,47063096)T;(ii)
F2(x1,x2)=(ln(x12)-2ln(cos(x2)),x1tanx1/2+x2-2), ξ≈(0.9548041416,6.5849814845)T;(iii)
F3(x1,x2)=(x1+exp(x2)-cos(x2),3x1-x2-sin(x2)), ξ=(0,0)T;(iv)
F4(x1,x2,x3)=(cos(x2)-sin(x1),x3x1-1/x2,exp(x1)-x32), ξ≈(0.90956949,0.66122683,1.57583414)T.We have two comparative tables in which iterative methods are grouped according to the order of convergence of the involved schemes: in Table2, we can observe the behavior of the fourth-order methods and in Table 3, the schemes of order 5 can be found. We perform numerical tests for each of the selected systems of nonlinear equations either for the proposed methods MS4A, MS4B, MS5A, and MS5B or for consolidated methods JM4, SHM4, VM5, and SHM5.Table 2
Numerical results obtained with the fourth-order schemes.
Methods
ξ
|
x
i
(
k
+
1
)
-
x
i
(
k
)
|
|
F
(
x
(
k
+
1
)
)
i
|
Iteration
ACOC
F
1 x(0)=(3,-2)T
MS4A
3.470631
2.7
e
-
1012
2.7
e
-
12
6
4.0000
−2.470631
2.7
e
-
1012
0
MS4B
3.470631
1.1
e
-
993
2.7
e
-
12
6
4.0000
−2.470631
1.1
e
-
993
0
SHM4
3.470631
1.0
e
-
1006
2.7
e
-
12
6
4.0000
−2.470631
1.0
e
-
1006
0
JM4
3.470631
2.0
e
-
1015
2.7
e
-
12
6
4.0000
−2.470631
2.0
e
-
1015
0
F
2 x(0)=(2,6)T
MS4A
0.954804
4.0
e
-
1691
1.7
e
-
2009
7
3.9192
6.584981
4.6
e
-
1691
1.8
e
-
11
MS4B
0.954804
2.1
e
-
1115
2.9
e
-
2008
7
3.9821
6.584981
1.3
e
-
1114
1.8
e
-
11
SHM4
0.954804
3.8
e
-
1427
6.1
e
-
2008
7
3.9601
6.584981
5.2
e
-
1427
1.8
e
-
11
JM4
0.954804
5.6
e
-
612
7.7
e
-
612
6
4.0266
6.584981
6.6
e
-
612
1.8
e
-
11
F
3 x(0)=(0.5,0.5)T
MS4A
0.000000
8.3
e
-
153
2.2
e
-
608
5
4.0000
0.000000
1.5
e
-
152
5.3
e
-
609
MS4B
0.000000
1.1
e
-
117
1.4
e
-
1868
5
3.9999
0.000000
1.7
e
-
117
8.5
e
-
1870
SHM4
0.000000
7.1
e
-
139
1.5
e
-
552
5
4.0000
0.000000
1.2
e
-
138
2.0
e
-
553
JM4
0.000000
1.4
e
-
662
4.5
e
-
662
6
4.0000
0.000000
3.0
e
-
662
1.7
e
-
662
F
4 x(0)=(0.8,0.5,1.4)T
MS4A
0.909569
1.0
e
-
466
7.4
e
-
1867
6
3.9937
0.661227
1.0
e
-
466
1.1
e
-
1863
1.575834
6.0
e
-
467
6.5
e
-
1864
MS4B
0.909569
2.4
e
-
297
1.8
e
-
1190
6
3.9996
0.661227
2.4
e
-
297
1.8
e
-
1185
1.575834
1.6
e
-
297
6.9
e
-
1186
SHM4
0.909569
1.2
e
-
382
3.8
e
-
1537
6
4.0013
0.661227
1.2
e
-
382
4.2
e
-
1527
1.575834
7.5
e
-
382
1.9
e
-
1527
JM4
0.909569
2.1
e
-
742
4.5
e
-
713
6
3.9870
0.661227
2.8
e
-
742
6.5
e
-
742
1.575834
3.5
e
-
743
4.2
e
-
742Table 3
Numerical results obtained with fifth-order schemes.
Methods
ξ
|
x
i
(
k
+
1
)
-
x
i
(
k
)
|
|
F
(
x
(
k
+
1
)
)
i
|
Iteration
ACOC
F
1 x(0)=(3,-2)T
MS5A
3.470631
1.4
e
-
850
2.7
e
-
12
5
5.0000
−2.470631
1.4
e
-
850
0
MS5B
3.470631
1.6
e
-
807
2.7
e
-
12
5
5.0000
−2.470631
1.6
e
-
807
0
SHM5
3.470631
2.6
e
-
805
2.7
e
-
12
5
5.0000
−2.470631
2.6
e
-
805
0
VM5
3.470631
4.1
e
-
2002
2.7
e
-
12
6
4.9999
−2.470631
4.1
e
-
2002
0
F
2 x(0)=(2,6)T
MS5A
0.954804
5.2
e
-
1386
1.7
e
-
2009
7
—
6.684981
6.7
e
-
1386
1.8
e
-
11
MS5B
0.954804
2.5
e
-
1649
2.5
e
-
2008
7
—
6.684981
5.4
e
-
1650
1.8
e
-
11
SHM5
0.954804
1.2
e
-
1001
2.0
e
-
1001
7
—
6.684981
7.3
e
-
1002
1.8
e
-
11
VM5
0.954804
2.9
e
-
813
5.6
e
-
2008
6
4.9994
6.684981
4.7
e
-
813
1.8
e
-
11
F
3 x(0)=(0.5,0.5)T
MS5A
0.000000
3.7
e
-
262
1.8
e
-
1306
5
5.0000
0.000000
5.7
e
-
262
5.3
e
-
1308
MS5B
0.000000
1.5
e
-
357
7.1
e
-
1784
5
4.9970
0.000000
2.9
e
-
357
1.8
e
-
1784
SHM5
0.000000
4.3
e
-
1880
1.3
e
-
1879
6
4.9978
0.000000
8.8
e
-
1880
4.6
e
-
1880
VM5
0.000000
8.9
e
-
323
5.1
e
-
1610
5
5.0000
0.000000
1.5
e
-
322
6.0
e
-
1611
F
4 x(0)=(0.8,0.5,1.4)T
MS5A
0.909569
2.2
e
-
429
2.3
e
-
1711
6
—
0.661227
1.1
e
-
428
1.2
e
-
1710
1.575834
1.1
e
-
428
1.6
e
-
1712
MS5B
0.909569
8.7
e
-
497
1.1
e
-
1985
6
—
0.661227
7.0
e
-
497
3.8
e
-
1985
1.575834
7.0
e
-
497
5.4
e
-
1985
SHM5
0.909569
1.7
e
-
1534
8.1
e
-
1534
7
—
0.661227
1.2
e
-
1533
3.7
e
-
1533
1.575834
1.1
e
-
1533
3.0
e
-
1533
VM5
0.909569
2.3
e
-
211
2.9
e
-
1060
5
4.9989
0.661227
2.3
e
-
211
4.8
e
-
1052
1.575834
1.5
e
-
211
2.1
e
-
1052The columns of the tables correspond from left to right to nonlinear system to be solved with the initial approximation, iterative method to be used, solutions or roots found, absolute value of the difference between the two last iterations (component by component),|xi(k+1)-xi(k)|, absolute value of each coordinate function evaluated at the last iteration, |F(x(k+1))i|, the number of iterations used in the process, and approximated computational order of convergence, ACOC, according to (see [8])(43)p≈ACOC=lnxk+1-xk/xk-xk-1lnxk-xk-1/xk-1-xk-2.Let us remark that the value of ACOC that is presented in Tables 2 and 3 is the last coordinate of vector ACOC when the variation between its values is small.Observing Tables2 and 3, we note that the approximated computational order of convergence confirms the theoretical results. The results in terms of accuracy in the approximation of the roots and number of iterations are kept in the range expected for the orders of convergence of the methods, as is apparent from a comparison with established methods. We believe that, in general, the proposed methods are competitive on each of the systems of nonlinear equations used.With respect to the extension to high-order methods, let us see in the following how the increasing order does not reduce the real region of good starting points, as it is usual in iterative methods. We have selected a rectangular region ofR2 that contains some of the solutions of system F1(x)=0 and have used the points of this region as starting ones for the fifth-order iterative methods and their partners of seventh and ninth order.So, in Figure2, the dynamical planes associated with known and proposed methods on F1(x) are showed. These planes have been generated by slightly modifying the routines described in [18]. In them, a mesh of 400×400 points has been used, 80 has been the maximum number of iterations involved, and 10-3 has been the tolerance used as a stopping criterium. Then, if a starting point of this mesh converges to one of the solutions of the system (marked with white stars), it is painted in the color assigned to the root which it has converged to. The color used is brighter when the number of iterations is lower. If it reaches the maximum number of iterations without converging to any of the roots, it is painted in black. At the sight of Figure 2, it can be concluded that the areas of converging starting points remain with slight variations, in spite of the increasing order of convergence; this makes this process especially interesting, as it can increase the speed of convergence to a solution with no need of being closer to it.Figure 2
Dynamical plane of increasing order methods onF1(x).
(a)
MS5A
(b)
MS7A
(c)
MS9A
(d)
MS5B
(e)
MS7B
(f)
MS9B


### 3.1. Molecular Interaction Problem

To solve the equation of molecular interaction (see [19])(44)uxx+uyy=u2,x,y∈0,1×0,1,ux,0=2x2-x+1,ux,1=2,u0,y=2y2-y+1,u1,y=2,we need to deal with a boundary value problem with a nonlinear partial differential equation of second order. To estimate its solution numerically, we have used central divided differences in order to transform the problem in a nonlinear system of equations, which is solved by using the proposed methods (of orders four and five) and the extensions of family MS5 up to order nine.The discretization process yields to the nonlinear system of equations,(45)ui+1,j-4ui,j+ui-1,j+ui,j+1+ui,j-1-h2ui,j2=0Mi=1,…,nx,j=1,…,ny,where ui,j denotes the estimation of the unknown u(xi,yj), xi=ih and yj=jk (with i=0,1,…,nx and j=0,1,…,ny) are the nodes in both variables, where h=1/nx and k=1/ny.In this case, we fixnx=ny=4, so a mesh of 5×5 is generated. As the boundary conditions give us the value of the unknown function at the nodes (x0,yj),(x4,yj) for all j and also at (xi,y0),(xi,y4) for all i, we have only nine unknowns that are renamed as(46)x1=u1,1,x2=u2,1,x3=u3,1,x4=u1,2,x5=u2,2,x6=u3,2,x7=u1,3,x8=u2,3,x9=u3,3.So, the system can be expressed as(47)Fx=Ax+ϕx-b=0,where(48)A=M-I0-IM-I0-IM,beingM=4-10-14-10-14,000000000000ϕx=h2x12,x22,…,x92T,and I is the 3×3 identity matrix and b=7/4,1,27/8,1,0,2,27/8,2,4T. Now, we will check the performance of the methods by means of some numerical tests, by using variable precision arithmetics of 2000 digits of mantissa. In Tables 4 and 5, we show the numerical results obtained for the problem of molecular interaction (45), with different initial estimations. We show, for the first three iterations, the residual of the function at the last iteration, F(x(k+1)), and the difference between the last iteration and the preceding one x(k+1)-x(k).Table 4
Numerical tests forx(0)=(1,…,1)T.
Method
∥

x
(
1
)
-
x
(
0
)

∥
∥

x
(
2
)
-
x
(
1
)

∥
∥

x
(
3
)
-
x
(
2
)

∥
∥

F
(
x
(
1
)
)

∥
∥

F
(
x
(
2
)
)

∥
∥

F
(
x
(
3
)
)

∥
MS4A
1.463
2.775
e
-
5
4.137
e
-
24
4.359
e
-
5
6.102
e
-
24
3.591
e
-
99
MS4B
1.463
1.004
e
-
4
2.39
e
-
21
1.603
e
-
4
3.536
e
-
21
1.335
e
-
87
MS5A
1.463
5.948
e
-
6
2.349
e
-
28
1.095
e
-
5
6.678
e
-
28
6.392
e
-
117
MS5B
1.463
9.502
e
-
7
1.903
e
-
30
3.374
e
-
6
2.708
e
-
30
4.787
e
-
126
MS7A
1.436
2.22
e
-
9
3.54
e
-
63
3.485
e
-
9
9.787
e
-
63
1.441
e
-
384
MS7B
1.463
3.642
e
-
10
1.924
e
-
66
1.141
e
-
9
2.795
e
-
66
3.04
e
-
405
MS9A
1.463
9.508
e
-
13
3.564
e
-
111
1.448
e
-
12
9.223
e
-
111
6.301
e
-
897
MS9B
1.463
1.79
e
-
13
1.416
e
-
115
4.545
e
-
13
2.084
e
-
115
1.772
e
-
933Table 5
Numerical tests forx(0)=(10,…,10)T.
Method
∥

x
(
1
)
-
x
(
0
)

∥
∥

x
(
2
)
-
x
(
1
)

∥
∥

x
(
3
)
-
x
(
2
)

∥
∥

F
(
x
(
1
)
)

∥
∥

F
(
x
(
2
)
)

∥
∥

F
(
x
(
3
)
)

∥
MS4A
25.2
0.5327
5.614
e
-
7
0.7255
8.027
e
-
7
1.314
e
-
30
MS4B
24.72
1.023
2.355
e
-
5
1.408
3.315
e
-
5
1.353
e
-
23
MS5A
25.2
0.5276
4.698
e
-
8
0.7184
7.861
e
-
8
5.271
e
-
37
MS5B
25.53
0.1917
9.957
e
-
11
0.2635
4.191
e
-
10
5.028
e
-
46
MS7A
25.63
0.08665
1.875
e
-
17
0.1176
4.203
e
-
17
7.064
e
-
111
MS7B
25.69
0.02205
1.189
e
-
21
0.03072
5.342
e
-
21
1.017
e
-
134
MS9A
25.7
0.01482
1.193
e
-
29
0.02016
3.181
e
-
29
5.147
e
-
245
MS9B
25.71
0.002606
4.285
e
-
36
0.003662
1.738
e
-
35
1.636
e
-
295We can observe in Tables4 and 5 that all the new methods converge to the solution of the problem that appears in Table 6. It can be noticed that the error of the test is lower when the order of the iterative method is higher, even at the first iterations. Indeed, if the initial estimation is far from the solution, the proposed methods converge with reasonable results. This is especially important as in real problems, where good initial estimations are not always known.Table 6
Approximated solution.
α
1
u
1,1
1.0259117…
u
2,1
1.2097139…
u
3,1
1.5167030…
u
1,2
1.2097139…
u
2,2
1.3877038…
u
3,2
1.6258725…
u
1,3
1.5167030…
u
2,3
1.6258725…
u
3,3
1.7642995…


## 3.1. Molecular Interaction Problem

To solve the equation of molecular interaction (see [19])(44)uxx+uyy=u2,x,y∈0,1×0,1,ux,0=2x2-x+1,ux,1=2,u0,y=2y2-y+1,u1,y=2,we need to deal with a boundary value problem with a nonlinear partial differential equation of second order. To estimate its solution numerically, we have used central divided differences in order to transform the problem in a nonlinear system of equations, which is solved by using the proposed methods (of orders four and five) and the extensions of family MS5 up to order nine.The discretization process yields to the nonlinear system of equations,(45)ui+1,j-4ui,j+ui-1,j+ui,j+1+ui,j-1-h2ui,j2=0Mi=1,…,nx,j=1,…,ny,where ui,j denotes the estimation of the unknown u(xi,yj), xi=ih and yj=jk (with i=0,1,…,nx and j=0,1,…,ny) are the nodes in both variables, where h=1/nx and k=1/ny.In this case, we fixnx=ny=4, so a mesh of 5×5 is generated. As the boundary conditions give us the value of the unknown function at the nodes (x0,yj),(x4,yj) for all j and also at (xi,y0),(xi,y4) for all i, we have only nine unknowns that are renamed as(46)x1=u1,1,x2=u2,1,x3=u3,1,x4=u1,2,x5=u2,2,x6=u3,2,x7=u1,3,x8=u2,3,x9=u3,3.So, the system can be expressed as(47)Fx=Ax+ϕx-b=0,where(48)A=M-I0-IM-I0-IM,beingM=4-10-14-10-14,000000000000ϕx=h2x12,x22,…,x92T,and I is the 3×3 identity matrix and b=7/4,1,27/8,1,0,2,27/8,2,4T. Now, we will check the performance of the methods by means of some numerical tests, by using variable precision arithmetics of 2000 digits of mantissa. In Tables 4 and 5, we show the numerical results obtained for the problem of molecular interaction (45), with different initial estimations. We show, for the first three iterations, the residual of the function at the last iteration, F(x(k+1)), and the difference between the last iteration and the preceding one x(k+1)-x(k).Table 4
Numerical tests forx(0)=(1,…,1)T.
Method
∥

x
(
1
)
-
x
(
0
)

∥
∥

x
(
2
)
-
x
(
1
)

∥
∥

x
(
3
)
-
x
(
2
)

∥
∥

F
(
x
(
1
)
)

∥
∥

F
(
x
(
2
)
)

∥
∥

F
(
x
(
3
)
)

∥
MS4A
1.463
2.775
e
-
5
4.137
e
-
24
4.359
e
-
5
6.102
e
-
24
3.591
e
-
99
MS4B
1.463
1.004
e
-
4
2.39
e
-
21
1.603
e
-
4
3.536
e
-
21
1.335
e
-
87
MS5A
1.463
5.948
e
-
6
2.349
e
-
28
1.095
e
-
5
6.678
e
-
28
6.392
e
-
117
MS5B
1.463
9.502
e
-
7
1.903
e
-
30
3.374
e
-
6
2.708
e
-
30
4.787
e
-
126
MS7A
1.436
2.22
e
-
9
3.54
e
-
63
3.485
e
-
9
9.787
e
-
63
1.441
e
-
384
MS7B
1.463
3.642
e
-
10
1.924
e
-
66
1.141
e
-
9
2.795
e
-
66
3.04
e
-
405
MS9A
1.463
9.508
e
-
13
3.564
e
-
111
1.448
e
-
12
9.223
e
-
111
6.301
e
-
897
MS9B
1.463
1.79
e
-
13
1.416
e
-
115
4.545
e
-
13
2.084
e
-
115
1.772
e
-
933Table 5
Numerical tests forx(0)=(10,…,10)T.
Method
∥

x
(
1
)
-
x
(
0
)

∥
∥

x
(
2
)
-
x
(
1
)

∥
∥

x
(
3
)
-
x
(
2
)

∥
∥

F
(
x
(
1
)
)

∥
∥

F
(
x
(
2
)
)

∥
∥

F
(
x
(
3
)
)

∥
MS4A
25.2
0.5327
5.614
e
-
7
0.7255
8.027
e
-
7
1.314
e
-
30
MS4B
24.72
1.023
2.355
e
-
5
1.408
3.315
e
-
5
1.353
e
-
23
MS5A
25.2
0.5276
4.698
e
-
8
0.7184
7.861
e
-
8
5.271
e
-
37
MS5B
25.53
0.1917
9.957
e
-
11
0.2635
4.191
e
-
10
5.028
e
-
46
MS7A
25.63
0.08665
1.875
e
-
17
0.1176
4.203
e
-
17
7.064
e
-
111
MS7B
25.69
0.02205
1.189
e
-
21
0.03072
5.342
e
-
21
1.017
e
-
134
MS9A
25.7
0.01482
1.193
e
-
29
0.02016
3.181
e
-
29
5.147
e
-
245
MS9B
25.71
0.002606
4.285
e
-
36
0.003662
1.738
e
-
35
1.636
e
-
295We can observe in Tables4 and 5 that all the new methods converge to the solution of the problem that appears in Table 6. It can be noticed that the error of the test is lower when the order of the iterative method is higher, even at the first iterations. Indeed, if the initial estimation is far from the solution, the proposed methods converge with reasonable results. This is especially important as in real problems, where good initial estimations are not always known.Table 6
Approximated solution.
α
1
u
1,1
1.0259117…
u
2,1
1.2097139…
u
3,1
1.5167030…
u
1,2
1.2097139…
u
2,2
1.3877038…
u
3,2
1.6258725…
u
1,3
1.5167030…
u
2,3
1.6258725…
u
3,3
1.7642995…


## 4. Conclusions

As far as we know, the weight functions procedure has been used only for designing iterative schemes to solve nonlinear equations. In this paper, by using matrix functions, the mentioned procedure is applied to obtain iterative methods for solving nonlinear systems. Fourth- and fifth-order methods are obtained and a technique for designing iterative schemes of any order is presented.By using different academic test problems and the discretization of a partial differential equation modeling the molecular interaction problem, we compare our methods with several known ones such as Jarratt’s method and Sharma’s method, some of them optimal in the context of nonlinear equations. The numerical tests confirm the theoretical results and show that our methods are more competitive than those used for comparing. In addition, the calculus of the computational efficiency index of the different schemes allows us to ensure that methods MS5A and MS5B are the best for nonlinear systems with big size.Finally, with respect to the extension to high-order methods, we show that the increasing order does not reduce the real region of good starting points, as it is usual in iterative schemes.



---

*Source: 289029-2015-08-05.xml*