Open access peer-reviewed chapter

Design, Analysis, and Applications of Iterative Methods for Solving Nonlinear Systems

By Alicia Cordero, Juan R. Torregrosa and Maria P. Vassileva

Submitted: October 22nd 2015Reviewed: May 5th 2016Published: October 19th 2016

DOI: 10.5772/64106

Downloaded: 920

Abstract

In this chapter, we present an overview of some multipoint iterative methods for solving nonlinear systems obtained by using different techniques such as composition of known methods, weight function procedure, and pseudo-composition, etc. The dynamical study of these iterative schemes provides us valuable information about their stability and reliability. A numerical test on a specific problem coming from chemistry is performed to compare the described methods with classical ones and to confirm the theoretical results.

Keywords

  • system of nonlinear equations
  • iterative methods
  • order of convergence
  • weight function procedure
  • stability
  • basin of attraction

1. Introduction

The problem of solving equations and systems of nonlinear equations is among the most important in theory and practice, not only of applied mathematics, but also in many branches of science, engineering, physics, computer science, astronomy, finance, etc. A glance at the literature shows a high level of contemporary interest.

The search for solutions of systems of nonlinear equations is an old, frequent, and important problem for many applications in mathematics and engineering (for example, see Refs. [15]).

The main goal of this chapter is to describe different methods for approximating a solution ξ of a system of nonlinear equations F(x) = 0, where F : Ω ⊆ ℝn → ℝn is a sufficiently differentiable function on the convex set Ω ⊆ ℝn. The most commonly used techniques are iterative methods, which, from an initial guess a sequence of iterates is built, were converging to the solution of the problem under some conditions. Although not as many as in the case of scalar equations, some publications have appeared in the recent years, proposing different iterative methods for solving nonlinear systems (see, for example, Refs. [68], among others). They have made several modifications to the classical methods to accelerate the convergence and to reduce the number of operations and functional evaluations per step of the iterative method. Newton’s method is the most used iterative technique for solving these kind of problems (see Ref. [9]), whose iterative expression is

x(k+1)=x(k)[F'(x(k))]1F(x(k)),k=0,1,E1

where F′ (x(k)) denotes the Jacobian matrix associated to function F on x(k).

We remember the concepts of order of convergence and efficiency index of an iterative scheme.

Definition 1.1. Let {x(k)}k≥0 be a sequence inn convergent to ξ. Then, the convergence is said to be linear, if there exist M, 0 < M < 1, and k0 ∈ ℕ such that x(k+1)ξMx(k)ξ,kk0, and of order p, p > 1, if there exist M, M > 0, and k0 ∈ ℕ such that x(k+1)ξMx(k)ξp,kko.

On the other hand, Ostrowski [10] introduced the efficiency index of an iterative method as p1/d, where p is the order of convergence and d is the number of functional evaluations per iteration. In the multidimensional case, it is more useful for the efficiency index to be defined as p1/(d+op), where op is the number of products-quotients per iteration.

The most direct technique is to adapt the methods designed for solving nonlinear equations to the multidimensional case. This process is easy only if, in the denominators of the iterative expression, does not appear any evaluation of the nonlinear function that describes the system. This is the case of Newton’s schemes and Newton-type methods coming from quadrature formulas, such as those described in Refs. [1118]. In Refs. [19] and [20], the authors designed a general procedure called pseudo-composition that allows to obtain predictor-corrector methods with high order of convergence. These multipoint schemes use any method as a predictor and a corrector step, where Gaussian quadrature is used.

On the other hand, other multipoint schemes have been developed by using different techniques: Adomian decomposition [2123], the replacement of the second derivative in one-point schemes by some approximation that yields multipoint iterative methods [8, 9], the Steffensen-type methods adapted to multidimensional case (see, for instance, Refs. [17, 24], among others). In all these papers, the references therein are also important.

A generally used technique for constructing iterative schemes is the composition of known methods. This technique was introduced by Traub [9]: if two iterative methods with orders of convergence p1 and p2, respectively, are composed, the resulting scheme has order p1p2. The use of this technique greatly increases the number of functional evaluations of F and F′. Therefore, it is necessary to avoid as much as possible such new evaluations by means of approximation techniques. In Ref. [25], the authors composed twice Newton’s scheme with itself ‘frozening’ the derivatives and, by means of undetermined coefficients method, obtained the following fifth-order iterative scheme

z(k)=y(k)5[F'(x(k))]1F(y(k)),x(k+1)=z(k)+15[F'(x(k))]1[15F(y(k))F(z(k))],

where y(k) is a Newton’s step. Let us observe that it only needs three functional evaluations and one Jacobian evaluation. Moreover, all the linear systems involved have the same matrix of coefficients. As a consequence, the efficiency index of this method is the best one, as far as we know. A similar procedure is used by Cordero et al. [26], getting Newton-Jarratt type methods of fifth and sixth orders. In addition, the following method described by Arroyo et al. [27] belongs to the class of Jarratt-type methods, but it has order of convergence five,

x(k+1)=y(k)+[F'(x(k))5F'(y(k))]1[F'(x(k))5F'(y(k))][F'(x(k))]1F(y(k)),

where y(k) is a Newton’s step.

In recent years, the technique of weight functions has also been developed, mainly for scalar equations. Weight functions are introduced in the iterative scheme to increase the order of convergence without increasing the number of functional evaluations. Among others, Sharma et al. [28] constructed a fourth-order scheme by using this procedure and, more recently, Artidiello et al. [7, 29] presented different families of high-order iterative methods by using matrix weight functions.

As we have previously mentioned, most of the iterative methods for nonlinear equations are not directly extendable to systems. However, in Refs. [6, 30], the authors present a general procedure to transform any scalar iterative method to the multidimensional case.

On the other hand, the dynamical analysis of an iterative method is becoming a trend in recent publications on iterative methods for scalar equations because it allows us to classify the different iterative formulas, not only from the point of view of its order of convergence, but also analyzing how these formulas behave as a function of the initial estimate that is taken. Another advantage of this analysis is to select the more stable elements of a parametric family whose members have the same order of convergence (see, for example, Ref. [31]). A first step in this direction on nonlinear systems was given by Cordero et al. [32], and a deeper analysis was made by Cordero et al. [33], studying the behavior of several methods on particular polynomial systems. In any case, the dynamical study in the multidimensional case is an emerging research topic with a promising future.

The rest of the chapter is organized as follows. In Section 2, we describe different techniques for designing iterative methods for nonlinear systems. In Section 3, we give some touches about the dynamic study of an iterative method for the scalar case and its extension to the multidimensional case. In order to check the introduced methods and compare them with other classical ones, in Section 4, we apply them on different test problems.

2. Design of the methods

In this section, we introduce three techniques for designing iterative methods for solving nonlinear systems of equations: pseudo-composition, weight function procedure, and a technique for extending scalar methods to the multidimensional case, in a non-trivial way.

The convergence results are going to be demonstrated by means of the n-dimensional Taylor expansion of the functions involved. Let F: Ωnnbe a sufficiently Frèchet differentiable function in Ω. By using the notation introduced by Cordero et al. [26], the qth derivative of F at un,q1, is the q-linear function F(q)(u):n×n×nnsuch that F(q)(u)(v1,,vq)n. It is easy to observe that: and F(q)(u)(v1,…,vq−1,·) ∈,(ℝn) and F(q)(u)(vσ1,,vσq)=F(q)(u)(v1,,vq), for all permutation σ of {1, 2, …, q}. We will use the notation:

F(q)(u)(v1,v2,,vq)=F(q)(u)v1v2vqandF(q)(u)vq1F(p)v(p)=F(q)(u)F(p)(u)vq+p1.E96321

For ξ+hnlying in a neighborhood of a solution ξ of the nonlinear system F(x) = 0 and assuming that the Jacobian matrix F′(ξ) is nonsingular, Taylor’s expansion can be applied, obtaining

F(ξ+h)=F'(ξ)[h+q=1p1Cqhq]+O[hp],

where Cq=(1q!)[F'(ξ)]1F(q)(ξ),q ≥ 2 We observe that CqFqnsince F(q)(ξ) ∈(ℝn × ⋯ × ℝn → ℝn). In addition, we can express the Jacobian matrix of F, F′ as

F'(ξ+h)=F'(ξ)[I+q=1p1qCqhq1]+O[hp],

where I is the identity matrix. Therefore, qCqhq1(ℝn). From this expansion we can conjecture that

[F(ξ+h)]1=[I+X2h+X3h2+X4h3+][F(ξ)]1+O[hp]

and taking into account that [F'(ξ+h)]1F'(ξ+h)=F'(ξ+h)[F'(ξ+h)]1=I, we obtain

X2=2C2,X3=4C223C3,X4=8C23+6C2C3+6C3C24C4, E5258525

We denote ek=x(k)ξthe error in the kth iteration. The equation ek+1=Lekp+O[ekp+1], where L is a p-linear function L(ℝn×…×ℝn, ℝn), is called error equation, and p is the order of convergence.

2.1. Pseudo-composition technique

We use the generic formulas of the Gaussian quadrature and develop families of predictor-corrector iterative methods, variants of Newton’s scheme, for solving nonlinear systems. Starting with any method of order p as a predictor and correcting over Gaussian quadrature, we will show that the final order of the obtained method will depend, among other things, on the order of the last two steps of the predictor. Let

y(k)=ξ+j=q5q1Mjekj+O[ek5q],z(k)=ξ+j=p5p1Njekj+O[ek5p],E800

be the penultimate and last steps of any iterative method with orders of convergence p and q, respectively. Taking this scheme as a predictor, we introduce the Gaussian quadrature as a corrector and we get four cases with the following iterative formulas:

(a)x(k+1)=y(k)2K1F(y(k)),(b)x(k+1)=z(k)2K1F(z(k)),(c)x(k+1)=y(k)2K1F(z(k)),(d)x(k+1)=z(k)2K1F(y(k)),E2

where for all cases K=i=1mωiF(ηi(k)),ηi(k)=12[(1+τi)z(k)+(1τi)y(k)],and ωi and τi are the weights and nodes, respectively, of the orthogonal polynomial of degree m, which defines the corresponding Gaussian quadrature. Then, ηi(k)is calculated by using the points obtained in the last two steps of the predictor.

To simplify the calculations, we use the following notation: j=q5q1Mjekj=A1(q)and j=q5q1Njekj=A2(p), where the subscripts in parentheses denote the value of the smallest power assumed by j in the sum. By using this notation, ηi(k)can be expressed as
ηi(k)=ξ+12(R+τiS)(q), where R=A2(p)+A1(q)and S=A2(p)A1(q). By expansion of F(y(k)),F(z(k))and F'(ηi(k))in the Taylor series around ξ, we obtain

F(y(k))=F(ξ)[A1(q)+C2A12(2q)+C3A13(3q)+C4A14(4q)]+O[ek5q],F(z(k))=F(ξ)[A1(p)+C2A12(2p)+C3A13(3p)+C4A14(4p)]+O[ek5p],F(ηi(k))=F(ξ)[I+B+Cτi+Dτi2+Eτi3]+O[ek4q],E4563217
whereB=C1R+34C3R2+12C4R3,C=C2S+34C3(RS+SR)+12C4(R2S+RSR+SR2), D=34C3S2+12C4(SRS+S2R+R2S)and E=12C4S3. We also introduce the following notation: i=1mωi=σand 1σi=1mωiτij=σjwith j = 1, 2, …, which will allow us to simplify the analysis of the convergence conditions of the described methods.

Now, we develop the expression K=i=1mωiF(ηi(k))appearing in Eq. (2) and we obtain K=σF(ξ)[I+K1(q)+K2(2q)+K3(3q)]+O[ek4q],whereK1(q)=C2(R+σ1S)(q), K2(2q)=34C3(R2+σ1(RS+SR)+σ2S2)(2q)andK3(3q)=[12C4(R3+σ1(R2S+RSR+SR2)+σ2(RS2+SRS+S2R)+σ2S3)(3q). Recalling that K1K=I, we get K1=σ1[I+K1(q)+K2(2q)+K3(3q)][F(ξ)]1+O[ek4q], where K1(q)=K1(q),K2(2q)=(K12K2)(2q), and K3(3q)=(K13+K2K1+K1K2K3)(3q). Therefore, considering case (a), we obtain for L=2K1F(y(k))the following expression:

L=2σA1(q)+2σ[(C2A1+K1)A1](2q)+2σ[(C2A12+K1C2A1+K2)A1](3q)+2σ[(C4A13+K1C3A12+K2C2A1+K3)A1](4q)+O[ek4q].E1478521

Since x(k+1)=y(k)L,the error equation can be expressed as

ek+1=A1(q)2σA1(q)2σ[(C2A1+K1)A1](2q)2σ[(C2A12+K1C2A1+K2)A1](3q)2σ[(C4A13+K1C3A12+K2C2A1+K3)A1](4q)+O[ek4q]

We note that if σ = 2 we obtain order of convergence at least 2q. The possibility of obtaining a convergence order greater than 2q depends on the expression (C2A1+K1)A1. We develop it and get (C2A1+K1)A1=σ1C2(A12)(2q)(1+σ1)C2(A2A1)(p+q). Then, σ1 = 0, the error equation is:

ek+1=C2(A2A1)(p+q)2σ[(C2A12+K1C2A1+K2)A1](3q)2σ[(C4A13+K1C3A12+K2C2A1+K3)A1](4q)+O[ek4q].E3

It is clear that the higher order of convergence in case (a) is min{p+q, 3q}. Finally, it is easy to prove that in case (b), the order of convergence of the resulting method is at least p. If σ = 2 the convergence order is p + q unless σ1 = 2 or C2 = 0, being in case the order of convergence is at least 2q + p. In cases (c) and (d) the convergence order is q, which is lower than the order of the predictor. It should be noticed that in case (b), two further functional evaluations are required and a new linear system per step must be solved. This causes the obtained method to be inefficient, from the standpoint of computational efficiency. In addition, we would like to remark that the order of convergence can be greater than p + q depending on expressions A1 and A2, which represent the errors of penultimate and last steps of the predictor method, and also of σ1 in case (b). The comments above allow us to state the following result, whose proof can be found by Cordero et al. [19], which establishes the order of convergence of the schemes that are obtained by using any method as a predictor of order p and eventually correcting it by the use of Gaussian quadrature, in case (a).

Theorem 2.1. Let F:Ωnnbe sufficiently differentiable function in Ω and F′(x) continuous and nonsingular at ξΩ, solution of the nonlinear system. Let y(k) and z(k) be the penultimate and last steps of orders q and p, respectively, of a certain iterative method. Taking this scheme as a predictor we get a new approximation x(k + 1) of ξ given by Eq. (2). Then,

  1. the methods of the obtained families have an order of convergence at least q,

  2. if σ = 2 is satisfied, then the order of convergence is at least 2q,

  3. if, also, σ1 = 0 the order of convergence is min{p + q, 3q}.

In Table 1, we show the values of σ and σ1 for some Gaussian quadrature.

In terms of computational efficiency, the most efficient methods are those which use fewer nodes and few functional evaluations, so we only consider case (a). Also, Theorem 2.1 shows that the order of convergence does not depend on the number of nodes, it only depends on the order of convergence of the penultimate and last step of the predictor method. Therefore, it is computationally more efficient to use one or two nodes. We note that the Gauss-Chebyshev quadrature does not fulfill the second condition of Theorem 2.1. Then, its order of convergence is q. The method obtained by using Gauss-Radau quadrature of one node does not fulfill the third condition but it verifies the second one; hence, its order of convergence is at least 2q. The remaining quadrature with nodes 1, 2, and 3 satisfies the conditions of Theorem 2.1 and the order of convergence is at least min{p + q, 3q}.

Number of nodesGaussian Quadrature
ChebyshevLegendreLobattoRadau
σσ1σσ1σσ1σσ1
1π020202−1
2π0202020
3π0202020

Table 1.

Quadrature formula used.

If we use case (a) and the Gauss-Legendre quadrature with 1 node or Gauss-Lobatto quadrature with one node such as corrector, we obtain the midpoint method, where K=2F(y(k)+z(k)2). In case of using Gauss-Radau quadrature with one node, we obtain Newton’s method (K=F(y(k))). Finally, if we use Gauss-Lobatto quadrature with two nodes or Gauss-Radau quadrature with two nodes such as corrector, we obtain trapezoidal method with K=2F(y(k)+z(k))and Noor’s scheme where K=8[F(y(k))+3F(y(k)+2z(k)3)], respectively.

2.2. Weight function procedure

The different methods obtained in the previous section are not optimal in the sense of Kung-Traub conjecture [34], when they are applied to scalar equations. By using the weight functions technique, we can increase the order of convergence of the designed methods without adding new functional evaluations.

We denote by X=n×nthe Banach space of all n × n real square matrices. The weight function in this context is a Frèchet differentiable function H:XXsatisfying:

  1. H(u)(v)=H1uv,Hbeing the first derivative of H,H:X(X),H1and (X)denotes the space of linear mapping from X to itself.

  2. H(u,v)(w)=H2uvw,Hbeing the second derivative of H,H:X×X(X),H2.

Then, the Taylor expansion of H around the identity matrix I, of size n × n, gives

H(ψ(k))H0(I)+H1(ψ(k)I)+12H2(ψ(k)I)2

Now, by using a relaxed Newton’s method as a predictor and a weight function procedure in the corrector step, we design the following families of two-point schemes:

z(k)=x(k)β[F(xk)]1F(xk),x(k+1)=z(k)2H(u(k))[i=1nωiF'(ηi(k))]1F(x(k))E4

where β is a non-zero real parameter and u(k)=1i=1mωi[F(x(k))]1i=1mωiF(ηi(k)). Here, we have used the same notations previously introduced. The following result establishes the order of convergence of this iterative scheme under some conditions of function H.

Theorem 2.2. Let ξ ∈ Ω be a zero of a sufficiently differentiable function F:Ωnn. Let us also suppose that the initial estimation x(0) is close enough to the solution ξ and F′(ξ) is nonsingular. The iterative methods (Eq. 4) have order of convergence four if

β=4(1+σ1)3(1+2σ1+σ2)

and a sufficiently differentiable function H is chosen satisfying the conditions

H0(I)=σ2I,H1=σ(1+2σ1+4σ12+3σ2)8(1+σ1)2H2=3σ(14σ12σ12+4σ134σ28σ1σ2+2σ12σ23σ22)8(1+σ1)4

and H‴ is a bounded operator, where σ and σj j = 1, 2, depend on the Gaussian quadrature used (see Ref. [35]).

Table 2 shows the values of σ1, σ2, and β, depending on the weights and the nodes of the orthogonal polynomials used in the Gaussian quadrature and also the corresponding values of the weight function and its derivatives, which are used in the iterative scheme.

QuadratureNo. of nodesσσ1σ2βH0(I)H1H2
Gauss-Chebyshev 1π004/3(π/2)Iπ/83π/8
Gauss-Legendre 12004/3I1/43/4
Gauss-Lobatto 12004/3I1/43/4
22012/3I−1/26
Gauss-Radau 12−100(1/2)I
2 2 0 1/3 1 I 0 2 

Table 2.

Nodes, H0(I), H1 and H2, for different Gaussian quadrature.

Let us remark that the class of methods designed allows the use of any Gaussian quadrature. In fact, when the technique of pseudo-composition was defined in Section 2.1 based also in Gaussian quadrature, the Chebyshev orthogonal polynomials could not be employed, as they did not verify the hypothesis of the main theorem. Nevertheless, with this new design, all the orthogonal polynomials can be applied and all of them derive optimal methods in the scalar case. In practice, we only use quadrature formulas with one and two nodes because, according to Theorem 2.2, the order of convergence is independent of the number of nodes.

For the one-dimensional case, according to the Kung-Traub’s conjecture [34], the obtained fourth-order methods are optimal (in case of Gauss-Radau with one node, classical Newton’s method is obtained). Other new methods are obtained in the rest of cases.

QuadratureNo. of nodesNameH(t)
Gauss-Chebyshev1GC1π16512t+15t2t2
Gauss-Legendre1GLe118(94t+3t2)
Gauss-Lobatto2GLo292132t+3t2
Gauss-Radau2GR2t22t+2

Table 3.

Notation and weight functions for different Gaussian quadrature formulas.

In Table 3, we show the weight functions used for each iterative scheme coming from the respective Gaussian quadrature rules. Let us note that the iterative method coming from Gauss-Lobatto with one node is the same as the one resulting from the application of Gauss-Legendre, also with one node. In this case, both coincide with the fourth-order procedure recently published by Sharma et al. [28]. Let us note that other weight functions should derive in other new schemes. For example, the expression of the method obtained by using Gauss-Legendre quadrature with one node is

y(k)=x(k)43[F'(x(k))]1F(x(k)),z(k)=12(x(k)+y(k)),x(k+1)=x(k)98[F'(x(k))]1F(x(k))+(12I38[F'(x(k))]1[F'(z(k))]1)[F'(x(k))]1F(x(k)).

2.3. Divided difference operator

In this section, we present a design published by Cordero et al. [30] of some families of parametric iterative methods for solving nonlinear equations by means of some known schemes and afterwards extend one of them to systems of nonlinear equations. For this purpose, we use Ostrowski’s [10] and Chun’s [36] methods with iterative schemes xk+1=ykf(xk)f(xk)2f(yk)f(yk)f(xk)and xk+1=ykf(xk)+2f(yk)f(xk)f(yk)f(xk), respectively, where yk is a Newton’s step. We propose a new family as a generalization of the previous methods in the following form:

yk=xkαf(xk)f(xk)xk+1=yk[f(xk)a1f(xk)+a2f(yk)+b1f(xk)+b2f(yk)f(xk)]f(yk)f(xk),E5

where α, a1, a2, b1, and b2 are real parameters. In the following result, we show which values of these parameters are necessary to guarantee at least order of convergence 4.

Theorem 2.3. Let f:Ibe a sufficiently differentiable function in an open interval I, such that ξ ∈ I is a simple root of the nonlinear equation f(x) = 0. If α=1,a2=a12(b22),b1=11a1and for all a1 and b2with, a10then sequence {xk}k0obtained from (Eq. 5) converges to ξ with local order of convergence at least four. In this case, the error equation is

ek+1=((5a1(b22)2)c23c2c3)ek4+O[ek5],E56324456

where ek=xkξand cq=(1q!)f(q)(ξ)f'(ξ),q2.

It is easy to prove this result by using any symbolic software as Wolfram Mathematica. To extend family (Eq. 5) to multivariate case, we need to rewrite its iterative expression in such a way that no functional evaluations of f remain at the denominator, as they will become vectors in the multidimensional case. So, let us consider that the first step of (Eq. 5) can be rewritten as f(xk)=1α(xkyk)f(xk). By using this, we can rewrite quotient f(yk)f(xk)as

f(yk)f(xk)=1αf[xk,yk]f(xk),E9000

where f[xk,yk]=f(xk)f(yk)xkykis the first-order divided difference. By using this transformation, the proposed family (Eq. 5) is fully extensible to several variables in the following way

y(k+1)=x(k)α[F'(x(k))]1F(x(k)),x(k+1)=y(k)(G1(x(k),y(k))+G2(x(k),y(k)))[F'(x(k))]1F(y(k)),G1(x(k),y(k))=[(a1+a2)Iαa2[F'(x(k))]1[x(k),y(k);F]]1,G2(x(k),y(k))=(b1+b2)Iαb2[F'(x(k))]1[x(k),y(k);F],E6

where [x(k),y(k);F]denotes the divided difference operator of F on x(k) and y(k), I is the identity matrix, and F(x(k))is the Jacobian matrix of the system.

Since the analysis of the local convergence is based on Taylor series expansion around the solution, we need to obtain the corresponding development of the divided difference operator. Let us denote by [x(k),y(k);F]the divided difference operator defined by Ortega and Rheinboldt [37] as the function [,;F]:Ω×Ωn×n(n)that satisfies [x,y;F](xy)=F(x)F(y),x,yΩ. To achieve this, we use the Genocchi-Hermite formula (see [37])

[x,x+h;F]=01F(x+th)dtE963214

and, by developing F(x+th)in Taylor series around x, we obtain

01F(x+th)dt=F(x)+12F(x)h+16F(x)h2+O[h3].E7

Assuming that

F(x)=F'(ξ)(ek+C2ek2+C3ek3+C4ek4)+O[ek5],F'(x)=F'(ξ)(I+2C2ek+3C3ek2+4C4ek3)+O[ek4],F''(x)=F'(ξ)(2C2+6C3ek+12C4ek2)+O[ek3],F'''(x)=F'(ξ)(6C3+24C4ek)+O[ek2],E8

and replacing these developments in the formula of Genocchi-Hermite, denoting the second point of the divided difference by y = x + h and the error at the first step by ek,y=yξ, we have [x,y;F]=F'(ξ)[I+C2(ek,yek)+C3ek2]+O[ek3]. In particular, if y is the approximation of the solution provided by Newton's method, i.e., h=xy=[F'(x)]1F(x), we obtain [x,y;F]=F'(ξ)[I+C2ek+(C22+C3)ek2]+O[ek3]. These tools allow us to prove the following result.

Theorem 2.4. Let f:Ωnnbe a sufficiently differentiable function in a convex set Ω, and ξΩbe a solution of the nonlinear system of equations F(x) = 0. Then, the sequence {x(k)}k0obtained by using expression (6) converges to ξ with local order of convergence at least four if α=1,a2=a12(b22),b1=11a1and for all a1 and b2with a1 ≠ 0. The error equation is

ek+1=((5a1(b22)2)C23C2C3)ek4+O[ek5],E456321

where ek=x(k)ξand Cq=(1q!)[F'(ξ)]1F(q)(ξ),q2.

Proof: By using Taylor expansion around ξ, we obtain:

F(x(k))=F'(ξ)(ek+C2ek2+C3ek3+C4ek4)+O[ek5],F'(x(k))=F'(ξ)(I+2C2ek+3C3ek2+4C4ek3)+O[ek4].

Let us consider [F'(x(k))]1=(I+X2ek+X3ek2+X4ek3)[F'(ξ)]1+O[ek4]. Forcing [F'(x(k))]1F'(x(k))=I, we get X2 = −2C2, X3 = 2C22−3C3 and X4=4C4+6C3C24C22+6C2C3. These expressions allow us to obtain for first step of iterative formula (6)

y(k)=ξ+(1α)ekα(A2ek2+A3ek3+A4ek4)+O[ek5],E9

where A2=C2X2,A3=C3C2X2X3and A4=C4C3X2C2X3X4. By using these results and the Taylor series expansion around ξ, we obtain

F(y(k))=F'(ξ)(B1ek+B2ek2+B3ek3+B4ek4)+O[ek5],

where B1=β, B2=(α + β2)C2, B3=−αA3 + 2αβC2A3 + 3αβ2C3C2 + β4C4, B4=−αA4 + α2C23−2αβC2A3 + 3αβ2C3C2 + β4C4 and β = 1 − α. We calculate the Taylor expansion of [x(k), y(k);F] by using Eq. (9), [ x(k), y(k);F]= F'(ξ)(I + D2ek + D3ek2 + D4ek3) + O[ek4], where D2= (2−α)C2, D3= αC22 + (3−3α + α2)C3 and D4= 2αC2C3 + α(3−2α)C3C2−(4−6α + 4α2α3)C4.

Then,

M=(a1+a2)Iαa2[F'(x(k))]1[x(k),y(k);F]=a1+E2ek+E3ek2+E4ek3+O[ek4],

where E2 = αa2C2, E3 = αC23 + α(α−3)C3 and E4 = 6αC2C3−2αC23−4C4 + 5α(2−α)C3C2 + (4−6α + 4α2α3)C3C4.

Thus, we obtain G1(x(k), y(k)) as the inverse of matrix M: G1(x(k), y(k))=I + Y2ek+Y3ek2+Y4ek3+O[ek4], where Y2=αa2a1C2, Y3=αa2a12[(αa23)C22+(α3)C3], Y4=αa2a13(8a1+3αa1a2+3αa2α2a23)C23and G2(x(k), y(k)) = b1+F2ek+F3ek2+F4ek3+F5ek4+O[ek5] where F2=αb2C2, F3=−αb2[2C22−(α−3)C3] and F4=b2[α(6−4α+α2)C4−6α(2−α)C3C2+4(α+1)C23−6(α+1)C2C3].

Finally, we obtain the error equation of the proposed method

ek+1=H1ek+H2ek2+H3ek3+H4ek4+O[ek5],E1003

where H1=1a1(1+a1(b11))(α1). If α = 1, then H1 = 0 and the error equation takes the form: ek+1=H2'ek2+H3'ek3+H4'ek4+O[ek5]where H2'=1a1(1+a1(b11))C2. We note that if b1=11/a1, then H2'=0We introduce this value of b1 and obtain the new form of the error equation ek+1=H3''ek3+H4''ek4+O[ek5]where H3''=(1/a12)[a2a12(b22)]C22. Finally, if a2=a12(b22), the error equation is:

ek+1=[(a1(b22)25)C23+C2C3]ek4+O[ek5]

and the proof is finished.□

By using the same hypothesis as in Theorem 2.4, the iterative scheme of the class (Eq. 6) takes the form:

y(k+1)=x(k)[F'(x(k))]1F(x(k)),x(k+1)=y(k)(G1(x(k),y(k))+G2(x(k),y(k)))[F'(x(k))]1F(y(k)),G1(x(k),y(k))=1a1[(1+a1b22a1)Ia1(b22)[F'(x(k))]1[x(k),y(k);F]]1G2(x(k),y(k))=1a1[(a1+a1b21)Ib2[F'(x(k))]1[x(k),y(k);F]].E10

In the following we propose some particular cases:

1. When a1 = 1, the iterative expression, being G(x(k),y(k))=G1(x(k),y(k))+G2(x(k),y(k))takes the form:

G(x(k),y(k))=[(b21)I(b22)[F'(x(k))]1[x(k),y(k);F]]1+b2Ib2[F'(x(k))]1[x(k),y(k);F]

and we have a family of schemes with interesting particular cases, among others,

a) If b2 = 2, Chun's method transferred to systems is obtained

x(k+1)=y(k)(I2[F'(x(k))]1[x(k),y(k);F])[F'(x(k))]1F(y(k))E1000

b) If b2 = 0 we get Ostrowski's scheme transferred to systems

x(k+1)=y(k)(I+2[F'(x(k))]1[x(k),y(k);F])1[F'(x(k))]1F(y(k)).E1001

2. When b2 = 0 for any a1 ≠ 0, the iterative expression of the parametric family is

G(x(k),y(k))=(32a1)I+2(a11)[F'(x(k))]1[x(k),y(k);F](12a1)I+2a1[F'(x(k))]1[x(k),y(k);F].

If we express −2(a1 − 1) = β we get the King's family transferred to systems

x(k+1)=y(k)(1+β)I+β[F'(x(k))]1[x(k),y(k);F](β1)I+(β2)[F'(x(k))]1[x(k),y(k);F][F'(x(k))]1F(y(k)).

3. When b2 = 1 for any a1 ≠ 0 and a1 ≠ 1 the iterative form, the last step of the method (Eq. 10) takes the form:

x(k+1)=y(k)(1a1[(1a1)Ia1[F'(x(k))]1[x(k),y(k);F]]1+1a1[(2a11)I[F'(x(k))]1[x(k),y(k);F]])[F'(x(k))]1F(y(k)).

3. Dynamic studies of some methods

In the last years, a new branch of the analysis of iterative methods for solving nonlinear equations or systems has taken relevance: the dynamical analysis of the rational functions associated with the fixed point operator associated with the iterative scheme on polynomials. By using complex or real dynamics techniques, the stability and reliability of a method can be checked. Indeed, if a parametric family of iterative procedures is considered, this kind of analysis allows selecting those elements of the class with better stability and also to know which ones behave chaotically even on the most simple functions, as low degree polynomials.

The use of these dynamical tools is very frequent on scalar iterative methods; see, for example, Refs. [3845] and the references therein, but in the multidimensional case, it is a starting area of research.

Now, let us recall some basic concepts on complex dynamics. Given a rational function, whereis the Riemann sphere, the orbit of a point z0is defined as {z0,R(z0),R2(z0),,Rn(z0),}. A point z*is called a fixed point of R(z)if it verifies that R(z*)=z*. Moreover, z* is called a periodic point of period p > 1 if it is a point such that Rp(z*)=z*but Rk(z*)z0, for each k < p. Moreover, a point z* is called pre-periodic if it is not periodic but there exists a k > 0 such that Rk(z*)is periodic.

There exist different types of fixed points depending on their associated multiplierTaking the associated multiplier into account, a fixed point z* is called: superattracting if, attracting if, repulsive if, and parabolic if. The fixed point operator of any iterative method on an arbitrary polynomial p(z) is a rational function. The fixed points of this rational function that do not correspond to the roots of the polynomial p(z) are called strange fixed points. On the other hand, a critical point z* is a point satisfying |R'(z*)|=0.

The basin of attraction of an attractor α is defined asThe Fatou set of the rational function R,is the set of points zwhose orbits tend to an attractor (fixed point, periodic orbit or infinity). Its complement inis the Julia set,. That means that the basin of attraction of any fixed point belongs to the Fatou set and the boundaries of these basins of attraction belong to the Julia set.

Some other concepts are a key fact in this kind of analysis, as the immediate basin of attraction of an attracting fixed point α (considered as a periodic point of period 1), as the connected component of the basin containing α. This concept is directly related with the existence of critical points, as can be seen in the following classical result, due to Fatou and Julia.

Theorem 3.1. Let R be a rational function. The immediate basins of attraction of any attracting periodic point hold, at least, a critical point.

If a scaling theorem can be established, the qualitative behavior of the class of iterative schemes is analyzed on a generic quadratic polynomial p(z)=z2c, its dynamics being analytically conjugated by affine transformations. Then, the fixed points of their associated rational function (being or not roots of the polynomial) are calculated and it is studied if they are stable (that is, if the successive iterations of the method converge to them) or if they are unstable, repelling the iterations near them. Finally, the calculation of critical points and their use as starting points of the iterative process will allow us, by applying Theorem 3.1, to find all the values of the parameter (that is, methods of the family) that do not converge to the roots of the polynomial, resulting in unstable behavior (attracting periodic orbits, chaos, etc.).

This kind of analysis has been developed by Cordero et al. [46] on the fourth-order biparametric Ostrowski-Chun family, whose iterative expression was shown in Eq. (5). In this manuscript, the strange fixed points and free independent critical points of the fixed point rational operator

Op(z,a1,b2)=z4(z+1)2(z2+4z+5)a1(b22b2(z3+4z2+5z+4)+2(z+1)2(z+2))z4(a1(b22)25)+z3(5a1(b22)14)2z2(2a1(b22)+7)z(a1(b22)+6)1E4568956

have been identified. As the behavior of the elements of the family depends on two parameters, only real values have been considered in the plane (a1, b2) in order to calculate the multipliers of the strange fixed points: z = 1 and

ex1=112(A3B6(C+D)),ex2=112(A3B+6(C+D)),ex3=112(A+3B6(C+D)),ex4=112(A+3B+6(C+D)),

where A, B, C, and D depend on both parameters a1 and b2. In Figure 1, the stability function of strange fixed point z = 1 is represented, |Op'(1,a1,b2)|. The orange region corresponds to real values of parameters a1 and b2 where z = 1 is attracting, meanwhile gray region includes values of a1 and b2 where this strange fixed point is repulsive, and therefore, the corresponding iterative methods have better behavior. Similar studies can be made on the rest of strange fixed points, searching for a more stable position of plane (a1, b2).

Figure 1.

Stability function of z = 1 as fixed point of operator Op(z, a1, b2).

On the other hand, as the dynamics of critical points could lead to a Fatou component, their associated parameter planes are represented to see the behavior of the method when the initial estimate is a critical point. In this case, z = −1 and cri,i=1,2,3,4,are free critical points of Op(z,a1, b2) such that cr1=1cr2and cr3=1cr4, that will be used to calculate the different parameter planes associated with the fixed point operator. As z = −1 is a pre-image of z = 1, only two critical points can be considered as free and independent: cr2 and cr4.

The parameter space of a free critical point is obtained by associating each point of the parameter plane with real values of free parameters a1 and b2, so every point on the plane represents a different member of the iterative family. These parameter planes (one of them can be seen in Figure 2) have been created using a vectorized version of the MATLAB code programs presented by Chicharro et al. [47], with 800 × 800 different combinations of a1 and b2. Black points correspond to the parameter values for which the associated iterative method does not converge to the conjugated values of the zeros of polynomial p(z) with a tolerance of 10−3 after 500 iterations, taking the same free independent critical point as the initial estimate.

Figure 2.

Parameter plane of operator Op(z, a1, b2).

Points shown in red in Figure 2 (and also simultaneously red in the rest of parameter planes corresponding to other free independent critical points of the rational function) correspond to the most stable methods of the family. In these terms, Figure 3 shows the dynamical plane associated with one of these stable elements of the family: each point with 800 × 800 mesh in the complex plane is an initial estimation for the iterative method. This point has an associated color depending on the convergence of the method after 80 iterations: if it converges to a fixed point that corresponds to a root of p(z), then a color (orange and blue) is assigned (see Figure 3a). If not, then it diverges or converges to any other element (attracting strange fixed point, periodic orbit, etc.), and it is painted in black (see Figure 3b). In the former case, the orbit of an initial point in a black region is marked in yellow.

Figure 3.

Dynamical plane of operator Op(z, a1, b2).

Selected values of parameters a1 and b2 in these red regions were tested numerically, not only on scalar equations, but also on multidimensional problems. These tests showed that good performance of the member of the family observed in the dynamical study could also be noticed in the multidimensional case.

However, very recently real multivariate dynamical tools have revealed also to be a reliable implement to analyze the stability of iterative methods, specially designed for solving nonlinear systems. In order to study the stability of the fixed points of vectorial rational functions associated with iterative schemes for solving nonlinear systems, a new tool was presented by Cordero et al. [33]. In it, the authors checked the consistence of this tool by applying it on known methods as Newton's and Traub's schemes and also on a family of parametric iterative procedures

y(k)=x(k)θ[F'(x(k))]1F(x(k)),z(k)=x(k)[F'(x(k))]1(θF(x(k))+F(y(k))),x(k+1)=x(k)[F'(x(k))]1(θF(x(k))+F(y(k))+F(z(k))),

Presented by Hueso et al. [48] and denoted by HMT.

In order to analyze, under a multidimensional point of view, the qualitative behavior of this family, the Cordero et al. [33] introduced some concepts of multidimensional real discrete dynamics that we will recall in the following; some of them are the natural extensions of the defined concepts in complex dynamics, but others will be defined, as being different from those.

Let us denote by G(x) the vectorial fixed-point function associated with the iterative method or family on polynomial p(x). The dynamical behavior of the orbit of a point of ncan be classified depending on its asymptotic behavior. In this way, a point x* is a fixed point of G if G(x*) = x*.

Theorem 3.2. [43, page 558] Let G:nnbe a function in. Assume x* is a period-k point. Let λ1,λ2,,λnbe the eigenvalues of G'(x*).

  1. If all the eigenvalues λjhave |λj|<1, then x* is attracting.

  2. If one eigenvalue λj0has |λj0|>1, then x* is unstable, that is, repelling or saddle.

  3. If all the eigenvalues λjhave |λj|>1, then x* is repelling.

Moreover, a fixed point is called hyperbolic if all the eigenvalues λjof G′(x*) have |λj|1. Indeed, if there is an eigenvalue λisuch that |λi|<1and also there exists an eigenvalue λjsuch that |λj|>1, the hyperbolic point is called saddle point.

To avoid the calculation of spectrum of G′(x*), the authors proposed by Cordero et al. [33] a result that, being consistent with the previous theorem, provided a practical tool to classify the stability of fixed points in many multidimensional cases.

Proposition 3.1. Let x* be a fixed point of G. then,

  1. if, for all {1,2,,n}, then x* is attracting.

  2. if, for all i,jϵ{1,2,,n}, then x* is superattracting.

  3. if, for all jϵ{1,2,,n}, then x* is unstable and lies at the Julia set.

gi(x),i=1,2,,nbeing the coordinate functions of the fixed point multivariate function G.

Let us remark that if the order of convergence of the iterative method is at least two, then the roots of the nonlinear function are superattracting fixed points of the vectorial rational function corresponding to the iterative scheme. If a fixed point is not a root of the nonlinear function, it is called strange fixed point and its character can be analyzed in the same manner.

The concept of critical point can be defined following the idea of multivariate convergence of iterative methods.

Definition 3.1. A fixed point x* is a critical point of G if its coordinate functions gi(x)satisfy, for all i,j ϵ{1,2,,n}.

From this definition, a superattracting fixed point will also be a critical point of the operator and, from numerical point of view, the iterative scheme will have, at least, quadratic order of convergence. A critical point that is not root of p(x) will be called free critical point.

The stability of this family is studied on two systems of real quadratic polynomials, p(x) = 0 and q(x) = 0, where

p(x1,x2)={x121,x221,andq(x1,x2)={x121,x221.

The components of the iteration function G are, in case of p(x),

g1(x1,x2)=8θ2x12(x121)3+θ4(x121)416x14(3x14+6x121)128x17,
g2(x1,x2)=8θ2x22(x221)3+θ4(x221)416x24(3x24+6x221)128x27.

It was proved by Cordero et al. [33] that the number of fixed points of the vectorial rational function G(x) associated with HMT iterative method on p(x) is 64, four of them (those corresponding to the roots of p(x)) are superattractive for any value of θ. Moreover, there are 48 unstable strange fixed points and the character of the final 12 real strange fixed points is simultaneously attractive for two ranges of values of θ, [−0.3847551, −0.3838109] and [0.3838109, 0.3847551], being superattractive if θ0.3838109or θ0.3838109.

In Figure 4, we can see the dynamical plane of HMT method for θ0.3838109, where those 12 attracting strange fixed points appear as white circles, with small basins of attraction. The roots of p(x) appear as white stars, in their own basins of attraction.

Figure 4.

Dynamical plane of multivariate HMT method on p(x) for θ ≈ −0.3838109.

The parametric plots, as the extension of parameter planes in multivariable case, were revealed to be an interesting procedure that allowed to detect the most stable and unstable elements of a family of iterative methods. By using the information that gives us the iteration of the elements of the family on the different free critical points, we can know the global behavior of the family depending on the value of the parameter. Orbits of each free critical point are showed in Figure 5. In each one of these pictures, a different free critical point is used as initial guess of each member of the class of iterative schemes, taking values of parameter θ in [−5,5]. The values of θ corresponding to the members of the family are placed in the abscissa axis and the ordinate axis corresponds to 0.1, 0.2, 0.3, or 0.4 if the iterative process has converged to each one of the solutions of the quadratic polynomial real system, p(x), respectively. Moreover, the ordinate of a point is −0.1 if the process diverges and it is null in other cases.

Figure 5.

Different parameter plots of HMT family on p(x).

Figure 6.

4-Periodic orbits in HMT family.

The goal of these graphics is to show the elements of the family HMT (that is, the values of parameter θ) that present unstable behavior (attracting strange fixed points, periodic orbits, …). These plots were obtained by using 20,000 subintervals, a maximum of 40 iterations and an error estimation of 10−3, when iterates tend to a fixed point.

Thanks to this analysis, some non-desirable behavior was detected, as attracting periodic orbits of different periods, as can be observed in Figure 6, where two periodic orbits of period four are showed in yellow.

4. Numerical results

In this section, we are going to check the numerical performance of some elements of family (Eq. 10) compared with the multidimensional version of some other known methods of the same order and also with Newton's one.

Definition 4.1. Let ξ be a zero of function F and suppose that x(k−1), x(k) and x(k+1) are three consecutive iterations close to ξ. Then, the computational order of convergence p can be approximated using the formula (see Ref. [12]).

pln(x(k+1)x(k)/x(k)x(k1))ln(x(k)x(k1)/x(k1)x(k2))

Numerical computations have been carried out in MATLAB, with variable precision arithmetic that uses floating point representation of 1000 decimal digits of mantissa. If n > 1, every iterate x(k+1) is obtained from the previous one, x(k), by adding one term of the form A−1b. Matrix A and vector b are different according to the method used, but in any case the inverse calculation −A−1b is computed by solving the linear system Ay = −b, by using Gaussian elimination with partial pivoting. Nevertheless, if several systems with the same matrix A must be solved in iteration, LU factorization is made once and the corresponding triangular systems are solved by substitution. The stopping criterion used is F(x(k))<10700or x(k+1)x(k)<10700.

4.1. Molecular interaction problem

We are going to solve the equation of molecular interaction,

uxx+uyy=u2,(x,y)[0,1]×[0,1],

with the boundary conditions u(x,0)=2x2x+1,u(x,1)=2,u(0,y)=2y2y+1and u(1, y) = 2, for all (x, y) in their domain.

For approximating its solution, we are going to use central divided differences procedure to transform the problem in a nonlinear system of equations. This system is solved by means of the proposed methods of order four and another known ones.

The discretization process yields the nonlinear system of equations,

ui+1,j4ui,j+ui1,j+ui,j+1+ui,j1h2ui,j2=0,i=1,2,,nx,j=1,2,,ny,

where ui,j denotes the estimation of the unknown u(xi,yj), where xi = ih with i = 1,2,…, nx and yj = jk with j = 1,2, …, ny are the nodes in both variables, with h=1nx,k=1nyand nx = ny.

For checking purposes, we consider nx = ny = 4, getting a mesh of 5 × 5 points. Applying the boundary conditions we have only nine unknowns, which we rename as:

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.E4563

Then, the nonlinear system associated with the partial differential equation, including the boundary conditions, can be expressed as

F(x)=Ax+φ(x)b=0,E4564

where

A=(MI0IMI0IM),M=(410141014)E4565

I being the 3 × 3 identity matrix, φ(x)=h2(x12,x22,,x92)Tand b=(74,1,278,1,0,2,278,2,4)T.

In this case,

F'(x)=A+2h2diag(x1,x2,x9).E4566

For solving this system, we apply the extension for systems of Ostrowski's method (OM), Chun's scheme (CM), Jarratt's method (JM), Newton's method (NM), and three elements of family (Eq. 10) obtained by selecting stable values of the parameters.

MA:a1=54,b2=0,MB:a1=1,b2=1,MC:a1=1,b2=3.E4567

In Table 4, the approximated computational order of convergence ACOC, the number of iterations, the difference between the two last iterations and the residual of the function at the last iteration are shown, for each one of the methods. In all cases, the initial estimation is a null vector and the Euclidean norm is used in the calculation of the residuals.

MethodACOCIterationx(k+1)x(k)F(x(k))
NM1.999991.482e-4136.448e-828
JM3.995451.482e-4131.976e-1007
OM3.996451.482e-4131.618e-1007
CM3.995951.998e-3531.618e-1007
MA4.051955.362e-5101.707e-2007
MB3.996057.123e-3621.409e-1449
MC3.996053.110e-3623.811e-1451

Table 4.

Numerical results for molecular interaction problem.

All the checked schemes provide the same solution of the nonlinear system. In Table 4, we can observe that all the fourth-order methods have a similar performance, but we note that the lowest error corresponds to method MA, duplicating the number of exact digits with respect to the other ones.

5. Conclusions

Many problems in science and engineering are modeled in such a way that, for their solution, it is necessary to solve systems of nonlinear equations. Therefore, designing iterative methods for solving these types of problems is an important task and it is a fruitful area of research. In this chapter, a review of the different techniques for constructing iterative methods is presented. Moreover, it is shown that real discrete dynamics tools are useful for analyzing the stability of the designed methods, selecting those with good dynamical behavior. In the numerical section, a chemical problem is used for testing the presented methods and the theoretical results are confirmed.

Acknowledgments

This research was partially supported by Ministerio de Economía y Competitividad of Spain, MTM2014-52016-C2-2-P, and by Ministry of Higher Education Science and Technology of Dominican Republic, FONDOCYT 2014-1C1-088.

© 2016 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution 3.0 License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

How to cite and reference

Link to this chapter Copy to clipboard

Cite this chapter Copy to clipboard

Alicia Cordero, Juan R. Torregrosa and Maria P. Vassileva (October 19th 2016). Design, Analysis, and Applications of Iterative Methods for Solving Nonlinear Systems, Nonlinear Systems - Design, Analysis, Estimation and Control, Dongbin Lee, Tim Burg and Christos Volos, IntechOpen, DOI: 10.5772/64106. Available from:

chapter statistics

920total chapter downloads

More statistics for editors and authors

Login to your personal dashboard for more detailed statistics on your publications.

Access personal reporting

Related Content

This Book

Next chapter

Nonlinear State and Parameter Estimation Using Iterated Sigma Point Kalman Filter: Comparative Studies

By Marwa Chaabane, Imen Baklouti, Majdi Mansouri, Nouha Jaoua, Hazem Nounou, Mohamed Nounou, Ahmed Ben Hamida and Marie‐France Destain

Related Book

First chapter

Recent Advances in Fragment Molecular Orbital-Based Molecular Dynamics (FMO-MD) Simulations

By Yuto Komeiji, Yuji Mochizuki, Tatsuya Nakano and Hirotoshi Mori

We are IntechOpen, the world's leading publisher of Open Access books. Built by scientists, for scientists. Our readership spans scientists, professors, researchers, librarians, and students, as well as business professionals. We share our knowledge and peer-reveiwed research papers with libraries, scientific and engineering societies, and also work with corporate R&D departments and government entities.

More About Us