Fourth- and Fifth-Order Methods for Solving Nonlinear Systems of Equations: An Application to the Global Positioning System

Abstract

Two iterative methods of order four and five, respectively, are presented for solving nonlinear systems of equations. Numerical comparisons are made with other existing second- and fourth-order schemes to solve the nonlinear system of equations of the Global Positioning System and some academic nonlinear systems.

1. Introduction

The search for solutions of nonlinear systems of equations is an old and difficult problem with wide applications in sciences and engineering. The best known method, for being very simple and effective, is Newton’s method. Its generalization to a system of equations was proposed by Ostrowski [1] and to Banach spaces by Kantorovič [2]. In the literature, several modifications have been made on classical methods in order to accelerate the convergence or to reduce the number of operations and evaluations of functions in each step of the iterative process. The extension of the variants of Newton’s method described by Weerakoon and Fernando in [3], by Özban in [4] and Gerlach in [5], to the functions of several variables has been developed in [6–9]. In [6, 7], families of variants of Newton’s method of third-order have been designed by using open and closed formulas of quadrature, including the families of the methods defined by Frontini and Sormani in [8]. Using the generic formula of the interpolatory quadrature, in [9] a family of methods is obtained with order of convergence 2𝑑+1, under certain conditions, where 𝑑 is the order up to which the partial derivatives of each coordinate function evaluated in the solution are canceled. Indeed, Darvishi and Barati improved in [10] the method from Frontini and Sormani, getting a fourth-order scheme. In addition to multistep methods based on interpolatory quadrature, other schemes have been developed by using different techniques, as extension to several variables of one-dimensional schemes (see [11]), Adomian decomposition (see [12, 13], e.g.), the one proposed by Darvishi and Barati in [14, 15] with super cubic convergence, and the methods proposed by Cordero et al. in [16] with orders of convergence four and five. Another procedure to develop iterative methods for nonlinear systems is the replacement of the second derivative by some approximation. In [17], Traub presented a family of multipoint methods based on approximating the second derivative that appears in the iterative formula of Chebyshev’s scheme, and more recently, Babajee et al. in [18] designed two Chebyshev-like methods free from second derivatives. Recently, Sharma et al. [19] designed a fourth-order scheme by using weight-function technique. Another well-known acceleration technique is the composition of two iterative methods of orders 𝑝₁ and 𝑝₂, respectively, obtaining a method of order 𝑝₁𝑝₂ (see [17]). New evaluations of the Jacobian matrix and the nonlinear function are usually needed in order to increase the order of convergence.

Now, we are going to introduce the problem and some necessary concepts in order to develop the modified methods and to analyze their convergence. Let us consider the problem of finding a real zero of a function 𝐹∶𝐷⊆ℝ^𝑛→ℝ^𝑛, that is, a solution 𝑥∈𝐷 of the nonlinear system 𝐹(𝑥)=0 of 𝑛 equations with 𝑛 unknowns. The best known iterative method is the classical Newton method given by 𝑥^(𝑘+1)=𝑥^(𝑘)−𝐹^′𝑥^(𝑘)−1𝐹𝑥^(𝑘),𝑘=0,1,…, (1) where 𝐹^′(𝑥^(𝑘)) is the Jacobian matrix of the function 𝐹 evaluated in the 𝑘th iteration 𝑥^(𝑘). Traub, in [17], introduced a variant of Newton’s method of convergence order three. We are going to describe it because our methods combine Traub with Newton's method. Traub's scheme consists of the composition of Newton's method with itself, but with a frozen Jacobian matrix, its iterative expression is 𝑥^(𝑘+1)=𝑥^(𝑘)−𝐹^′𝑥^(𝑘)−1𝐹𝑥^(𝑘)+𝐹𝑦^(𝑘),𝑘=0,1,…, (2) where 𝑦^(𝑘) is the 𝑘th iteration of Newton’s method.

On the other hand, recently Sharma et al. in [19] have developed a fourth-order method for solving nonlinear systems of equations. The algorithm is composed of two weighted Newton steps, and it is given by 𝑦^(𝑘)=𝑥^(𝑘)−23𝐹^′𝑥^(𝑘)−1𝐹𝑥^(𝑘),𝑥^(𝑘+1)=𝑥^(𝑘)−129−𝐼+4𝐹^′𝑦^(𝑘)−1𝐹^′𝑥^(𝑘)+34𝐹^′𝑥^(𝑘)−1𝐹^′𝑦^(𝑘)×𝐹^′𝑥^(𝑘)−1𝐹𝑥^(𝑘),𝑘=0,1,…. (3)

In the following, we remember some known notions and results that we need in order to analyze the convergence of the new methods. Let 𝐹∶𝐷⊆ℝ^𝑛→ℝ^𝑛 be sufficiently Frechet differentiable in 𝐷. By using the notation introduced in [20], the 𝑞th derivative of 𝐹 at 𝑢∈ℝ^𝑛, 𝑞≥1 is the 𝑞-linear function 𝐹^(𝑞)(𝑢)∶ℝ^𝑛×ℝ^𝑛×⋯×ℝ^𝑛→ℝ^𝑛 such that 𝐹^(𝑞)(𝑢)(𝑣₁,𝑣₂,…,𝑣_𝑞)∈ℝ^𝑛. It is easy to observe that(1)𝐹^(𝑞)(𝑢)(𝑣₁,𝑣₂,…,𝑣_𝑞−1,⋅)∈ℒ(ℝ^𝑛), (2)𝐹^(𝑞)(𝑢)(𝑣_𝜎(1),𝑣_𝜎(2),…,𝑣_𝜎(𝑞))=𝐹^(𝑞)(𝑢)(𝑣₁,𝑣₂,…,𝑣_𝑞) for all permutation 𝜎 of {1,2,…,𝑞},where ℒ(ℝ^𝑛) is the set of lineal operators of ℝ^𝑛 in ℝ^𝑛.

From the above properties, we can use the following notation:(1)𝐹^(𝑞)(𝑢)(𝑣₁,𝑣₂,…,𝑣_𝑞)=𝐹^(𝑞)(𝑢)𝑣₁⋅⋯⋅𝑣_𝑞, (2)𝐹^(𝑞)(𝑢)𝑣^𝑞−1𝐹^(𝑝)𝑣^𝑝=𝐹^(𝑞)(𝑢)𝐹^(𝑝)(𝑢)𝑣^{𝑞+𝑝−1}.

On the other hand, for 𝑥+ℎ∈ℝ^𝑛 lying in a neighborhood of a solution 𝑥 of 𝐹(𝑥)=0, we can apply Taylor's expansion, and assuming that the Jacobian matrix 𝐹^′(𝑥) is nonsingular, we have 𝐹(𝑥+ℎ)=𝐹^′(𝑥)ℎ+^𝑝−1_𝑞=2𝐶_𝑞ℎ^𝑞+𝑂ℎ^𝑝, (4) where 𝐶_𝑞=(1/𝑞!)[𝐹^′(𝑥)]⁻¹𝐹^(𝑞)(𝑥), 𝑞≥2. We observe that 𝐶_𝑞ℎ^𝑞∈ℝ^𝑛 since 𝐹^(𝑞)(𝑥)∈ℒ(ℝ×⋯×ℝ,ℝ) and [𝐹^′(𝑥)]⁻¹∈ℒ(ℝ^𝑛).

In addition, we can express 𝐹^′ as 𝐹^′(𝑥+ℎ)=𝐹^′(𝑥)𝐼+^𝑝−1_𝑞=2𝑞𝐶_𝑞ℎ^𝑞−1+𝑂ℎ^𝑝=𝐹^′(𝑥)𝐷(ℎ)+𝑂ℎ^𝑝, (5) where 𝐼 is the identity matrix. Therefore, 𝑞𝐶_𝑞ℎ^𝑞−1∈ℒ(ℝ^𝑛). From the previous equation, we obtain 𝐹^′(𝑥+ℎ)⁻¹=𝐷(ℎ)⁻¹𝐹^′(𝑥)⁻¹+𝑂ℎ^𝑝. (6)

Then, if the inverse of 𝐷(ℎ) is 𝐷(ℎ)⁻¹=𝐼+𝑋₂ℎ+𝑋₃ℎ²+𝑋₄ℎ³+⋯, (7)provided that 𝑋_𝑖, 𝑖=2,3,… verify 𝐷(ℎ)𝐷(ℎ)⁻¹=𝐷(ℎ)⁻¹𝐷(ℎ)=𝐼. (8)

Solving the system involved in (8), we have that 𝑋₂=−2𝐶₂,𝑋₃=4𝐶²₂−3𝐶₃,𝑋₄=−8𝐶³₂+6𝐶₂𝐶₃+6𝐶₃𝐶₂−4𝐶₄,⋮ (9)

We denote 𝑒^(𝑘)=𝑥^(𝑘)−𝑥 as the error in the 𝑘th iteration. The equation𝑒^(𝑘+1)=𝐿𝑒^(𝑘)𝑝+𝑂𝑒^(𝑘)𝑝+1, (10) where 𝐿 is a 𝑝-linear function 𝐿∈ℒ(ℝ×⋯×ℝ,ℝ), is called the error equation, and 𝑝 is the order of convergence. Observe that 𝑒^(𝑘)𝑝 is (𝑒^(𝑘),𝑒^(𝑘),…,𝑒^(𝑘)).

In [7], the concept of computational order of convergence was introduced as follows.

Definition 1. Let 𝑥 be a zero of a function 𝐹, and suppose that 𝑥^(𝑘−1), 𝑥^(𝑘), and 𝑥^(𝑘+1) are three consecutive iterations close to 𝑥. Then, the computational order of convergence 𝑝 can be approximated using the formula 𝑥𝑝≈ACOC=ln^(𝑘+1)−𝑥^(𝑘)/𝑥^(𝑘)−𝑥^(𝑘−1)ln𝑥^(𝑘)−𝑥^(𝑘−1)/𝑥^(𝑘−1)−𝑥^(𝑘−2). (11)

In addition, in order to compare different methods, we use the efficiency index, 𝐼=𝑝^1/𝑑, where 𝑝 is the order of convergence and 𝑑 is the total number of new functional evaluations (per iteration) required by the method. This is the most used index, but not the only one. In [17], Traub uses a computational index defined as 𝐶𝐼=𝑝^1/𝑜𝑝, where 𝑜𝑝 is the number of products/quotients per iteration. We recall that the number of products and quotients that we need for solving 𝑚 linear systems with the same matrix of coefficients, by using 𝐿𝑈 factorization, is 13𝑛³+𝑚𝑛²−13𝑛, (12) where 𝑛 is the size of each system. We will use these indices in order to compare the different iterative methods. Kung and Traub in [21] conjectured that the order of convergence of any multipoint method without memory for solving nonlinear equations cannot exceed the bound 2^𝑑−1 (called the optimal order). Ostrowski's method [1], Jarrett's scheme [22], and King's procedure [23] are some of the optimal one-dimensional fourth-order methods. We have adapted the definition of optimal order of convergence to the case of iterative methods to solve nonlinear systems. The extension to several variables of the conjecture of Kung and Traub could be done in the following way [24].

Conjecture 2. Given a multipoint iterative method to solve nonlinear systems of equations which requires 𝑑=𝑘₁+𝑘₂ functional evaluations per step such that 𝑘₁ of them correspond to the functional evaluations of the Jacobian matrix and 𝑘₂ to evaluations of the nonlinear function. We conjecture that the optimal order for this method is 2^{𝑘1+𝑘2−1} if 𝑘₁≤𝑘₂.

In this paper, we propose two new and competitive iterative methods of orders four and five, respectively, that improve other known methods.

The rest of this paper is organized as follows: in Section 2, we make an introduction to the Global Positioning System (GPS), focusing on the way that the receiver calculates the user position using the ephemeris data of the artificial satellites. In Section 3, we present our new iterative methods and analyze its convergence order, and by using the idea of a technique presented in [25], it is also proved that, in general, if we combine two methods of orders 𝑝 and 𝑞, respectively, with 𝑝≥𝑞, in the same way that we do it in our method of order five, the order of convergence of the resultant method is 𝑝+𝑞. In Section 4, we show an application of this analysis in order to solve the nonlinear system of the GPS and several academic nonlinear systems of equations. A comparison is established among the new methods and Newton and Sharma’s methods in terms of convergence order, approximated computational convergence order (ACOC), and computational and efficiency indices, 𝐶𝐼 and 𝐼, respectively.

2. Basics on Global Positioning System

This section introduces the basic concept of how a GPS receiver determines its position. From the satellite constellation, the equations required for solving the user position conform a nonlinear system of equations. In addition, some practical considerations (i.e., the inaccuracy of the user clock) will be included in these equations. These equations are usually solved through a linearization and a fixed point iteration method. The obtained solution is in a Cartesian coordinate system, and after that the result will be converted into a spherical coordinate system. However, the Earth is not a perfect sphere; therefore, once the user position is estimated, the shape of the Earth must be taken into consideration. The user position is then translated into the Earth-based coordinate system. In this paper, we are going to focus our attention in solving the nonlinear system of equations of the GPS giving the results in a Cartesian coordinate system. We can find further information about GPS in [26].

2.1. Basic GPS Concepts

The position of a point in space can be found by using the distances measured from this point to some known position in space. We are going to use an example to illustrate this point.

Figure 1 shows a two-dimensional case. In order to determine the user position 𝑈, three satellites 𝑆₁, 𝑆₂, and 𝑆₃ and three distances are required. The trace of a point with constant distance to a fixed point is a circle in the two-dimensional case. Two satellites and two distances give two possible solutions because two circles intersect at two points. A third circle is needed to uniquely determine the user position. For similar reasons in a three-dimensional case, four satellites and four distances are needed. The equal-distance trace to a fixed point is a sphere in a three-dimensional case. Two spheres intersect to make a circle. This circle intersects another sphere, and this intersection produces two points. In order to determine which point is the user position, one more satellite should be needed. In GPS, the position of the satellite is known from the ephemeris data transmitted by the satellite. By measuring the distance from the receiver to the satellite, the position of the receiver can be determined. In the above discussion, the distance measured from the user to the satellite is assumed to be very accurate, and there is no bias error. However, the distance measured between the receiver and the satellite has a constant unknown bias, because the user clock usually is different from the GPS clock. In order to solve this bias error, one more satellite is required. Therefore, in order to find the user position, five satellites are needed. If one uses four satellites and the measured distance with bias error to measure a user position, two possible solutions can be obtained. Theoretically, one cannot determine the user position. However, one of the solutions is close to the Earth's surface, and the other one is in the space. In fact, as we will see in Section 4, in this memory, we have used four satellites, and sometimes we have found the solution in the space. Since the user position is usually close to the surface of the earth, it can be uniquely determined. Therefore, the general statement is that four satellites can be used to determine a user position, even though the distance measured has a bias error. The method of solving the user position discussed in the next subsections is through iteration. The initial position is often selected at the center of the Earth. In the following discussion, four satellites are considered as the minimum number required for finding the user position.

Figure 1 

Two-dimensional user position.

2.2. Basic Equations for Finding User Position

In this section, the basic equations for determining the user position will be presented. Assume that the distance measured is accurate, and under this condition, three satellites should be sufficient. Let us suppose that there are three known points at locations 𝑟₁ or (𝑥₁,𝑦₁,𝑧₁), 𝑟₂ or (𝑥₂,𝑦₂,𝑧₂), and 𝑟₃ or (𝑥₃,𝑦₃,𝑧₃) and an unknown point at 𝑟_𝑢 or (𝑥_𝑢,𝑦_𝑢,𝑧_𝑢). If the distances between the three known points to the unknown point can be measured as 𝜌₁, 𝜌₂, and 𝜌₃, these distances can be written as 𝜌₁=𝑥₁−𝑥_𝑢²+𝑦₁−𝑦_𝑢²+𝑧₁−𝑧_𝑢²,𝜌₂=𝑥₂−𝑥_𝑢²+𝑦₂−𝑦_𝑢²+𝑧₂−𝑧_𝑢²,𝜌₃=𝑥₃−𝑥_𝑢²+𝑦₃−𝑦_𝑢²+𝑧₃−𝑧_𝑢². (13)

Because there are three unknowns and three equations, the values of 𝑥_𝑢, 𝑦_𝑢, and 𝑧_𝑢 can be determined from these equations. Theoretically, there should be two sets of solutions as they are second-order equations. These equations can be solved by linearizing them and making an iterative approach. The solution of these equations will be discussed later in Section 2.4. In GPS operation, the positions of the satellites are given. This information can be obtained from the data transmitted from the satellites. The distances from the user (the unknown position) to the satellites must be measured simultaneously at a certain time instance. Each satellite transmits a signal with a time reference associated with it. By measuring the time of the signal traveling from the satellite to the user, the distance between the user and the satellite can be found. The distance measurement is discussed in the next section.

2.3. Measurement of Pseudorange

Every satellite sends a signal at a certain time 𝑡_𝑠𝑖. The receiver will receive the signal at a later time 𝑡_𝑢. The distance between the user and the satellite 𝑖 can be determined as 𝜌_𝑖𝑇=𝑐𝑡_𝑢−𝑡_𝑠𝑖, (14) where 𝑐 is the speed of light, 𝜌_𝑖𝑇 is often referred to as the true value of pseudorange from user to satellite 𝑖, 𝑡_𝑠𝑖 is referred to as the true time of transmission from satellite 𝑖, and 𝑡_𝑢 is the true time of reception. From a practical point of view, it is difficult, if not impossible, to obtain the correct time from the satellite or the user. The actual satellite clock time 𝑡^′_𝑠𝑖 and actual user clock time 𝑡^′_𝑢 are related to the true time as 𝑡^′_𝑠𝑖=𝑡_𝑠𝑖+Δ𝑏_𝑖,𝑡^′_𝑢=𝑡_𝑢+𝑏_𝑢𝑡, (15) where Δ𝑏_𝑖 is the satellite clock error and 𝑏_𝑢𝑡 is the user clock bias error. Besides the clock error, there are other factors affecting the pseudorange measurement. The measured pseudorange 𝜌_𝑖 can be written as 𝜌_𝑖=𝜌_𝑖𝑇+Δ𝐷_𝑖−𝑐Δ𝑏_𝑖−𝑏_𝑢𝑡+𝑐Δ𝑇_𝑖+Δ𝐼_𝑖+𝑣_𝑖+Δ𝑣_𝑖, (16) where Δ𝐷_𝑖 is the satellite position error effect on the range, Δ𝑇_𝑖 is the tropospheric delay error, Δ𝐼_𝑖 is the ionospheric delay error, 𝑣_𝑖 is the receiver measurement noise error, and Δ𝑣_𝑖 is the relativistic time correction. Some of these errors can be corrected; for example, the tropospheric delay can be modeled, and the ionospheric error can be corrected in a two-frequency receiver. The errors will cause inaccuracy of the user position. However, the user clock error cannot be corrected through receiver information. Thus, it will remain as an unknown. So, the system of (13) must be modified as 𝜌₁=𝑥₁−𝑥_𝑢²+𝑦₁−𝑦_𝑢²+𝑧₁−𝑧_𝑢²+𝑏_𝑢,𝜌₂=𝑥₂−𝑥_𝑢²+𝑦₂−𝑦_𝑢²+𝑧₂−𝑧_𝑢²+𝑏_𝑢,𝜌₃=𝑥₃−𝑥_𝑢²+𝑦₃−𝑦_𝑢²+𝑧₃−𝑧_𝑢²+𝑏_𝑢,𝜌₄=𝑥₄−𝑥_𝑢²+𝑦₄−𝑦_𝑢²+𝑧₄−𝑧_𝑢²+𝑏_𝑢, (17) where 𝑏_𝑢 is the user clock bias error expressed in distance, which is related to the quantity 𝑏_𝑢𝑡 by 𝑏_𝑢=𝑐𝑏_𝑢𝑡. In the system of (17), four equations are needed to solve four unknowns 𝑥_𝑢, 𝑦_𝑢, 𝑧_𝑢, and 𝑏_𝑢. Thus, in a GPS receiver, a minimum of four satellites is required to solve the user position.

2.4. Solution of User Position from Pseudoranges

One common way to solve the system of (17) is to linearize them. The system can be written in a simplified form as 𝜌_𝑖=𝑥_𝑖−𝑥_𝑢²+𝑦_𝑖−𝑦_𝑢²+𝑧_𝑖−𝑧_𝑢²+𝑏_𝑢, (18) (with 𝑖=1,2,3,4, and 𝑥_𝑢, 𝑦_𝑢, 𝑧_𝑢 and 𝑏_𝑢) are the unknowns. The pseudorange 𝜌_𝑖 and the positions of the satellites 𝑥_𝑖, 𝑦_𝑖, 𝑧_𝑖 are known. By differentiating (18), 𝛿𝜌_𝑖=𝑥_𝑖−𝑥_𝑢𝛿𝑥_𝑢+𝑦_𝑖−𝑦_𝑢𝛿𝑦_𝑢+𝑧_𝑖−𝑧_𝑢𝛿𝑧_𝑢𝑥_𝑖−𝑥_𝑢²+𝑦_𝑖−𝑦_𝑢²+𝑧_𝑖−𝑧_𝑢²+𝛿𝑏_𝑢=𝑥_𝑖−𝑥_𝑢𝛿𝑥_𝑢+𝑦_𝑖−𝑦_𝑢𝛿𝑦_𝑢+𝑧_𝑖−𝑧_𝑢𝛿𝑧_𝑢𝜌_𝑖−𝑏_𝑢+𝛿𝑏_𝑢. (19)

In (19), 𝛿𝑥_𝑢, 𝛿𝑦_𝑢, 𝛿𝑧_𝑢, and 𝛿𝑏_𝑢 can be considered as the only unknowns. The quantities 𝑥_𝑢, 𝑦_𝑢, 𝑧_𝑢, and 𝑏_𝑢 are treated as known values because one can assume some initial values for these quantities. From these initial values, a new set of 𝛿𝑥_𝑢, 𝛿𝑦_𝑢, 𝛿𝑧_𝑢, and 𝛿𝑏_𝑢 can be calculated. These values are used to modify the original 𝑥_𝑢, 𝑦_𝑢, 𝑧_𝑢, and 𝑏_𝑢 to find another new set of solutions. This new set of 𝑥_𝑢, 𝑦_𝑢, 𝑧_𝑢, and 𝑏_𝑢 can be considered again as known quantities. This process continues until the absolute values of 𝛿𝑥_𝑢, 𝛿𝑦_𝑢, 𝛿𝑧_𝑢, and 𝛿𝑏_𝑢 are very small and within a certain predetermined limit. The final values of 𝑥_𝑢, 𝑦_𝑢, 𝑧_𝑢, and 𝑏_𝑢 are the desired solution. This method is often referred to as an iteration method of fixed point. With 𝛿𝑥_𝑢, 𝛿𝑦_𝑢, 𝛿𝑧_𝑢, and 𝛿𝑏_𝑢 as unknowns, the above equation becomes a set of linear equations. This procedure is often referred to as linearization. The expression (19) can be written in matrix form as 𝛿𝜌₁𝛿𝜌₂𝛿𝜌₃𝛿𝜌₄=𝛼₁₁𝛼₁₂𝛼₁₃1𝛼₂₁𝛼₂₂𝛼₂₃1𝛼₃₁𝛼₃₂𝛼₃₃1𝛼₄₁𝛼₄₂𝛼₄₃1𝛿𝑥_𝑢𝛿𝑦_𝑢𝛿𝑧_𝑢𝛿𝑏_𝑢, (20) where 𝛼_𝑖1=𝑥_𝑖−𝑥_𝑢𝜌_𝑖−𝑏_𝑢,𝛼_𝑖2=𝑦_𝑖−𝑦_𝑢𝜌_𝑖−𝑏_𝑢,𝛼_𝑖3=𝑧_𝑖−𝑧_𝑢𝜌_𝑖−𝑏_𝑢. (21)

The solution of (20) is 𝛿𝑥_𝑢𝛿𝑦_𝑢𝛿𝑧_𝑢𝛿𝑏_𝑢=𝛼₁₁𝛼₁₂𝛼₁₃1𝛼₂₁𝛼₂₂𝛼₂₃1𝛼₃₁𝛼₃₂𝛼₃₃1𝛼₄₁𝛼₄₂𝛼₄₃1⁻¹𝛿𝜌₁𝛿𝜌₂𝛿𝜌₃𝛿𝜌₄. (22)

This process obviously does not provide the needed solutions directly. However, the desired solutions can be obtained from it. In order to find the desired position solution, this procedure must be used repetitively in an iterative way. A quantity is often used to determine whether the desired result is reached, and this quantity can be defined as 𝛿𝜐=𝛿𝑥²_𝑢+𝛿𝑦²_𝑢+𝛿𝑧²_𝑢+𝛿𝑏²_𝑢. (23)

When 𝛿𝜐 is lower than a certain predetermined threshold, the iteration will stop. Sometimes, the clock bias 𝑏_𝑢 is not included in (23). In this paper, we use as stopping criterion the quantity ||𝑥^(𝑘+1)−𝑥^(𝑘)||+||𝐹(𝑥^(𝑘+1))|| because it is stronger than (23). As we can verify in [27] the above iterative method used to calculate via software, the receiver position in the GPS is Newton's method, a well-known method of second-order of convergence. In this work, we improve the GPS software by means of two methods of order four and five, respectively, that converge to the solution with less number of iterations and better 𝐼 or 𝐶𝐼 than Newton scheme.

3. Description of the Methods and Convergence Analysis

3.1. A Fourth-Order Method

In this section, we display a new method for solving nonlinear systems that we call 𝑀4, obtained by combining (Newton and Traub's method). Its iterative expression is 𝑧^(𝑘)=𝑥^(𝑘)−𝐹^′𝑥^(𝑘)−1𝐹𝑥^(𝑘)+𝐹𝑦^(𝑘),𝑥^(𝑘+1)=𝑦^(𝑘)−𝐹^′𝑧^(𝑘)−1𝐹𝑦^(𝑘), (24) where 𝑦^(𝑘) is the 𝑘th iteration of Newton’s method. In the next result, we are going to prove that the convergence order of the method is four.

Theorem 3. Let 𝐹∶𝐷⊆ℝ^𝑛→ℝ^𝑛 be sufficiently differentiable at each point of an open neighborhood 𝐷 of  𝑥∈ℝ^𝑛 that is a solution of the nonlinear system 𝐹(𝑥)=0. Let one suppose that 𝐹^′(𝑥) is continuous and nonsingular in 𝑥. Then, the sequence {𝑥^(𝑘)}_𝑘≥0 obtained by using the iterative expression (24) converges to 𝑥 with order four. The error equation is 𝑒^(𝑘+1)=−𝐶³₂𝑒^(𝑘)4+𝑂𝑒^(𝑘)5, (25) where 𝐶_𝑘=(1/𝑘!)[𝐹^′(𝑥)]⁻¹𝐹^(𝑘)(𝑥),𝑘=2,3,…, and 𝑒^(𝑘)=𝑥^(𝑘)−𝑥.

Proof. Taylor expansion of 𝐹(𝑥^(𝑘)) and 𝐹^′(𝑥^(𝑘)) around 𝑥 gives 𝐹𝑥^(𝑘)=𝐹^′(𝑥)𝑒^(𝑘)+𝐶₂𝑒^(𝑘)2+𝐶₃𝑒^(𝑘)3+𝐶₄𝑒^(𝑘)4+𝑂𝑒^(𝑘)5,𝐹^′𝑥^(𝑘)=𝐹^′(𝑥)𝐼+2𝐶₂𝑒^(𝑘)+3𝐶₃𝑒^(𝑘)2+4𝐶₄𝑒^(𝑘)3+5𝐶₅𝑒^(𝑘)4+𝑂𝑒^(𝑘)5, (26) where 𝐶_𝑘=(1/𝑘!)[𝐹^′(𝑥)]⁻¹𝐹^(𝑘)(𝑥),𝑘=2,3,…, and 𝑒^(𝑘)=𝑥^(𝑘)−𝑥. From (26), we obtain 𝐹^′𝑥^(𝑘)−1=𝐼+𝑋₂𝑒^(𝑘)+𝑋₃𝑒^(𝑘)2+𝑋₄𝑒^(𝑘)3×𝐹^′(𝑥)⁻¹+𝑂𝑒^(𝑘)4, (27) where 𝑋₂=−2𝐶₂, 𝑋₃=4𝐶²₂−3𝐶₃, and 𝑋₄=−8𝐶³₂+6𝐶₂𝐶₃+6𝐶₃𝐶₂−4𝐶₄. Taylor's expansion of 𝑦^(𝑘) is 𝑦^(𝑘)−𝑥=𝐶₂𝑒^(𝑘)2+2𝐶₃−2𝐶²₂𝑒^(𝑘)3+4𝐶³₂−4𝐶₂𝐶₃−3𝐶₃𝐶₂+3𝐶₄𝑒^(𝑘)4+𝑂𝑒^(𝑘)5. (28)
On the other hand, we have that 𝐹𝑦^(𝑘)=𝐹^′(𝑥)𝑦^(𝑘)−𝑥+𝐶₂𝑦^(𝑘)−𝑥²+𝑂𝑒^(𝑘)5, (29) and by operating, we get 𝐹𝑦^(𝑘)=𝐹^′(𝑥)𝐶₂𝑒^(𝑘)2+2𝐶₃−2𝐶²₂𝑒^(𝑘)3+5𝐶³₂−4𝐶₂𝐶₃−3𝐶₃𝐶₂+3𝐶₄𝑒^(𝑘)4+𝑂𝑒^(𝑘)5. (30)
Analogously, we obtain the expression of 𝑧^(𝑘)−𝑥. Given that the 𝑘th iteration of Traub's scheme is 𝑧^(𝑘)=𝑦^(𝑘)−𝐹^′𝑥^(𝑘)−1𝐹𝑦^(𝑘), (31) then 𝑧^(𝑘)−𝑥=2𝐶²₂𝑒^(𝑘)3+−9𝐶³₂+4𝐶₂𝐶₃+3𝐶₃𝐶₂𝑒^(𝑘)4+𝑂𝑒^(𝑘)5. (32)
Besides, the expression of 𝐹^′(𝑧^(𝑘)) is 𝐹^′𝑧^(𝑘)=𝐹^′(𝑥)𝐼+2𝐶₂𝑧^(𝑘)−𝑥+𝑂𝑒^(𝑘)5, (33) or, equivalently from (32), 𝐹^′𝑧^(𝑘)=𝐹^′(𝑥)𝐼+4𝐶³₂𝑒^(𝑘)3+−18𝐶⁴₂+8𝐶²₂𝐶₃+6𝐶₂𝐶₃𝐶₂𝑒^(𝑘)4+𝑂𝑒^(𝑘)5. (34)
So, 𝐹^′𝑧^(𝑘)−1=𝐼+𝑊₂𝑒^(𝑘)+𝑊₃𝑒^(𝑘)2+𝑊₄𝑒^(𝑘)3+𝑊₅𝑒^(𝑘)4×𝐹^′(𝑥)⁻¹+𝑂𝑒^(𝑘)5, (35) and provided that [𝐹^′(𝑧^(𝑘))]⁻¹𝐹^′(𝑧^(𝑘))=𝐼, solving the linear system of equations involved, we have 𝑊₂=𝑊₃=0, 𝑊₄=−4𝐶³₂, and 𝑊₅=18𝐶⁴₂−8𝐶²₂𝐶₃−6𝐶₂𝐶₃𝐶₂, so 𝐹^′𝑧^(𝑘)−1=𝐼−4𝐶³₂𝑒^(𝑘)3+18𝐶⁴₂−8𝐶²₂𝐶₃−6𝐶₂𝐶₃𝐶₂𝑒^(𝑘)4×𝐹^′(𝑥)⁻¹+𝑂𝑒^(𝑘)5. (36)
Then, 𝐹^′𝑧^(𝑘)−1𝐹𝑦^(𝑘)=𝐶₂𝑒^(𝑘)2+2𝐶₃−𝐶²₂𝑒^(𝑘)3+5𝐶³₂−4𝐶₂𝐶₃−3𝐶₃𝐶₂+3𝐶₄𝑒^(𝑘)4+𝑂𝑒^(𝑘)5. (37)
Finally, by replacing (28) and (37) in the iterative expression (24), we obtain the error equation 𝑒^(𝑘+1)=𝑥^(𝑘+1)−𝑥=𝑦^(𝑘)−𝑥−𝐹^′𝑧^(𝑘)−1𝐹𝑦^(𝑘)=−𝐶³₂𝑒^(𝑘)4+𝑂𝑒^(𝑘)5, (38) and the proof of the theorem is completed.

3.2. A Fifth-Order Method

In this section, we show a new method for solving nonlinear systems that we call 𝑀5, which is obtained by combining again Newton and Traub's methods but in a different way. Its iterative expression is 𝑧^(𝑘)=𝑥^(𝑘)−𝐹^′𝑥^(𝑘)−1𝐹𝑥^(𝑘)+𝐹𝑦^(𝑘),𝑥^(𝑘+1)=𝑧^(𝑘)−𝐹^′𝑦^(𝑘)−1𝐹𝑧^(𝑘), (39) where 𝑦^(𝑘) is the 𝑘th iteration of Newton's method. We prove in the next result that the convergence order of this method is five.

Theorem 4. Let 𝐹∶𝐷⊆ℝ^𝑛→ℝ^𝑛 be sufficiently differentiable at each point of an open neighborhood 𝐷 of  𝑥∈ℝ^𝑛 that is a solution of the nonlinear system 𝐹(𝑥)=0. Let one suppose that 𝐹^′(𝑥) is continuous and nonsingular in 𝑥. Then, the sequence {𝑥^(𝑘)}_𝑘≥0 obtained by using the iterative expression (39) converges to 𝑥 with order five. The error equation is 𝑒^(𝑘+1)=4𝐶⁴₂𝑒^(𝑘)5+𝑂𝑒^(𝑘)6, (40) where 𝐶_𝑘=(1/𝑘!)[𝐹^′(𝑥)]⁻¹𝐹^(𝑘)(𝑥),𝑘=2,3,…, and 𝑒^(𝑘)=𝑥^(𝑘)−𝑥.

Proof. Following the procedure used in Theorem 3, we have that 𝐹𝑦^(𝑘)=𝐹^′(𝑥)𝐶₂𝑒^(𝑘)2+2𝐶₃−2𝐶²₂𝑒^(𝑘)3+5𝐶³₂−4𝐶₂𝐶₃−3𝐶₃𝐶₂+3𝐶₄𝑒^(𝑘)4+−6𝐶²₃−12𝐶⁴₂+10𝐶²₂𝐶₃+8𝐶₂𝐶₃𝐶₂+6𝐶₃𝐶²₂−6𝐶₂𝐶₄−4𝐶₄𝐶₂+4𝐶₅𝑒^(𝑘)5+𝑂𝑒^(𝑘)6,𝑧^(𝑘)−𝑥=2𝐶²₂𝑒^(𝑘)3+−9𝐶³₂+4𝐶₂𝐶₃+3𝐶₃𝐶₂𝑒^(𝑘)4+30𝐶⁴₂−18𝐶²₂𝐶₃+6𝐶²₃−12𝐶₃𝐶²₂−14𝐶₂𝐶₃𝐶₂+6𝐶₂𝐶₄+4𝐶₄𝐶₂𝑒^(𝑘)5+𝑂𝑒^(𝑘)6. (41)(42)
On the other hand, 𝐹𝑧^(𝑘)=𝐹^′(𝑥)𝑧^(𝑘)−𝑥+𝑂𝑒^(𝑘)4,𝐹^′𝑦^(𝑘)−1=𝐼+𝑋₂𝑦^(𝑘)−𝑥+𝑋₃𝑦^(𝑘)−𝑥²×𝐹^′(𝑥)⁻¹+𝑂𝑒^(𝑘)5. (43)(44)
Finally, by replacing (42) and (44) in the iterative expression (39), we obtain the error equation 𝑒^(𝑘+1)=𝑥^(𝑘+1)−𝑥=𝑧^(𝑘)−𝑥−𝐹^′𝑦^(𝑘)−1𝐹𝑧^(𝑘)=4𝐶⁴₂𝑒^(𝑘)5+𝑂𝑒^(𝑘)6, (45) and the proof of the theorem is completed.

3.3. Pseudocomposition

In [25] a technique called pseudocomposition that uses a known method as a predictor and the Gaussian quadrature as a corrector was introduced. The order of convergence of the resulting scheme depends, among other factors, on the order of the last two steps of the predictor. Following this idea, we generalize the procedure used to design method 𝑀5.

Then, we can establish the next result.

Theorem 5. Let 𝐹∶𝐷⊆ℝ^𝑛→ℝ^𝑛 be sufficiently differentiable at each point of an open neighborhood 𝐷 of  𝑥∈ℝ^𝑛 that is a solution of the nonlinear system 𝐹(𝑥)=0. Let one suppose that 𝐹^′(𝑥) is continuous and nonsingular in 𝑥. Let 𝑦^(𝑘) be the 𝑘th iteration of an iterative method of order 𝑞 and 𝑧^(𝑘) the 𝑘th iteration of an iterative method of order 𝑝, with 𝑝≥𝑞. The sequence {𝑥^(𝑘)}_𝑘≥0 obtained by the iterative expression 𝑥^(𝑘+1)=𝑧^(𝑘)−𝐹^′𝑦^(𝑘)−1𝐹𝑧^(𝑘), (46) converges to 𝑥 with the order of convergence 𝑝+𝑞.

Proof. Taylor’s expansions of 𝑦^(𝑘) and 𝑧^(𝑘) are 𝑦^(𝑘)=𝑥+𝑀₁𝑒^(𝑘)𝑞+𝑀₂𝑒^(𝑘)𝑞+1+⋯+𝑀_𝑝+1𝑒^{(𝑘)𝑞+𝑝}+𝑂𝑒^{(𝑘)𝑝+𝑞+1},𝑧^(𝑘)=𝑥+𝑁₁𝑒^(𝑘)𝑝+𝑁₂𝑒^(𝑘)𝑝+1+⋯+𝑁_𝑞+1𝑒^{(𝑘)𝑝+𝑞}+𝑂𝑒^{(𝑘)𝑝+𝑞+1}. (47)
Taylor’s expansion of 𝐹(𝑧^(𝑘)) around 𝑥 gives 𝐹𝑧^(𝑘)=𝐹^′(𝑥)𝑧^(𝑘)−𝑥+𝐶₂𝑧^(𝑘)−𝑥²+𝑂𝑒^{(𝑘)2(𝑝+𝑞)+1}=𝐹^′(𝑥)𝑁₁𝑒^(𝑘)𝑝+𝑁₂𝑒^(𝑘)𝑝+1+⋯+𝑁_𝑞+1𝑒^{(𝑘)𝑝+𝑞}+𝐶₂𝑁₁𝑒^(𝑘)𝑝+𝑁₂𝑒^(𝑘)𝑝+1+⋯+𝑁_𝑞+1𝑒^{(𝑘)𝑝+𝑞2}+𝑂𝑒^{(𝑘)2(𝑝+𝑞)+1}, (48) where 𝐶_𝑘=(1/𝑘!)[𝐹^′(𝑥)]⁻¹𝐹^(𝑘)(𝑥),𝑘=2,3,…, and 𝑒^(𝑘)=𝑥^(𝑘)−𝑥. On the other hand, we have that 𝐹^′𝑦^(𝑘)−1=𝐼+𝑋₂𝑦^(𝑘)−𝑥+𝑋₃𝑦^(𝑘)−𝑥²×𝐹^′(𝑥)⁻¹+𝑂𝑒^{(𝑘)2𝑞+1}=𝐼+𝑋₂𝑀₁𝑒^(𝑘)𝑞+𝑀₂𝑒^(𝑘)𝑞+1+⋯+𝑀_𝑝+1𝑒^{(𝑘)𝑞+𝑝}+𝑋₃𝑅₁𝑒^(𝑘)2𝑞×𝐹+⋯_′(𝑥)⁻¹+𝑂𝑒^(𝑘)𝑚, (49) where 𝑚=min{𝑝+𝑞+1,2𝑞+1}=𝑝+𝑞+1. Finally, by replacing (47), (48), and (49) in the iterative expression (46), we obtain the error equation 𝑒^(𝑘+1)=𝑥^(𝑘+1)−𝑥=𝑧^(𝑘)−𝑥−𝐹^′𝑦^(𝑘)−1𝐹𝑧^(𝑘)=2𝐶₂𝑀₁𝑁₁𝑒^{(𝑘)𝑝+𝑞}+𝑂𝑒^{(𝑘)𝑝+𝑞+1}. (50)
Then, the convergence order of the method that results from this combination of a method of order 𝑞 with another of order 𝑝 with 𝑝≥𝑞 is 𝑝+𝑞.

4. Numerical Results

Numerical computations have been carried out using variable precision arithmetic, with 2000 digits of mantissa, in MATLAB 7.1. The stopping criterion has been ||𝑥^(𝑘+1)−𝑥^(𝑘)||+||𝐹(𝑥^(𝑘+1))||<10⁻²⁵⁰, and therefore, we check that the iterate sequence converges to an approximation of the solution of the nonlinear system. For every method, we count the number of iterations needed to reach the wished tolerance, and we calculate the approximated computational order of convergence ACOC, the efficiency index 𝐼, the computational index 𝐶𝐼, and an error estimation made with the last values of ||𝑥^(𝑘+1)−𝑥^(𝑘)|| and ||𝐹(𝑥^(𝑘+1))||.

4.1. Numerical Results Obtained with Academic Nonlinear Systems

Now, we are going to compare 𝑀4, 𝑀5, and Newton (𝑁) and Sharma’s (𝑆) schemes with some nonlinear academic systems in order to prove the effectiveness and the computational order of convergence of the methods developed in this work. The test systems used are as follows:(a)𝐹(𝑥₁,𝑥₂)=(exp(𝑥²₁)−exp(√2𝑥₁),𝑥₁−𝑥₂),𝑥=(√2,√2)^𝑇, (b)𝐹(𝑥₁,𝑥₂)=(𝑥₁+exp(𝑥₂)−cos(𝑥₂),3𝑥₁−𝑥₂−sin(𝑥₂)),𝑥=(0,0)^𝑇, (c)𝐹(𝑥)=(𝑓₁(𝑥),𝑓₂(𝑥),…,𝑓_𝑛(𝑥)), where 𝑥=(𝑥₁,𝑥₂,…,𝑥_𝑛)^𝑇 and 𝑓_𝑖∶ℝ^𝑛→ℝ,𝑖=1,2,…,𝑛 such that 𝑓_𝑖(𝑥)=𝑥_𝑖𝑥_𝑖+1𝑓−1,𝑖=1,2,…𝑛−1,_𝑛(𝑥)=𝑥_𝑛𝑥₁−1. (51) When 𝑛 is odd, the exact zeros of 𝐹 are 𝑥₁=(1,1,…,1) and 𝑥₂=(−1,−1,…,−1).

In Table 1, we can find a comparative among the different numerical methods for the nonlinear systems (a), (b), and (c). As we can see, the approximated computational orders of convergence are the expected ones, and methods 𝑀4 and 𝑀5 are clearly very competitive in terms of error estimation.

Table 1 

Comparative of the iterative methods with nonlinear systems (a) to (c).

The efficiency index, 𝐼, and the computational index, 𝐶𝐼, of the different methods are as follows: 𝐼_𝑁=2^{1/(𝑛+𝑛2)},𝐼_𝑆=4^{1/(𝑛+2𝑛2)},𝐼_𝑀4=4^{1/(2𝑛+2𝑛2)},𝐼_𝑀5=5^{1/(3𝑛+2𝑛2)},𝐶𝐼_𝑁=2^{3/(𝑛3+3𝑛2−𝑛)},𝐶𝐼_𝑆=4^{3/(2𝑛3+15𝑛2−2𝑛)},𝐶𝐼_𝑀4=4^{3/(2𝑛3+6𝑛2−2𝑛)},𝐶𝐼_𝑀5=5^{3/(2𝑛3+6𝑛2−2𝑛)}. (52)

In Figures 2 and 3, we show these efficiency indices for 𝑛=2,3,…,10. It can be concluded that our methods improve Sharma's scheme in terms of 𝐶𝐼, although the classical efficiency of Sharma's procedure is better for 𝑛≤6. 𝑀4 and 𝑀5 are competitive, obtaining better error estimation than 𝑁 and 𝑆, with the same number of iterations.

Figure 2 

Comparative among the efficiency indices of the methods.

Figure 3 

Comparative among the computational indices of the methods.

4.2. Numerical Results for the GPS Problem

In order to test the proposed schemes on the problem of a user position of a GPS device, we have requested to the Cartographic Institute of Valencia to provide us with data of known geocentric coordinates.

Concretely, the Cartographic Institute of Valencia provided us with the following:(i)an example of a fixed-point GPS in the geocentric coordinates: 𝑥=4984687,426, 𝑦=−41199,155, and 𝑧=3966605,952. It is a point located in Alcoy (Alicante, Spain),(ii)observations from that fixed point (file *.09o) for a day,(iii)positions of the satellites for that day:*.09n and *.sp3 files,(iv)description of RINEX format (*.09o file): http://www.igs.org/components/formats.html(v)description of the ephemeris file and satellite positions sp3: http://igscb.jpl.nasa.gov/igscb/data/format/sp3c.txt(vi)link to other libraries for analysis calculations: http://www.ngs.noaa.gov/gps-toolbox/exist.htm.

With these data, we obtain the positions of the visible satellites in the instant that corresponds to the provided data. With these coordinates, we calculate the approximated pseudoranges for every satellite, and then we are able to build the nonlinear system of equations of GPS (18) using four of the satellites, with which we check the iterative methods of Newton, Sharma, 𝑀4, and 𝑀5.

In Table 2, we can find a comparative among the iterative methods 𝑁, 𝑆, 𝑀4, and 𝑀5 for the nonlinear system of the GPS. We recall that the coordinates of the center of the Earth and 𝑏_𝑢=0, that is, 𝑥⁽⁰⁾=(0,0,0,0)^𝑇, are usually used as initial estimation. Despite this, we have also tested the methods with some other initial conditions. We denote that 𝑥^∗≈(4984687.426,−41199.155,3966605.952,and0.116𝑒−8)^𝑇 as the Earth’s solution and 𝑥^𝑠≈(−39720114.893,−16748760.539,−23938190.113,and−0.159)^𝑇 as the exterior space solution.

Table 2 

Comparative of the iterative methods with the nonlinear system of equations of the GPS.

As we can see, for this particular system of equations, Newton's method does not converge to the user position for all the initial estimations, so does 𝑀5, but 𝑀4 is a good method in all senses, very competitive in respect of known methods.

5. Conclusions

In this paper, we have gone in depth on an emerging line of investigation, the GPS receivers software improvement. Concretely, GPS receivers currently use Newton's method to solve the nonlinear system (18) and to calculate their exact position with the information obtained from signals received from the GPS constellation of satellites. We propose two different combinations of Newton method and Traub's methods, obtaining two methods of fourth- (𝑀4) and fifth-order (𝑀5). Using the idea presented in [25], called pseudocomposition, it is proved that combining in a particular way two methods of order 𝑝 and 𝑞, respectively, with 𝑝≥𝑞, the order of convergence of the resulting scheme is 𝑝+𝑞. We have numerically compared the different methods, and we have concluded that 𝑀4 and 𝑀5 are very competitive in terms of the error estimation.

Acknowledgments

This research was supported by Ministerio de Ciencia y Tecnología MTM2011-28636-C02-02 and FONDOCYT 2011-1-B1-33 República Dominicana. The authors would also like to thank the work of the anonymous referee.