Constructing Frozen Jacobian Iterative Methods for Solving Systems of Nonlinear Equations, Associated with ODEs and PDEs Using the Homotopy Method

Abstract

A homotopy method is presented for the construction of frozen Jacobian iterative methods. The frozen Jacobian iterative methods are attractive because the inversion of the Jacobian is performed in terms of LUfactorization only once, for a single instance of the iterative method. We embedded parameters in the iterative methods with the help of the homotopy method: the values of the parameters are determined in such a way that a better convergence rate is achieved. The proposed homotopy technique is general and has the ability to construct different families of iterative methods, for solving weakly nonlinear systems of equations. Further iterative methods are also proposed for solving general systems of nonlinear equations.

1. Introduction

Iterative methods for solving nonlinear equations and systems of nonlinear equations always represent an attractive research topic. Newton’s method [1,2,3] is a classical iterative method for solving systems of nonlinear equations and offers quadratic convergence. Many authors have developed further iterative methods for solving either nonlinear equations or systems of nonlinear equations [4,5,6,7]. Here, we focus on the idea of the frozen Jacobian: we observe that this idea is computationally attractive for designing iterative methods because the inversion of the Jacobian (regarding LUfactorization) is performed a single time, for a single instance of the iterative method. Many researchers have proposed frozen Jacobian iterations (see [8,9,10,11,12] and references therein). In our view, the use of perturbation techniques in connection with the homotopy method can be very effective for developing new frozen Jacobian iterative methods, when solving systems of nonlinear equations. The homotopy method can be found in the work of Poincaré [13]. After a careful review of the relevant literature, we come to know that much work has been done in the direction of designing iterative methods for solving nonlinear scalar equations, using the homotopy perturbation method, but not for systems of nonlinear equations. The reason behind this fact is that usually, the constructed iterative methods, for solving nonlinear equations using different homotopy methods, contain higher order derivatives. From a computational point of view, the evaluation of higher order Fréchet derivatives is expensive for a general system of nonlinear equations. However, when considering systems of weakly nonlinear equations, associated with wide classes of ordinary and partial differential equations, we find that all higher order Fréchet derivatives possess a diagonal structure: as a consequence, the resulting computational cost is mild.

The core idea of homotopy (or homotopy continuation) was developed for solving hard problems numerically. The nonlinear problem is approximated around some initial guess to construct a linear model or a model that can be solved easily. Once we have a linear model, we construct a homotopy function that maps linear model to a nonlinear model with the help of auxiliary parameters. The embedding of auxiliary parameters in the homotopy function helps us to solve stiff problems iteratively in a way that the solution of the less hard problem is used as an initial guess for the next harder problem. The gradual change in the parameters assists us in solving a hard problem softly in the mathematical sense. The above-described scenario is one of the pictures of the homotopy methods or, precisely speaking, homotopy continuation methods. The other picture is to use the homotopy method to approximate solution of nonlinear problems analytically [14,15,16,17,18,19,20]. Many authors have used [7,21] homotopy methods for constructing iterative methods for solving nonlinear equations. For the construction of iterative methods, they build a homotopy between a linear model of the nonlinear equation and the nonlinear equation, by expanding it around some initial guess. The constructed homotopy contains auxiliary parameters. The auxiliary parameters can be used for different purposes. Some authors use the homotopy parameters to enhance the convergence by changing the dynamics of the iterative method and give estimates for the optimal parameters [22].

Now, we present in detail our proposals. Noor [21] proposed the following homotopy iterative method for solving nonlinear equations. Let $\mathrm{f}\left(\mathrm{x}\right)=0$ be a scalar nonlinear equation, and let ${\mathrm{x}}_0$ be an initial guess. We perform a linearization of $\mathrm{f}\left(\mathrm{x}\right)$ around a point γ:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \mathrm{f}\left(\mathrm{x}\right)& =\mathrm{f}\left(\mathsf{\gamma }\right)+{\mathrm{f}}^{\prime }\left(\mathsf{\gamma }\right)(\mathrm{x}-\mathsf{\gamma })+\mathrm{g}\left(\mathrm{x}\right)=0\hfill \\ \hfill \mathrm{g}\left(\mathrm{x}\right)& =\mathrm{f}\left(\mathrm{x}\right)-\mathrm{f}\left(\mathsf{\gamma }\right)-{\mathrm{f}}^{\prime }\left(\mathsf{\gamma }\right)(\mathrm{x}-\mathsf{\gamma })\hfill \end{array}\end{array}$$

(1)

where $\mathrm{f}\left(\mathsf{\gamma }\right)+{\mathrm{f}}^{\prime }\left(\mathsf{\gamma }\right)(\mathrm{x}-\mathsf{\gamma })$ is the linear approximation of $\mathrm{f}\left(\mathrm{x}\right)$ around γ and $\mathrm{g}\left(\mathrm{x}\right)$ is the nonlinear part of $\mathrm{f}\left(\mathrm{x}\right)$. Equation (1) can be written as:

$$\begin{array}{c}\hfill \mathrm{x}=\mathrm{c}+\mathrm{N}\left(\mathrm{x}\right)\end{array}$$

(2)

where $\mathrm{c}=\mathsf{\gamma }-\frac{\mathrm{f}\left(\mathsf{\gamma }\right)}{{\mathrm{f}}^{\prime }\left(\mathsf{\gamma }\right)}$ and $\mathrm{N}\left(\mathrm{x}\right)=-\frac{\mathrm{g}\left(\mathrm{x}\right)}{{\mathrm{f}}^{\prime }\left(\mathsf{\gamma }\right)}$. When assuming that L is a linear function, according to [21], a homotopy can be defined between $\mathrm{L}\left(\mathrm{v}\right)-\mathrm{L}\left({\mathrm{x}}_0\right)$ and $\mathrm{L}\left(\mathrm{v}\right)+\mathrm{N}\left(\mathrm{v}\right)-\mathrm{c}$:

$$\begin{array}{c}\hfill \mathrm{H}(\mathrm{v},\mathrm{q},\mathrm{T})=(1-\mathsf{\alpha }\mathrm{q})\left(\mathrm{L}\left(\mathrm{v}\right)-\mathrm{L}\left({\mathrm{x}}_0\right)\right)+\mathrm{q}\mathsf{\alpha }\left(\mathrm{L}\left(\mathrm{v}\right)+\mathrm{N}\left(\mathrm{v}\right)-\mathrm{c}\right)-\mathsf{\alpha }{\mathrm{q}}^2(1-\mathrm{q})\mathrm{T}=0\end{array}$$

(3)

where $\mathrm{H}(\mathrm{v},\mathrm{q},\mathrm{T}):\mathbb{R}\times [0,1]\times \mathbb{R}⟶\mathbb{R}$, $\mathrm{q}\in [0,1]$ is an embedded parameter, α is a real free parameter and T is a suitable operator. In Equation (3), $\mathrm{L}\left(\mathrm{v}\right)-\mathrm{L}\left({\mathrm{x}}_0\right)$ and $\mathrm{L}\left(\mathrm{v}\right)+\mathrm{N}\left(\mathrm{v}\right)-\mathrm{c}$ are called homotopic, i.e., we established homotopy between the linear model and the nonlinear model. When $\mathrm{q}=0$, we have to solve the linear model $\mathrm{H}(\mathrm{v},0,\mathrm{T})=\mathrm{L}\left(\mathrm{v}\right)-\mathrm{L}\left({\mathrm{x}}_0\right)$, and for $\mathrm{q}=1$, we have as the target the nonlinear system in Equation (2), that is $\mathrm{H}(\mathrm{v},1,\mathrm{T})=\mathsf{\alpha }\left(\mathrm{N}\left(\mathrm{v}\right)-\mathrm{c}\right)$. When we change q from zero to one, the linear problem continuously deforms into the original nonlinear problem. Further, we assume that the solution of Equation (3) can be expressed as power series in q:

$$\begin{array}{c}\hfill \mathrm{v}={\mathrm{x}}_0+\mathrm{q}{\mathrm{x}}_1+{\mathrm{q}}^2{\mathrm{x}}_2+\cdots \end{array}$$

(4)

Notice that the solution of target nonlinear problem can be approximated as:

$$\begin{array}{c}\hfill \mathrm{x}=\underset{\mathrm{q}\to 1}{lim}\mathrm{v}={\mathrm{x}}_0+{\mathrm{x}}_1+{\mathrm{x}}_2+\cdots \end{array}$$

(5)

In 2010, Yongyan et al. [22] proposed a two-parameter homotopy method for the construction of iterative methods, for solving scalar nonlinear equations. The authors in [22] use the homotopy auxiliary parameters to change the dynamics of the proposed iterative methods for solving nonlinear equations and provide the optimal values of auxiliary parameters, for enhancing the convergence speed without altering the convergence order of the iterative method. It is worth noting that they used auxiliary parameters to change the path of converging sequences of successive approximations to obtain fast convergence. Inspired by their mathematical model for homotopy, we introduce a generalization of the two-parameter homotopy for solving a system of weakly nonlinear equations. However, our primary motivation is to use the auxiliary parameters to increase the convergence order of the iterative methods. Secondly, our proposal is for the construction of iterative methods for solving a system of nonlinear equations instead of a single nonlinear equation. In summary, we use the homotopy techniques to develop mathematical models that we use for designing higher order iterative methods with the help of new parameters.

2. Parametric Homotopy

Denote a system of nonlinear equations in a general from as:

$$\begin{array}{c}\hfill \mathbf{F}\left(\mathbf{y}\right)=\mathbf{0}\end{array}$$

(6)

where $\mathbf{y}={[{\mathrm{y}}_1,{\mathrm{y}}_2,\cdots ,{\mathrm{y}}_{\mathrm{n}}]}^{\mathrm{T}}$ and $\mathbf{F}\left(\mathbf{y}\right)={[{\mathrm{F}}_1\left(\mathbf{y}\right),{\mathrm{F}}_2\left(\mathbf{y}\right),\cdots ,{\mathrm{F}}_{\mathrm{n}}\left(\mathbf{y}\right)]}^{\mathrm{T}}$. We are interested in finding a simple root ${\mathbf{y}}^{*}$ of Equation (6), i.e., $\mathbf{F}\left({\mathbf{y}}^{*}\right)=\mathbf{0}$ and det$\left({\mathbf{F}}^{\prime }\left({\mathbf{y}}^{*}\right)\right)\ne 0$. Let ${\mathbf{y}}_{\mathsf{\vartheta }}$ be an initial guess. Hence, we construct the linear approximation of Equation (6) as:

$$\begin{array}{c}\hfill \mathbf{F}\left(\mathbf{y}\right)\approx \mathbf{L}\left(\mathbf{y}\right)=\mathbf{F}\left({\mathbf{y}}_{\mathsf{\vartheta }}\right)+{\mathbf{F}}^{\prime }\left({\mathbf{y}}_{\mathsf{\vartheta }}\right)\left(\mathbf{y}-{\mathbf{y}}_{\mathsf{\vartheta }}\right)\end{array}$$

(7)

Now, we define a parametric homotopy between $\mathbf{L}\left(\mathbf{y}\right)$ and $\mathbf{F}\left(\mathbf{y}\right)$, taking inspiration from the work in [22]. Therefore:

$$\begin{array}{c}\hfill \left(1+{\displaystyle \sum _{\mathrm{i}=1}^{\mathrm{m}}}{\mathsf{\alpha }}_{\mathrm{i}}{\mathrm{q}}^{\mathrm{i}}\right)\mathbf{L}\left(\mathbf{y}\right)-{\mathsf{\alpha }}_0\mathrm{q}\mathbf{F}\left(\mathbf{y}\left(\mathrm{q}\right)\right)=\mathbf{0}\end{array}$$

(8)

where the quantities ${\mathsf{\alpha }}_{\mathrm{i}}$, $\mathrm{i}=0,\dots ,\mathrm{m}$, are parameters, ${\mathsf{\alpha }}_0\ne 0$, $1+{\displaystyle \sum _{\mathrm{i}=1}^{\mathrm{m}}}{\mathsf{\alpha }}_{\mathrm{i}}=0$, $\mathbf{y}\left(\mathrm{q}\right)={\mathbf{x}}_0+\mathrm{q}{\mathbf{x}}_1+{\mathrm{q}}^2{\mathbf{x}}_2+\cdots$ and $\mathbf{y}\left(1\right)={\mathbf{y}}^{*}$.

3. Construction of Iterative Methods

We start from the case where $\mathrm{m}=3$ and when Equation (8) can be written as:

$$\begin{array}{c}\hfill \left(1+{\mathsf{\alpha }}_1\mathrm{q}+{\mathsf{\alpha }}_2{\mathrm{q}}^2-(1+{\mathsf{\alpha }}_1+{\mathsf{\alpha }}_2){\mathrm{q}}^3\right)\mathbf{L}\left(\mathbf{y}\right)-{\mathsf{\alpha }}_0\mathrm{q}\mathbf{F}\left(\mathbf{y}\left(\mathrm{q}\right)\right)=\mathbf{0}\end{array}$$

(9)

The Taylor series expansion of $\mathbf{F}(·)$ around ${\mathbf{x}}_0$ gives:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \mathbf{F}\left(\mathbf{y}\right)& =\mathbf{F}\left({\mathbf{x}}_0\right)+{\mathbf{F}}^{\prime }\left({\mathbf{x}}_0\right)\left(\mathbf{y}-{\mathbf{x}}_0\right)+1/2{\mathbf{F}}^{″}\left({\mathbf{x}}_0\right){\left(\mathbf{y}-{\mathbf{x}}_0\right)}^2+1/6{\mathbf{F}}^{‴}\left({\mathbf{x}}_0\right){\left(\mathbf{y}-{\mathbf{x}}_0\right)}^3+\mathrm{O}{\left(\mathbf{y}-{\mathbf{x}}_0\right)}^4\hfill \\ \hfill \mathbf{F}\left(\mathbf{y}\right)& ={\mathbf{C}}_0+{\mathbf{C}}_1\mathbf{z}+1/2{\mathbf{C}}_2{\mathbf{z}}^2+1/6{\mathbf{C}}_3{\mathbf{z}}^3+\mathrm{O}{\left(\mathbf{z}\right)}^4\hfill \end{array}\end{array}$$

(10)

where $\mathbf{z}=\mathbf{y}-{\mathbf{x}}_0$ and ${\mathbf{C}}_{\mathrm{j}}={\mathbf{F}}^{\left(\mathrm{j}\right)}\left({\mathbf{x}}_0\right)$, i.e., ${\mathbf{C}}_{\mathrm{j}}$ is a j-th order Fréchet derivative. Using $\mathbf{y}\left(\mathrm{q}\right)={\mathbf{x}}_0+{\displaystyle \sum _{\mathrm{i}=1}^{\infty }}{\mathrm{q}}^{\mathrm{i}}{\mathbf{x}}_{\mathrm{i}}$ in Equation (10), we get:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \mathbf{F}\left(\mathbf{y}\right)& ={\mathbf{C}}_0+{\mathbf{C}}_1{\mathbf{x}}_1\mathrm{q}+\left({\mathbf{C}}_1{\mathbf{x}}_2+1/2{\mathbf{C}}_2{{\mathbf{x}}_1}^2\right){\mathrm{q}}^2+\left({\mathbf{C}}_2{\mathbf{x}}_2{\mathbf{x}}_1+{\mathbf{C}}_1{\mathbf{x}}_3+1/6{\mathbf{C}}_3{{\mathbf{x}}_1}^3\right){\mathrm{q}}^3+({\mathbf{C}}_1{\mathbf{x}}_4\hfill \\ & +{\mathbf{C}}_2{\mathbf{x}}_3{\mathbf{x}}_1+1/2{\mathbf{C}}_2{{\mathbf{x}}_2}^2+1/2{\mathbf{C}}_3{{\mathbf{x}}_1}^2{\mathbf{x}}_2){\mathrm{q}}^4+({\mathbf{C}}_2{\mathbf{x}}_4{\mathbf{x}}_1+{\mathbf{C}}_2{\mathbf{x}}_3{\mathbf{x}}_2+1/2{\mathbf{C}}_3{{\mathbf{x}}_1}^2{\mathbf{x}}_3\hfill \\ & +1/2{\mathbf{C}}_3{\mathbf{x}}_1{{\mathbf{x}}_2}^2){\mathrm{q}}^5+\mathrm{O}\left({\mathrm{q}}^6\right)\hfill \end{array}\end{array}$$

(11)

The linear approximation $\mathbf{L}\left(\mathbf{y}\right)$ can be written as:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \mathbf{L}\left(\mathbf{y}\right)& ={\mathbf{F}}^{\prime }\left({\mathbf{y}}_{\mathsf{\vartheta }}\right){\mathbf{x}}_0-{\mathbf{F}}^{\prime }\left({\mathbf{y}}_{\mathsf{\vartheta }}\right){\mathbf{y}}_{\mathsf{\vartheta }}+\mathbf{F}\left({\mathbf{y}}_{\mathsf{\vartheta }}\right)+{\mathbf{F}}^{\prime }\left({\mathbf{y}}_{\mathsf{\vartheta }}\right){\mathbf{x}}_1\mathrm{q}+{\mathbf{F}}^{\prime }\left({\mathbf{y}}_{\mathsf{\vartheta }}\right){\mathbf{x}}_2{\mathrm{q}}^2+{\mathbf{F}}^{\prime }\left({\mathbf{y}}_{\mathsf{\vartheta }}\right){\mathbf{x}}_3{\mathrm{q}}^3+{\mathbf{F}}^{\prime }\left({\mathbf{y}}_{\mathsf{\vartheta }}\right){\mathbf{x}}_4{\mathrm{q}}^4\hfill \\ & +{\mathbf{F}}^{\prime }\left({\mathbf{y}}_{\mathsf{\vartheta }}\right){\mathbf{x}}_5{\mathrm{q}}^5+\mathrm{O}\left({\mathrm{q}}^6\right)\hfill \end{array}\end{array}$$

(12)

Simplification of Equation (9) gives:

$$\begin{array}{c}\hfill \begin{array}{cc}& \mathbf{F}\left({\mathbf{y}}_{\mathsf{\vartheta }}\right)+{\mathbf{F}}^{\prime }\left({\mathbf{y}}_{\mathsf{\vartheta }}\right){\mathbf{x}}_0-{\mathbf{F}}^{\prime }\left({\mathbf{y}}_{\mathsf{\vartheta }}\right){\mathbf{y}}_{\mathsf{\vartheta }}+\left({\mathsf{\alpha }}_1{\mathbf{F}}^{\prime }\left({\mathbf{y}}_{\mathsf{\vartheta }}\right){\mathbf{x}}_0-{\mathsf{\alpha }}_1{\mathbf{F}}^{\prime }\left({\mathbf{y}}_{\mathsf{\vartheta }}\right){\mathbf{y}}_{\mathsf{\vartheta }}-{\mathsf{\alpha }}_0{\mathbf{C}}_0+{\mathsf{\alpha }}_1\mathbf{F}\left({\mathbf{y}}_{\mathsf{\vartheta }}\right)+{\mathbf{F}}^{\prime }\left({\mathbf{y}}_{\mathsf{\vartheta }}\right){\mathbf{x}}_1\right)\mathrm{q}\hfill \\ & +\left(-{\mathsf{\alpha }}_0{\mathbf{C}}_1{\mathbf{x}}_1+{\mathsf{\alpha }}_1{\mathbf{F}}^{\prime }\left({\mathbf{y}}_{\mathsf{\vartheta }}\right){\mathbf{x}}_1+{\mathsf{\alpha }}_2{\mathbf{F}}^{\prime }\left({\mathbf{y}}_{\mathsf{\vartheta }}\right){\mathbf{x}}_0-{\mathsf{\alpha }}_2{\mathbf{F}}^{\prime }\left({\mathbf{y}}_{\mathsf{\vartheta }}\right){\mathbf{y}}_{\mathsf{\vartheta }}+{\mathsf{\alpha }}_2\mathbf{F}\left({\mathbf{y}}_{\mathsf{\vartheta }}\right)+{\mathbf{F}}^{\prime }\left({\mathbf{y}}_{\mathsf{\vartheta }}\right){\mathbf{x}}_2\right){\mathrm{q}}^2\hfill \\ & +(-{\mathsf{\alpha }}_0{\mathbf{C}}_1{\mathbf{x}}_2-{\mathsf{\alpha }}_1{\mathbf{F}}^{\prime }\left({\mathbf{y}}_{\mathsf{\vartheta }}\right){\mathbf{x}}_0+{\mathsf{\alpha }}_1{\mathbf{F}}^{\prime }\left({\mathbf{y}}_{\mathsf{\vartheta }}\right){\mathbf{x}}_2+{\mathsf{\alpha }}_1{\mathbf{F}}^{\prime }\left({\mathbf{y}}_{\mathsf{\vartheta }}\right){\mathbf{y}}_0-{\mathsf{\alpha }}_2{\mathbf{F}}^{\prime }\left({\mathbf{y}}_{\mathsf{\vartheta }}\right){\mathbf{x}}_0+{\mathsf{\alpha }}_2{\mathbf{F}}^{\prime }\left({\mathbf{y}}_{\mathsf{\vartheta }}\right){\mathbf{x}}_1+{\mathsf{\alpha }}_2{\mathbf{F}}^{\prime }\left({\mathbf{y}}_{\mathsf{\vartheta }}\right){\mathbf{y}}_{\mathsf{\vartheta }}\hfill \\ & -{\mathsf{\alpha }}_1\mathbf{F}\left({\mathbf{y}}_{\mathsf{\vartheta }}\right)-{\mathsf{\alpha }}_2\mathbf{F}\left({\mathbf{y}}_{\mathsf{\vartheta }}\right)-{\mathbf{F}}^{\prime }\left({\mathbf{y}}_{\mathsf{\vartheta }}\right){\mathbf{x}}_0+{\mathbf{F}}^{\prime }\left({\mathbf{y}}_{\mathsf{\vartheta }}\right){\mathbf{x}}_3+{\mathbf{F}}^{\prime }\left({\mathbf{y}}_{\mathsf{\vartheta }}\right){\mathbf{y}}_{\mathsf{\vartheta }}-\mathbf{F}\left({\mathbf{y}}_{\mathsf{\vartheta }}\right)-1/2{\mathsf{\alpha }}_0{\mathbf{C}}_2{{\mathbf{x}}_1}^2){\mathrm{q}}^3+(-{\mathsf{\alpha }}_0{\mathbf{C}}_2{\mathbf{x}}_2{\mathbf{x}}_1\hfill \\ & -{\mathsf{\alpha }}_0{\mathbf{C}}_1{\mathbf{x}}_3-{\mathsf{\alpha }}_1{\mathbf{F}}^{\prime }\left({\mathbf{y}}_{\mathsf{\vartheta }}\right){\mathbf{x}}_1+{\mathsf{\alpha }}_1{\mathbf{F}}^{\prime }\left({\mathbf{y}}_{\mathsf{\vartheta }}\right){\mathbf{x}}_3-{\mathsf{\alpha }}_2{\mathbf{F}}^{\prime }\left({\mathbf{y}}_{\mathsf{\vartheta }}\right){\mathbf{x}}_1+{\mathsf{\alpha }}_2{\mathbf{F}}^{\prime }\left({\mathbf{y}}_{\mathsf{\vartheta }}\right){\mathbf{x}}_2-{\mathbf{F}}^{\prime }\left({\mathbf{y}}_{\mathsf{\vartheta }}\right){\mathbf{x}}_1+{\mathbf{F}}^{\prime }\left({\mathbf{y}}_{\mathsf{\vartheta }}\right){\mathbf{x}}_4\hfill \\ & -1/6{\mathsf{\alpha }}_0{\mathbf{C}}_3{{\mathbf{x}}_1}^3){\mathrm{q}}^4+\mathrm{O}\left({\mathrm{q}}^5\right)=\mathbf{0}\hfill \end{array}\end{array}$$

(13)

By equating to zero the coefficients of powers of q, we find five equations:

$$\begin{array}{cc}& {\mathbf{F}}^{\prime }\left({\mathbf{y}}_{\mathsf{\vartheta }}\right){\mathbf{x}}_0-{\mathbf{F}}^{\prime }\left({\mathbf{y}}_{\mathsf{\vartheta }}\right){\mathbf{y}}_{\mathsf{\vartheta }}+\mathbf{F}\left({\mathbf{y}}_{\mathsf{\vartheta }}\right)=\mathbf{0}\hfill \end{array}$$

(14)

$$\begin{array}{cc}& {\mathsf{\alpha }}_1{\mathbf{F}}^{\prime }\left({\mathbf{y}}_{\mathsf{\vartheta }}\right){\mathbf{x}}_0-{\mathsf{\alpha }}_1{\mathbf{F}}^{\prime }\left({\mathbf{y}}_{\mathsf{\vartheta }}\right){\mathbf{y}}_{\mathsf{\vartheta }}-{\mathsf{\alpha }}_0{\mathbf{C}}_0+{\mathsf{\alpha }}_1\mathbf{F}\left({\mathbf{y}}_{\mathsf{\vartheta }}\right)+{\mathbf{F}}^{\prime }\left({\mathbf{y}}_{\mathsf{\vartheta }}\right){\mathbf{x}}_1=\mathbf{0}\hfill \end{array}$$

(15)

$$\begin{array}{cc}& -{\mathsf{\alpha }}_0{\mathbf{C}}_1{\mathbf{x}}_1+{\mathsf{\alpha }}_1{\mathbf{F}}^{\prime }\left({\mathbf{y}}_{\mathsf{\vartheta }}\right){\mathbf{x}}_1+{\mathsf{\alpha }}_2{\mathbf{F}}^{\prime }\left({\mathbf{y}}_{\mathsf{\vartheta }}\right){\mathbf{x}}_0-{\mathsf{\alpha }}_2{\mathbf{F}}^{\prime }\left({\mathbf{y}}_{\mathsf{\vartheta }}\right){\mathbf{y}}_{\mathsf{\vartheta }}+{\mathsf{\alpha }}_2\mathbf{F}\left({\mathbf{y}}_{\mathsf{\vartheta }}\right)+{\mathbf{F}}^{\prime }\left({\mathbf{y}}_{\mathsf{\vartheta }}\right){\mathbf{x}}_2=\mathbf{0}\hfill \\ & {\mathbf{F}}^{\prime }\left({\mathbf{y}}_{\mathsf{\vartheta }}\right){\mathbf{x}}_3+{\mathsf{\alpha }}_1{\mathbf{F}}^{\prime }\left({\mathbf{y}}_{\mathsf{\vartheta }}\right){\mathbf{x}}_2+{\mathsf{\alpha }}_2{\mathbf{F}}^{\prime }\left({\mathbf{y}}_{\mathsf{\vartheta }}\right){\mathbf{x}}_1-\mathbf{F}\left({\mathbf{y}}_{\mathsf{\vartheta }}\right)-{\mathbf{F}}^{\prime }\left({\mathbf{y}}_{\mathsf{\vartheta }}\right){\mathbf{x}}_0+{\mathbf{F}}^{\prime }\left({\mathbf{y}}_{\mathsf{\vartheta }}\right){\mathbf{y}}_{\mathsf{\vartheta }}-{\mathsf{\alpha }}_1\mathbf{F}\left({\mathbf{y}}_{\mathsf{\vartheta }}\right)-{\mathsf{\alpha }}_1{\mathbf{F}}^{\prime }\left({\mathbf{y}}_{\mathsf{\vartheta }}\right){\mathbf{x}}_0\hfill \end{array}$$

(16)

$$\begin{array}{cc}& +{\mathsf{\alpha }}_1{\mathbf{F}}^{\prime }\left({\mathbf{y}}_{\mathsf{\vartheta }}\right){\mathbf{y}}_{\mathsf{\vartheta }}-{\mathsf{\alpha }}_2\mathbf{F}\left({\mathbf{y}}_{\mathsf{\vartheta }}\right)-{\mathsf{\alpha }}_2{\mathbf{F}}^{\prime }\left({\mathbf{y}}_{\mathsf{\vartheta }}\right){\mathbf{x}}_0+{\mathsf{\alpha }}_2{\mathbf{F}}^{\prime }\left({\mathbf{y}}_{\mathsf{\vartheta }}\right){\mathbf{y}}_{\mathsf{\vartheta }}-{\mathsf{\alpha }}_0{\mathbf{C}}_1{\mathbf{x}}_2-1/2{\mathsf{\alpha }}_0{\mathbf{C}}_2{{\mathbf{x}}_1}^2=\mathbf{0}\hfill \\ & {\mathbf{F}}^{\prime }\left({\mathbf{y}}_{\mathsf{\vartheta }}\right){\mathbf{x}}_4+{\mathsf{\alpha }}_1{\mathbf{F}}^{\prime }\left({\mathbf{y}}_{\mathsf{\vartheta }}\right){\mathbf{x}}_3+{\mathsf{\alpha }}_2{\mathbf{F}}^{\prime }\left({\mathbf{y}}_{\mathsf{\vartheta }}\right){\mathbf{x}}_2-{\mathbf{F}}^{\prime }\left({\mathbf{y}}_{\mathsf{\vartheta }}\right){\mathbf{x}}_1-{\mathsf{\alpha }}_1{\mathbf{F}}^{\prime }\left({\mathbf{y}}_{\mathsf{\vartheta }}\right){\mathbf{x}}_1-{\mathsf{\alpha }}_2{\mathbf{F}}^{\prime }\left({\mathbf{y}}_{\mathsf{\vartheta }}\right){\mathbf{x}}_1-{\mathsf{\alpha }}_0{\mathbf{C}}_1{\mathbf{x}}_3-{\mathsf{\alpha }}_0{\mathbf{C}}_2{\mathbf{x}}_2{\mathbf{x}}_1\hfill \end{array}$$

(17)

$$\begin{array}{cc}& -1/6{\mathsf{\alpha }}_0{\mathbf{C}}_3{{\mathbf{x}}_1}^3=\mathbf{0}\hfill \end{array}$$

(18)

By solving Equations (14)–(18), we obtain:

$$\begin{array}{cc}& {\mathbf{x}}_0={\mathbf{y}}_{\mathsf{\vartheta }}-{\mathbf{F}}^{\prime }{\left({\mathbf{y}}_{\mathsf{\vartheta }}\right)}^{-1}\mathbf{F}\left({\mathbf{y}}_{\mathsf{\vartheta }}\right)\hfill \\ & {\mathbf{x}}_1=-{\mathsf{\alpha }}_0{\mathbf{F}}^{\prime }{\left({\mathbf{y}}_{\mathsf{\vartheta }}\right)}^{-1}\mathbf{F}\left({\mathbf{x}}_0\right)\hfill \\ & {\mathbf{x}}_2={\mathsf{\alpha }}_0{\mathsf{\alpha }}_1{\mathbf{F}}^{\prime }{\left({\mathbf{y}}_{\mathsf{\vartheta }}\right)}^{-1}\mathbf{F}\left({\mathbf{x}}_0\right)+{\mathsf{\alpha }}_0^2{\mathbf{F}}^{\prime }{\left({\mathbf{y}}_{\mathsf{\vartheta }}\right)}^{-1}{\mathbf{F}}^{\prime }\left({\mathbf{x}}_0\right){\mathbf{F}}^{\prime }{\left({\mathbf{y}}_{\mathsf{\vartheta }}\right)}^{-1}\mathbf{F}\left({\mathbf{x}}_0\right)\hfill \\ & {\mathbf{x}}_3=-{\mathsf{\alpha }}_1^2{\mathsf{\alpha }}_0{\mathbf{F}}^{\prime }{\left({\mathbf{y}}_{\mathsf{\vartheta }}\right)}^{-1}\mathbf{F}\left({\mathbf{x}}_0\right)+{\mathsf{\alpha }}_2{\mathsf{\alpha }}_0{\mathbf{F}}^{\prime }{\left({\mathbf{y}}_{\mathsf{\vartheta }}\right)}^{-1}\mathbf{F}\left({\mathbf{x}}_0\right)+{\mathsf{\alpha }}_0^3{\mathbf{F}}^{\prime }{\left({\mathbf{y}}_{\mathsf{\vartheta }}\right)}^{-1}{\mathbf{F}}^{\prime }\left({\mathbf{x}}_0\right){\mathbf{F}}^{\prime }{\left({\mathbf{y}}_{\mathsf{\vartheta }}\right)}^{-1}{\mathbf{F}}^{\prime }\left({\mathbf{x}}_0\right){\mathbf{F}}^{\prime }{\left({\mathbf{y}}_{\mathsf{\vartheta }}\right)}^{-1}\mathbf{F}\left({\mathbf{x}}_0\right)\hfill \\ & +1/2{\mathsf{\alpha }}_0^3{\mathbf{F}}^{\prime }{\left({\mathbf{y}}_{\mathsf{\vartheta }}\right)}^{-1}{\mathbf{F}}^{″}\left({\mathbf{x}}_0\right){\mathbf{F}}^{\prime }{\left({\mathbf{y}}_{\mathsf{\vartheta }}\right)}^{-1}\mathbf{F}\left({\mathbf{x}}_0\right){\mathbf{F}}^{\prime }{\left({\mathbf{y}}_{\mathsf{\vartheta }}\right)}^{-1}\mathbf{F}\left({\mathbf{x}}_0\right)\hfill \end{array}$$

$$\begin{array}{cc}& {\mathbf{x}}_4={\mathsf{\alpha }}_1^3{\mathsf{\alpha }}_0{\mathbf{F}}^{\prime }{\left({\mathbf{y}}_{\mathsf{\vartheta }}\right)}^{-1}\mathbf{F}\left({\mathbf{x}}_0\right)-2{\mathsf{\alpha }}_1{\mathsf{\alpha }}_2{\mathsf{\alpha }}_0{\mathbf{F}}^{\prime }{\left({\mathbf{y}}_{\mathsf{\vartheta }}\right)}^{-1}\mathbf{F}\left({\mathbf{x}}_0\right)-{\mathsf{\alpha }}_1\left({\mathsf{\alpha }}_0^3\right){\mathbf{F}}^{\prime }{\left({\mathbf{y}}_{\mathsf{\vartheta }}\right)}^{-1}{\mathbf{F}}^{\prime }\left({\mathbf{x}}_0\right){\mathbf{F}}^{\prime }{\left({\mathbf{y}}_{\mathsf{\vartheta }}\right)}^{-1}{\mathbf{F}}^{\prime }\left({\mathbf{x}}_0\right)\hfill \\ & {\mathbf{F}}^{\prime }{\left({\mathbf{y}}_{\mathsf{\vartheta }}\right)}^{-1}\mathbf{F}\left({\mathbf{x}}_0\right)-3/2{\mathsf{\alpha }}_1{\mathsf{\alpha }}_0^3{\mathbf{F}}^{\prime }{\left({\mathbf{y}}_{\mathsf{\vartheta }}\right)}^{-1}{\mathbf{F}}^{″}\left({\mathbf{x}}_0\right){\mathbf{F}}^{\prime }{\left({\mathbf{y}}_{\mathsf{\vartheta }}\right)}^{-1}\mathbf{F}\left({\mathbf{x}}_0\right){\mathbf{F}}^{\prime }{\left({\mathbf{y}}_{\mathsf{\vartheta }}\right)}^{-1}\mathbf{F}\left({\mathbf{x}}_0\right)-{\mathsf{\alpha }}_0{\mathbf{F}}^{\prime }{\left({\mathbf{y}}_{\mathsf{\vartheta }}\right)}^{-1}\mathbf{F}\left({\mathbf{x}}_0\right)\hfill \\ & -{\mathsf{\alpha }}_0{\mathsf{\alpha }}_1{\mathbf{F}}^{\prime }{\left({\mathbf{y}}_{\mathsf{\vartheta }}\right)}^{-1}\mathbf{F}\left({\mathbf{x}}_0\right)-{\mathsf{\alpha }}_2{\mathsf{\alpha }}_0{\mathbf{F}}^{\prime }{\left({\mathbf{y}}_{\mathsf{\vartheta }}\right)}^{-1}\mathbf{F}\left({\mathbf{x}}_0\right)-{\mathsf{\alpha }}_1^2\left({\mathsf{\alpha }}_0^2\right){\mathbf{F}}^{\prime }{\left({\mathbf{y}}_{\mathsf{\vartheta }}\right)}^{-1}{\mathbf{F}}^{\prime }\left({\mathbf{x}}_0\right){\mathbf{F}}^{\prime }{\left({\mathbf{y}}_{\mathsf{\vartheta }}\right)}^{-1}\mathbf{F}\left({\mathbf{x}}_0\right)\hfill \\ & +{\mathsf{\alpha }}_0^4{\mathbf{F}}^{\prime }{\left({\mathbf{y}}_{\mathsf{\vartheta }}\right)}^{-1}{\mathbf{F}}^{\prime }\left({\mathbf{x}}_0\right){\mathbf{F}}^{\prime }{\left({\mathbf{y}}_{\mathsf{\vartheta }}\right)}^{-1}{\mathbf{F}}^{\prime }\left({\mathbf{x}}_0\right){\mathbf{F}}^{\prime }{\left({\mathbf{y}}_{\mathsf{\vartheta }}\right)}^{-1}{\mathbf{F}}^{\prime }\left({\mathbf{x}}_0\right){\mathbf{F}}^{\prime }{\left({\mathbf{y}}_{\mathsf{\vartheta }}\right)}^{-1}\mathbf{F}\left({\mathbf{x}}_0\right)\hfill \\ & +1/2{\mathsf{\alpha }}_0^4{\mathbf{F}}^{\prime }{\left({\mathbf{y}}_{\mathsf{\vartheta }}\right)}^{-1}{\mathbf{F}}^{\prime }\left({\mathbf{x}}_0\right){\mathbf{F}}^{\prime }{\left({\mathbf{y}}_{\mathsf{\vartheta }}\right)}^{-1}{\mathbf{F}}^{″}\left({\mathbf{x}}_0\right){\mathbf{F}}^{\prime }{\left({\mathbf{y}}_{\mathsf{\vartheta }}\right)}^{-1}\mathbf{F}\left({\mathbf{x}}_0\right){\mathbf{F}}^{\prime }{\left({\mathbf{y}}_{\mathsf{\vartheta }}\right)}^{-1}\mathbf{F}\left({\mathbf{x}}_0\right)\hfill \\ & -1/2{\mathsf{\alpha }}_0^4{\mathbf{F}}^{\prime }{\left({\mathbf{y}}_{\mathsf{\vartheta }}\right)}^{-1}{\mathbf{F}}^{″}\left({\mathbf{x}}_0\right){\mathbf{F}}^{\prime }{\left({\mathbf{y}}_{\mathsf{\vartheta }}\right)}^{-1}{\mathbf{F}}^{\prime }\left({\mathbf{x}}_0\right){\mathbf{F}}^{\prime }{\left({\mathbf{y}}_{\mathsf{\vartheta }}\right)}^{-1}\mathbf{F}\left({\mathbf{x}}_0\right){\mathbf{F}}^{\prime }{\left({\mathbf{y}}_{\mathsf{\vartheta }}\right)}^{-1}\mathbf{F}\left({\mathbf{x}}_0\right)\hfill \\ & -1/2{\mathsf{\alpha }}_0^4{\mathbf{F}}^{\prime }{\left({\mathbf{y}}_{\mathsf{\vartheta }}\right)}^{-1}{\mathbf{F}}^{″}\left({\mathbf{x}}_0\right){\mathbf{F}}^{\prime }{\left({\mathbf{y}}_{\mathsf{\vartheta }}\right)}^{-1}\mathbf{F}\left({\mathbf{x}}_0\right){\mathbf{F}}^{\prime }{\left({\mathbf{y}}_{\mathsf{\vartheta }}\right)}^{-1}{\mathbf{F}}^{\prime }\left({\mathbf{x}}_0\right){\mathbf{F}}^{\prime }{\left({\mathbf{y}}_{\mathsf{\vartheta }}\right)}^{-1}\mathbf{F}\left({\mathbf{x}}_0\right)\hfill \\ & -1/6{\mathsf{\alpha }}_0^4{\mathbf{F}}^{\prime }{\left({\mathbf{y}}_{\mathsf{\vartheta }}\right)}^{-1}{\mathbf{F}}^{‴}\left({\mathbf{x}}_0\right){\mathbf{F}}^{\prime }{\left({\mathbf{y}}_{\mathsf{\vartheta }}\right)}^{-1}\mathbf{F}\left({\mathbf{x}}_0\right){\mathbf{F}}^{\prime }{\left({\mathbf{y}}_{\mathsf{\vartheta }}\right)}^{-1}\mathbf{F}\left({\mathbf{x}}_0\right){\mathbf{F}}^{\prime }{\left({\mathbf{y}}_{\mathsf{\vartheta }}\right)}^{-1}\mathbf{F}\left({\mathbf{x}}_0\right).\hfill \end{array}$$

3.1. Iterative Method I

By adding ${\mathbf{x}}_0$ and ${\mathbf{x}}_1$, we get the following iterative method:

$$\begin{array}{c}\hfill \left\{\begin{array}{cc}& {\mathbf{u}}_0=\text{Initial}\text{guess}\hfill \\ & {\mathbf{F}}^{\prime }\left({\mathbf{u}}_0\right){\mathsf{\phi }}_1=\mathbf{F}\left({\mathbf{u}}_0\right)\hfill \\ & {\mathbf{u}}_1={\mathbf{u}}_0-{\mathsf{\phi }}_1\hfill \\ & {\mathbf{F}}^{\prime }\left({\mathbf{u}}_0\right){\mathsf{\phi }}_2=\mathbf{F}\left({\mathbf{u}}_1\right)\hfill \\ & {\mathbf{u}}_2={\mathbf{u}}_1-{\mathsf{\alpha }}_0{\mathsf{\phi }}_2\hfill \end{array}\right.\end{array}$$

(19)

If ${\mathsf{\alpha }}_0=1$, then the order of convergence of Equation (19) is three.

3.2. Iterative Method II

Summing to ${\mathbf{x}}_0$, ${\mathbf{x}}_1$ and ${\mathbf{x}}_2$, we deduce the following iterative method:

$$\begin{array}{c}\hfill \left\{\begin{array}{cc}& {\mathbf{u}}_0=\text{Initial}\text{guess}\hfill \\ & {\mathbf{F}}^{\prime }\left({\mathbf{u}}_0\right){\mathsf{\phi }}_1=\mathbf{F}\left({\mathbf{u}}_0\right)\hfill \\ & {\mathbf{u}}_1={\mathbf{u}}_0-{\mathsf{\phi }}_1\hfill \\ & {\mathbf{F}}^{\prime }\left({\mathbf{u}}_0\right){\mathsf{\phi }}_2=\mathbf{F}\left({\mathbf{u}}_1\right)\hfill \\ & {\mathbf{F}}^{\prime }\left({\mathbf{u}}_0\right){\mathsf{\phi }}_3={\mathbf{F}}^{\prime }\left({\mathbf{u}}_1\right){\mathsf{\phi }}_2\hfill \\ & {\mathbf{u}}_2={\mathbf{u}}_1-2{\mathsf{\phi }}_2+{\mathsf{\phi }}_3\hfill \end{array}\right.\end{array}$$

(20)

The convergence order of Equation (20) is four. Notice that we have the information ${\mathbf{F}}^{\prime }\left({\mathbf{u}}_1\right)$, and we want to use this information in a computationally-efficient way to enhance the convergence order. This means we can add more steps of the form ${\mathbf{F}}^{\prime }\left({\mathbf{u}}_0\right){\mathsf{\phi }}_3={\mathbf{F}}^{\prime }\left({\mathbf{u}}_1\right){\mathsf{\phi }}_2$ in the iterative structure of the iterative scheme. Thus, we can enhance the convergence order of the iterative method Equation (20) by including more parameters owing to the inclusion of additional steps. We construct a model of the iterative method Equation (20) as:

$$\begin{array}{c}\hfill \left\{\begin{array}{cc}& {\mathbf{u}}_0=\text{Initial}\text{guess}\hfill \\ & {\mathbf{F}}^{\prime }\left({\mathbf{u}}_0\right){\mathsf{\phi }}_1=\mathbf{F}\left({\mathbf{u}}_0\right)\hfill \\ & {\mathbf{u}}_1={\mathbf{u}}_0-{\mathsf{\phi }}_1\hfill \\ & {\mathbf{F}}^{\prime }\left({\mathbf{u}}_0\right){\mathsf{\phi }}_2=\mathbf{F}\left({\mathbf{u}}_1\right)\hfill \\ & {\mathbf{F}}^{\prime }\left({\mathbf{u}}_0\right){\mathsf{\phi }}_3={\mathbf{F}}^{\prime }\left({\mathbf{u}}_1\right){\mathsf{\phi }}_2\hfill \\ & {\mathbf{F}}^{\prime }\left({\mathbf{u}}_0\right){\mathsf{\phi }}_4={\mathbf{F}}^{\prime }\left({\mathbf{u}}_1\right){\mathsf{\phi }}_3\hfill \\ & {\mathbf{u}}_2={\mathbf{u}}_1-13/4{\mathsf{\phi }}_2+7/2{\mathsf{\phi }}_3-5/4{\mathsf{\phi }}_4\hfill \end{array}\right.\end{array}$$

(21)

The order of convergence of the iterative method Equation (21) is five.

3.3. Iterative Method III

The summation of ${\mathbf{x}}_0$, ${\mathbf{x}}_1$, ${\mathbf{x}}_2$ and ${\mathbf{x}}_3$ induces the following iterative method:

$$\begin{array}{c}\hfill \left\{\begin{array}{cc}& {\mathbf{u}}_0=\text{Initial}\text{guess}\hfill \\ & {\mathbf{F}}^{\prime }\left({\mathbf{u}}_0\right){\mathsf{\phi }}_1=\mathbf{F}\left({\mathbf{u}}_0\right)\hfill \\ & {\mathbf{u}}_1={\mathbf{u}}_0-{\mathsf{\phi }}_1\hfill \\ & {\mathbf{F}}^{\prime }\left({\mathbf{u}}_0\right){\mathsf{\phi }}_2=\mathbf{F}\left({\mathbf{u}}_1\right)\hfill \\ & {\mathbf{F}}^{\prime }\left({\mathbf{u}}_0\right){\mathsf{\phi }}_3={\mathbf{F}}^{\prime }\left({\mathbf{u}}_1\right){\mathsf{\phi }}_2\hfill \\ & {\mathbf{F}}^{\prime }\left({\mathbf{u}}_0\right){\mathsf{\phi }}_4={\mathbf{F}}^{\prime }\left({\mathbf{u}}_1\right){\mathsf{\phi }}_3\hfill \\ & {\mathbf{F}}^{\prime }\left({\mathbf{u}}_0\right){\mathsf{\phi }}_5={\mathbf{F}}^{″}\left({\mathbf{u}}_1\right){\mathsf{\phi }}_2^2\hfill \\ & {\mathbf{u}}_2={\mathbf{u}}_1+(-2+{\mathsf{\alpha }}_0){\mathsf{\phi }}_2+(1-2{\mathsf{\alpha }}_0){\mathsf{\phi }}_3+{\mathsf{\alpha }}_0{\mathsf{\phi }}_4+(-2{\mathsf{\alpha }}_0-5/2){\mathsf{\phi }}_5\hfill \end{array}\right.\end{array}$$

(22)

The order of convergence of the iterative method Equation (22) is five. We extend the model of iterative method Equation (22), by including more free-parameters, so that we can enhance the convergence order.

$$\begin{array}{c}\hfill \left\{\begin{array}{cc}& {\mathbf{u}}_0=\text{Initial}\text{guess}\hfill \\ & {\mathbf{F}}^{\prime }\left({\mathbf{u}}_0\right){\mathsf{\phi }}_1=\mathbf{F}\left({\mathbf{u}}_0\right)\hfill \\ & {\mathbf{u}}_1={\mathbf{u}}_0-{\mathsf{\phi }}_1\hfill \\ & {\mathbf{F}}^{\prime }\left({\mathbf{u}}_0\right){\mathsf{\phi }}_2=\mathbf{F}\left({\mathbf{u}}_1\right)\hfill \\ & {\mathbf{F}}^{\prime }\left({\mathbf{u}}_0\right){\mathsf{\phi }}_3={\mathbf{F}}^{\prime }\left({\mathbf{u}}_1\right){\mathsf{\phi }}_2\hfill \\ & {\mathbf{F}}^{\prime }\left({\mathbf{u}}_0\right){\mathsf{\phi }}_4={\mathbf{F}}^{\prime }\left({\mathbf{u}}_1\right){\mathsf{\phi }}_3\hfill \\ & {\mathbf{F}}^{\prime }\left({\mathbf{u}}_0\right){\mathsf{\phi }}_5={\mathbf{F}}^{″}\left({\mathbf{u}}_1\right){\mathsf{\phi }}_2^2\hfill \\ & {\mathbf{F}}^{\prime }\left({\mathbf{u}}_0\right){\mathsf{\phi }}_6={\mathbf{F}}^{″}\left({\mathbf{u}}_1\right){\mathsf{\phi }}_2{\mathsf{\phi }}_3\hfill \\ & {\mathbf{u}}_2={\mathbf{u}}_1+{\mathsf{\beta }}_1{\mathsf{\phi }}_2+{\mathsf{\beta }}_2{\mathsf{\phi }}_3+{\mathsf{\beta }}_3{\mathsf{\phi }}_4+{\mathsf{\beta }}_4{\mathsf{\phi }}_5+{\mathsf{\beta }}_5{\mathsf{\phi }}_6\hfill \end{array}\right.\end{array}$$

(23)

If ${\mathsf{\beta }}_1=-3$, ${\mathsf{\beta }}_2=3$, ${\mathsf{\beta }}_3=-1$, ${\mathsf{\beta }}_4=-4$ and ${\mathsf{\beta }}_5=\frac{7}{2}$, then the order of convergence of the iterative method Equation (23) is six. Further iterative methods can be constructed by adding more ${\mathbf{x}}_{\mathrm{i}}^{\prime }\mathrm{s}$ and following essentially the same basic ideas.

4. Application to Systems of Nonlinear Equations Associated with Ordinary and Partial Differential Equations

According to our analysis, we observed that the parametric homotopy method has the potential to provide an effective model for the construction of a higher order iterative method for solving systems of nonlinear equations. The unknown parameters in the homotopy method and the extended models (bases on homotopy provided models) can be used to enhance the convergence order of the iterative method. The iterative method Equation (21) can be employed effectively for solving a general system of nonlinear equations. On the contrary, the iterative method Equation (23) is only efficient to solve weakly nonlinear systems associated with nonlinear ordinary equations (ODEs) and partial differential equations (PDEs). In fact, many ODEs and PDEs can be written as:

$$\begin{array}{c}\hfill \mathrm{L}\left(\mathrm{y}\right)+\mathrm{f}\left(\mathrm{y}\right)-\mathrm{w}\left(\mathrm{x}\right)=0\end{array}$$

(24)

where $\mathrm{L}(·)$ is a linear differential operator, $\mathrm{f}(·)$ is a nonlinear function, y is a dependent variable and x is an independent variable. Let $\mathbf{A}$ be the discrete approximation of linear differential operator $\mathrm{L}(·)$, $\mathrm{f}\left(\mathbf{y}\right)={[\mathrm{f}\left({\mathrm{y}}_1\right),\mathrm{f}\left({\mathrm{y}}_2\right),\cdots ,\mathrm{f}\left({\mathrm{y}}_{\mathrm{n}}\right)]}^{\mathrm{T}}$ and $\mathrm{w}\left(\mathbf{x}\right)={[\mathrm{w}\left({\mathrm{x}}_1\right),\mathrm{w}\left({\mathrm{x}}_2\right),\cdots ,\mathrm{w}\left({\mathrm{x}}_{\mathrm{n}}\right)]}^{\mathrm{T}}$. The nonlinear problem Equation (24) can be written in the form of a system of nonlinear equations:

$$\begin{array}{c}\hfill \mathbf{F}\left(\mathbf{y}\right)=\mathbf{A}\mathbf{y}+\mathrm{f}\left(\mathbf{y}\right)-\mathrm{w}\left(\mathbf{x}\right)=\mathbf{0}\end{array}$$

(25)

It should be noticed that there are discretizations of Equation (24) that do not correspond to Equation (25). Here, we explicitly assume that the numerical methods used for approximating Equation (24) lead to a representation as in Equation (25). The higher order Fréchet derivatives of Equation (25) can be computed as:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill & {\mathbf{F}}^{\prime }\left(\mathbf{y}\right)=\mathbf{A}+\text{diag}\left({\mathrm{f}}^{\prime }\left(\mathbf{y}\right)\right)\hfill \\ & {\mathbf{F}}^{″}\left(\mathbf{y}\right)=\text{diag}\left({\mathrm{f}}^{″}\left(\mathbf{y}\right)\right)\hfill \\ & {\mathbf{F}}^{‴}\left(\mathbf{y}\right)=\text{diag}\left({\mathrm{f}}^{‴}\left(\mathbf{y}\right)\right)\hfill \end{array}\end{array}$$

(26)

where $\text{diag}\left(·\right)$ represents the diagonal matrix, which keeps the input vector as its main diagonal. The iterative method Equation (23) for weakly nonlinear systems can be written as:

$$\begin{array}{c}\hfill \left\{\begin{array}{cc}& {\mathbf{u}}_0=\text{Initial}\text{guess}\hfill \\ & \mathbf{B}=\mathbf{A}+\text{diag}\left({\mathrm{f}}^{\prime }\left({\mathbf{u}}_0\right)\right)\hfill \\ & \mathbf{B}{\mathsf{\phi }}_1=\mathbf{A}{\mathbf{u}}_0+\mathrm{f}\left({\mathbf{u}}_0\right)-\mathrm{w}\left(\mathbf{x}\right)\hfill \\ & {\mathbf{u}}_1={\mathbf{u}}_0-{\mathsf{\phi }}_1\hfill \\ & \mathbf{B}{\mathsf{\phi }}_2=\mathbf{A}{\mathbf{u}}_1+\mathrm{f}\left({\mathbf{u}}_1\right)-\mathrm{w}\left(\mathbf{x}\right)\hfill \\ & \mathbf{B}{\mathsf{\phi }}_3=\mathbf{A}{\mathsf{\phi }}_2+{\mathrm{f}}^{\prime }\left({\mathbf{u}}_1\right)\odot {\mathsf{\phi }}_2\hfill \\ & \mathbf{B}{\mathsf{\phi }}_4=\mathbf{A}{\mathsf{\phi }}_3+{\mathrm{f}}^{\prime }\left({\mathbf{u}}_1\right)\odot {\mathsf{\phi }}_3\hfill \\ & \mathbf{B}{\mathsf{\phi }}_5={\mathrm{f}}^{″}\left({\mathbf{u}}_1\right)\odot {\mathsf{\phi }}_2^2\hfill \\ & \mathbf{B}{\mathsf{\phi }}_6={\mathrm{f}}^{″}\left({\mathbf{u}}_1\right)\odot {\mathsf{\phi }}_2\odot {\mathsf{\phi }}_3\hfill \\ & {\mathbf{u}}_2={\mathbf{u}}_1+{\mathsf{\beta }}_1{\mathsf{\phi }}_2+{\mathsf{\beta }}_2{\mathsf{\phi }}_3+{\mathsf{\beta }}_3{\mathsf{\phi }}_4+{\mathsf{\beta }}_4{\mathsf{\phi }}_5+{\mathsf{\beta }}_5{\mathsf{\phi }}_6\hfill \end{array}\right.\end{array}$$

(27)

where ⊙ is element-wise multiplication between two vectors.

5. Computational Cost of Equation (27)

The iterative method Equation (27) uses five function evaluations, one LU factorization, the solution of six lower and upper triangular systems, four matrix-vector multiplications, six vector-vector multiplications and three scalar-vector multiplications. For weakly nonlinear systems Equation (25), the Fréchet derivatives are diagonal matrices, but can be converted into vectors to perform element-wise vector-vector multiplications. In conclusion, the count of scalar function evaluations for a nonlinear function Equation (25) and its Fréchet derivatives is equal.

6. Convergence Analysis

We will only provide the proof of the convergence of iterative method Equation (23), and the convergence proofs of other iterative methods are similar. To perform the convergence analysis, we use the symbolic algebraic tools of Maple${}^{\circledR }$ 2015.

Theorem 1. 

Let $\mathbf{F}:\Gamma \subseteq {\mathbb{R}}^{\mathrm{n}}\to {\mathbb{R}}^{\mathrm{n}}$ be a sufficiently Fréchet differentiable function on an open convex neighborhood Γ of ${\mathbf{u}}^{*}\in {\mathbb{R}}^{\mathrm{n}}$ with $\mathbf{F}\left({\mathbf{u}}^{*}\right)=0$ and $det\left({\mathbf{F}}^{\prime }\left({\mathbf{u}}^{*}\right)\right)\ne 0$, where ${\mathbf{F}}^{\prime }\left(\mathbf{u}\right)$ denotes the Fréchet derivative of $\mathbf{F}\left(\mathbf{u}\right)$. Let ${\mathbf{D}}_1={\mathbf{F}}^{\prime }\left({\mathbf{u}}^{*}\right)$ and ${\mathbf{D}}_{\mathrm{s}}=\frac{1}{\mathrm{s}!}{{\mathbf{F}}^{\prime }\left({\mathbf{u}}^{*}\right)}^{-1}{\mathbf{F}}^{\left(\mathrm{s}\right)}\left({\mathbf{u}}^{*}\right)$, for $\mathrm{s}\ge 2$, where ${\mathbf{F}}^{\left(\mathrm{s}\right)}\left(\mathbf{u}\right)$ denotes the s-order Fréchet derivative of $\mathbf{F}\left(\mathbf{u}\right)$. Then, the iterative method Equation (23), with an initial guess in the neighborhood of ${\mathbf{u}}^{*}$, converges to ${\mathbf{u}}^{*}$ with local order of convergence at least six and error:

$$\begin{array}{c}\hfill {\mathbf{e}}_2=\mathbf{L}{{\mathbf{e}}_0}^6+\mathrm{O}\left({{\mathbf{e}}_0}^7\right)\end{array}$$

where ${\mathbf{e}}_0={\mathbf{u}}_0-{\mathbf{u}}^{*}$, ${{\mathbf{e}}_0}^{\mathrm{p}}=\stackrel{\mathrm{p}times}{\overbrace{({\mathbf{e}}_0,{\mathbf{e}}_0,\dots ,{\mathbf{e}}_0)}}$ and $\mathbf{L}=\left(2{\mathbf{D}}_2^2{\mathbf{D}}_3{\mathbf{D}}_2-9{\mathbf{D}}_2{\mathbf{D}}_3{\mathbf{D}}_2^2+14{\mathbf{D}}_3{\mathbf{D}}_2^3+16{\mathbf{D}}_2^6-88{\mathbf{D}}_2^3+114{\mathbf{D}}_2^5\right)$ is a six-linear function, i.e., $\mathbf{L}\in \mathbb{L}\stackrel{6times}{\overbrace{\left({\mathbb{R}}^{\mathrm{n}},{\mathbb{R}}^{\mathrm{n}},{\mathbb{R}}^{\mathrm{n}},\cdots ,{\mathbb{R}}^{\mathrm{n}}\right)}}$ with $\mathbf{L}{{\mathbf{e}}_0}^6\in {\mathbb{R}}^{\mathrm{n}}$.

Proof. 

Let $\mathbf{F}:\Gamma \subseteq {\mathbb{R}}^{\mathrm{n}}\to {\mathbb{R}}^{\mathrm{n}}$ be a sufficiently Fréchet differentiable function in Γ. The p-th Fréchet derivative of $\mathbf{F}$ at $\mathrm{v}\in {\mathbb{R}}^{\mathrm{n}}$, $\mathrm{p}\ge 1$, is a p-linear function ${\mathbf{F}}^{\left(\mathrm{p}\right)}\left(\mathrm{v}\right):\stackrel{\mathrm{p}\text{times}}{\overbrace{{\mathbb{R}}^{\mathrm{n}}{\mathbb{R}}^{\mathrm{n}}\cdots {\mathbb{R}}^{\mathrm{n}}}}$ with ${\mathbf{F}}^{\left(\mathrm{p}\right)}\left(\mathrm{v}\right)({\mathrm{h}}_1,{\mathrm{h}}_2,\cdots ,{\mathrm{h}}_{\mathrm{q}})\in {\mathbb{R}}^{\mathrm{n}}$. Taylor’s series expansion of $\mathbf{F}\left({\mathbf{u}}_0\right)$ around ${\mathbf{u}}^{*}$ is:

$$\begin{array}{cc}\hfill \mathbf{F}\left({\mathbf{u}}_0\right)& =\mathbf{F}\left({\mathbf{u}}^{*}+{\mathbf{u}}_0-{\mathbf{u}}^{*}\right)=\mathbf{F}({\mathbf{u}}^{*}+{\mathbf{e}}_0)\hfill \\ & =\mathbf{F}\left({\mathbf{u}}^{*}\right)+{\mathbf{F}}^{\prime }\left({\mathbf{u}}^{*}\right){\mathbf{e}}_0+\frac{1}{2!}{\mathbf{F}}^{″}\left({\mathbf{u}}^{*}\right){\mathbf{e}}_0^2+\frac{1}{3!}{\mathbf{F}}^{\left(3\right)}\left({\mathbf{u}}^{*}\right){\mathbf{e}}_0^3+\frac{1}{4!}{\mathbf{F}}^{\left(4\right)}\left({\mathbf{u}}^{*}\right){\mathbf{e}}_0^4+\cdots \hfill \\ & ={\mathbf{F}}^{\prime }\left({\mathbf{u}}^{*}\right)\left({\mathbf{e}}_0+\frac{1}{2!}{{\mathbf{F}}^{\prime }\left({\mathbf{u}}^{*}\right)}^{-1}{\mathbf{F}}^{\prime \prime }\left({\mathbf{u}}^{*}\right){\mathbf{e}}_0^2+\frac{1}{3!}{{\mathbf{F}}^{\prime }\left({\mathbf{u}}^{*}\right)}^{-1}{\mathbf{F}}^{\left(3\right)}\left({\mathbf{u}}^{*}\right){\mathbf{e}}_0^3+\frac{1}{4!}{{\mathbf{F}}^{\prime }\left({\mathbf{u}}^{*}\right)}^{-1}{\mathbf{F}}^{\left(4\right)}\left({\mathbf{u}}^{*}\right){\mathbf{e}}_0^4+\cdots \right)\hfill \\ & ={\mathbf{D}}_1\left({\mathbf{e}}_0+{\mathbf{D}}_2{\mathbf{e}}_0^2+{\mathbf{D}}_3{\mathbf{e}}_0^3+{\mathbf{D}}_4{\mathbf{e}}_0^4+\mathrm{O}\left({{\mathbf{e}}_0}^5\right)\right)\hfill \end{array}$$

where ${\mathbf{e}}_{\mathrm{k}}={\mathbf{u}}_{\mathrm{k}}-{\mathbf{u}}^{*}$. Computing the Fréchet derivative of $\mathbf{F}$ with respect to ${\mathbf{e}}_0$, we deduce:

$$\begin{array}{c}\hfill {\mathbf{F}}^{\prime }\left({\mathbf{u}}_0\right)={\mathbf{D}}_1\left(\mathbf{I}+2{\mathbf{D}}_2{\mathbf{e}}_0+3{\mathbf{D}}_3{{\mathbf{e}}_0}^2+4{\mathbf{D}}_4{{\mathbf{e}}_0}^3+\mathrm{O}\left({{\mathbf{e}}_0}^4\right)\right)\end{array}$$

(28)

where $\mathbf{I}$ is the identity matrix. Computing its inverse using the symbolic mathematical package Maple, we obtain:

$$\begin{array}{cc}\hfill {{\mathbf{F}}^{\prime }\left({\mathbf{u}}_0\right)}^{-1}& =(\mathbf{I}-2{\mathbf{D}}_2{\mathbf{e}}_0+\left(4{\mathbf{D}}_2^2-3{\mathbf{D}}_3\right){\mathbf{e}}_0^2+\left(6{\mathbf{D}}_3{\mathbf{D}}_2+6{\mathbf{D}}_2{\mathbf{D}}_3-8{\mathbf{D}}_2^3-4{\mathbf{D}}_4\right){\mathbf{e}}_0^3+(8{\mathbf{D}}_4{\mathbf{D}}_2+9{\mathbf{D}}_3^2\hfill \\ & +8{\mathbf{D}}_2{\mathbf{D}}_4-5{\mathbf{D}}_5-12{\mathbf{D}}_3{\mathbf{D}}_2^2-12{\mathbf{D}}_2{\mathbf{D}}_3{\mathbf{D}}_2-12{\mathbf{D}}_2^2{\mathbf{D}}_3+16{\mathbf{D}}_2^4){\mathbf{e}}_0^4+\mathrm{O}\left({{\mathbf{e}}_0}^5\right)){\mathbf{D}}_1^{-1}\hfill \end{array}$$

(29)

By using Maple, we find the following expressions:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\mathsf{\phi }}_1& ={\mathbf{e}}_0-{\mathbf{D}}_2{\mathbf{e}}_0^2+\left(2{\mathbf{D}}_2^2-2{\mathbf{D}}_3\right){\mathbf{e}}_0^3+\left(-4{\mathbf{D}}_2^3-3{\mathbf{D}}_4+4{\mathbf{D}}_2{\mathbf{D}}_3+3{\mathbf{D}}_3{\mathbf{D}}_2\right){\mathbf{e}}_0^4+(-6{\mathbf{D}}_3{\mathbf{D}}_2^2\hfill \\ & -6{\mathbf{D}}_2{\mathbf{D}}_3{\mathbf{D}}_2+8{\mathbf{D}}_2^4-8{\mathbf{D}}_2^2{\mathbf{D}}_3-4{\mathbf{D}}_5+6{\mathbf{D}}_2{\mathbf{D}}_4+6{\mathbf{D}}_3^2+4{\mathbf{D}}_4{\mathbf{D}}_2){\mathbf{e}}_0^5+(12{\mathbf{D}}_2^2{\mathbf{D}}_3{\mathbf{D}}_2\hfill \\ & +12{\mathbf{D}}_2{\mathbf{D}}_3{\mathbf{D}}_2^2+12{\mathbf{D}}_3{\mathbf{D}}_2^3-8{\mathbf{D}}_2{\mathbf{D}}_4{\mathbf{D}}_2-9{\mathbf{D}}_3^2{\mathbf{D}}_2-8{\mathbf{D}}_4{\mathbf{D}}_2^2-12{\mathbf{D}}_2{\mathbf{D}}_3^2-12{\mathbf{D}}_3{\mathbf{D}}_2{\mathbf{D}}_3\hfill \\ & -12{\mathbf{D}}_2^2{\mathbf{D}}_4+16{\mathbf{D}}_2^3{\mathbf{D}}_3-16{\mathbf{D}}_2^5+{\mathbf{D}}_6+8{\mathbf{D}}_2{\mathbf{D}}_5+9{\mathbf{D}}_3{\mathbf{D}}_4+8{\mathbf{D}}_4{\mathbf{D}}_3+5{\mathbf{D}}_5{\mathbf{D}}_2){\mathbf{e}}_0^6+\mathrm{O}\left({\mathbf{e}}_0^7\right)\hfill \end{array}\end{array}$$

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\mathbf{e}}_1& ={\mathbf{D}}_2{\mathbf{e}}_0^2+\left(-2{\mathbf{D}}_2^2+2{\mathbf{D}}_3\right){\mathbf{e}}_0^3+\left(4{\mathbf{D}}_2^3+3{\mathbf{D}}_4-4{\mathbf{D}}_2{\mathbf{D}}_3-3{\mathbf{D}}_3{\mathbf{D}}_2\right){\mathbf{e}}_0^4+(6{\mathbf{D}}_3{\mathbf{D}}_2^2+6{\mathbf{D}}_2{\mathbf{D}}_3{\mathbf{D}}_2-8{\mathbf{D}}_2^4\hfill \\ & +8{\mathbf{D}}_2^2{\mathbf{D}}_3+4{\mathbf{D}}_5-6{\mathbf{D}}_2{\mathbf{D}}_4-6{\mathbf{D}}_3^2-4{\mathbf{D}}_4{\mathbf{D}}_2){\mathbf{e}}_0^5+(-12{\mathbf{D}}_2^2{\mathbf{D}}_3{\mathbf{D}}_2-12{\mathbf{D}}_2{\mathbf{D}}_3{\mathbf{D}}_2^2-12{\mathbf{D}}_3{\mathbf{D}}_2^3\hfill \\ & +8{\mathbf{D}}_2{\mathbf{D}}_4{\mathbf{D}}_2+9{\mathbf{D}}_3^2{\mathbf{D}}_2+8{\mathbf{D}}_4{\mathbf{D}}_2^2+12{\mathbf{D}}_2{\mathbf{D}}_3^2+12{\mathbf{D}}_3{\mathbf{D}}_2{\mathbf{D}}_3+12{\mathbf{D}}_2^2{\mathbf{D}}_4-16{\mathbf{D}}_2^3{\mathbf{D}}_3+16{\mathbf{D}}_2^5-{\mathbf{D}}_6\hfill \\ & -8{\mathbf{D}}_2{\mathbf{D}}_5-9{\mathbf{D}}_3{\mathbf{D}}_4-8{\mathbf{D}}_4{\mathbf{D}}_3-5{\mathbf{D}}_5{\mathbf{D}}_2){\mathbf{e}}_0^6+\mathrm{O}\left({\mathbf{e}}_0^7\right)\hfill \end{array}\end{array}$$

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \mathbf{F}\left({\mathbf{u}}_1\right)& ={\mathbf{D}}_1{\mathbf{D}}_2{\mathbf{e}}_0^2+\left(2{\mathbf{D}}_1{\mathbf{D}}_3-2{\mathbf{D}}_1{\mathbf{D}}_2^2\right){\mathbf{e}}_0^3+\left(5{\mathbf{D}}_1{\mathbf{D}}_2^3+3{\mathbf{D}}_1{\mathbf{D}}_4-4{\mathbf{D}}_1{\mathbf{D}}_2{\mathbf{D}}_3-3{\mathbf{D}}_1{\mathbf{D}}_3{\mathbf{D}}_2\right){\mathbf{e}}_0^4+\hfill \\ & (6{\mathbf{D}}_1{\mathbf{D}}_3{\mathbf{D}}_2^2+8{\mathbf{D}}_1{\mathbf{D}}_2{\mathbf{D}}_3{\mathbf{D}}_2-12{\mathbf{D}}_1{\mathbf{D}}_2^4+10{\mathbf{D}}_1{\mathbf{D}}_2^2{\mathbf{D}}_3+4{\mathbf{D}}_1{\mathbf{D}}_5-6{\mathbf{D}}_1{\mathbf{D}}_2{\mathbf{D}}_4-6{\mathbf{D}}_1{\mathbf{D}}_3^2\hfill \\ & -4{\mathbf{D}}_1{\mathbf{D}}_4{\mathbf{D}}_2){\mathbf{e}}_0^5+(28{\mathbf{D}}_1{\mathbf{D}}_2^5+11{\mathbf{D}}_1{\mathbf{D}}_2{\mathbf{D}}_4{\mathbf{D}}_2-19{\mathbf{D}}_1{\mathbf{D}}_2^2{\mathbf{D}}_3{\mathbf{D}}_2-19{\mathbf{D}}_1{\mathbf{D}}_2{\mathbf{D}}_3{\mathbf{D}}_2^2+16{\mathbf{D}}_1{\mathbf{D}}_2{\mathbf{D}}_3^2\hfill \\ & -24{\mathbf{D}}_1{\mathbf{D}}_2^3{\mathbf{D}}_3+15{\mathbf{D}}_1{\mathbf{D}}_2^2{\mathbf{D}}_4-8{\mathbf{D}}_1{\mathbf{D}}_2{\mathbf{D}}_5-9{\mathbf{D}}_1{\mathbf{D}}_3{\mathbf{D}}_4-8{\mathbf{D}}_1{\mathbf{D}}_4{\mathbf{D}}_3-5{\mathbf{D}}_1{\mathbf{D}}_5{\mathbf{D}}_2+8{\mathbf{D}}_1{\mathbf{D}}_4{\mathbf{D}}_2^2\hfill \\ & +9{\mathbf{D}}_1{\mathbf{D}}_3^2{\mathbf{D}}_2+12{\mathbf{D}}_1{\mathbf{D}}_3{\mathbf{D}}_2{\mathbf{D}}_3-11{\mathbf{D}}_1{\mathbf{D}}_3{\mathbf{D}}_2^3-{\mathbf{D}}_1{\mathbf{D}}_6){\mathbf{e}}_0^6+\mathrm{O}\left({\mathbf{e}}_0^7\right)\hfill \end{array}\end{array}$$

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\mathsf{\phi }}_2& ={\mathbf{D}}_2{\mathbf{e}}_0^2+\left(-4{\mathbf{D}}_2^2+2{\mathbf{D}}_3\right){\mathbf{e}}_0^3+\left(-6{\mathbf{D}}_3{\mathbf{D}}_2+13{\mathbf{D}}_2^3-8{\mathbf{D}}_2{\mathbf{D}}_3+3{\mathbf{D}}_4\right){\mathbf{e}}_0^4+(-38{\mathbf{D}}_2^4-8{\mathbf{D}}_4{\mathbf{D}}_2\hfill \\ & +20{\mathbf{D}}_2{\mathbf{D}}_3{\mathbf{D}}_2+18{\mathbf{D}}_3{\mathbf{D}}_2^2-12{\mathbf{D}}_3^2+26{\mathbf{D}}_2^2{\mathbf{D}}_3-12{\mathbf{D}}_2{\mathbf{D}}_4+4{\mathbf{D}}_5){\mathbf{e}}_0^5+(-59{\mathbf{D}}_2^2{\mathbf{D}}_3{\mathbf{D}}_2-55{\mathbf{D}}_2{\mathbf{D}}_3{\mathbf{D}}_2^2\hfill \\ & -50{\mathbf{D}}_3{\mathbf{D}}_2^3+27{\mathbf{D}}_2{\mathbf{D}}_4{\mathbf{D}}_2+27{\mathbf{D}}_3^2{\mathbf{D}}_2+24{\mathbf{D}}_4{\mathbf{D}}_2^2+40{\mathbf{D}}_2{\mathbf{D}}_3^2+36{\mathbf{D}}_3{\mathbf{D}}_2{\mathbf{D}}_3+39{\mathbf{D}}_2^2{\mathbf{D}}_4-76{\mathbf{D}}_2^3{\mathbf{D}}_3\hfill \\ & +88{\mathbf{D}}_2^3+16{\mathbf{D}}_2^5-{\mathbf{D}}_6-16{\mathbf{D}}_2{\mathbf{D}}_5-18{\mathbf{D}}_3{\mathbf{D}}_4-16{\mathbf{D}}_4{\mathbf{D}}_3-10{\mathbf{D}}_5{\mathbf{D}}_2){\mathbf{e}}_0^6+\mathrm{O}\left({\mathbf{e}}_0^7\right)\hfill \end{array}\end{array}$$

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\mathsf{\phi }}_3& ={\mathbf{D}}_2{\mathbf{e}}_0^2+\left(-6{\mathbf{D}}_2^2+2{\mathbf{D}}_3\right){\mathbf{e}}_0^3+\left(-9{\mathbf{D}}_3{\mathbf{D}}_2+27{\mathbf{D}}_2^3-12{\mathbf{D}}_2{\mathbf{D}}_3+3{\mathbf{D}}_4\right){\mathbf{e}}_0^4+(42{\mathbf{D}}_2{\mathbf{D}}_3{\mathbf{D}}_2-104{\mathbf{D}}_2^4\hfill \\ & -18{\mathbf{D}}_3^2+36{\mathbf{D}}_3{\mathbf{D}}_2^2+54{\mathbf{D}}_2^2{\mathbf{D}}_3-12{\mathbf{D}}_4{\mathbf{D}}_2-18{\mathbf{D}}_2{\mathbf{D}}_4+4{\mathbf{D}}_5){\mathbf{e}}_0^5+(-163{\mathbf{D}}_2^2{\mathbf{D}}_3{\mathbf{D}}_2-149{\mathbf{D}}_2{\mathbf{D}}_3{\mathbf{D}}_2^2\hfill \\ & -128{\mathbf{D}}_3{\mathbf{D}}_2^3+57{\mathbf{D}}_2{\mathbf{D}}_4{\mathbf{D}}_2+54{\mathbf{D}}_3^2{\mathbf{D}}_2+48{\mathbf{D}}_4{\mathbf{D}}_2^2+84{\mathbf{D}}_2{\mathbf{D}}_3^2+72{\mathbf{D}}_3{\mathbf{D}}_2{\mathbf{D}}_3+16{\mathbf{D}}_2^6+81{\mathbf{D}}_2^2{\mathbf{D}}_4\hfill \\ & -208{\mathbf{D}}_2^3{\mathbf{D}}_3+88{\mathbf{D}}_2^3+258{\mathbf{D}}_2^5-{\mathbf{D}}_6-24{\mathbf{D}}_2{\mathbf{D}}_5-27{\mathbf{D}}_3{\mathbf{D}}_4-24{\mathbf{D}}_4{\mathbf{D}}_3-15{\mathbf{D}}_5{\mathbf{D}}_2){\mathbf{e}}_0^6+\mathrm{O}\left({\mathbf{e}}_0^7\right)\hfill \end{array}\end{array}$$

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\mathsf{\phi }}_4& ={\mathbf{D}}_2{\mathbf{e}}_0^2+\left(-8{\mathbf{D}}_2^2+2{\mathbf{D}}_3\right){\mathbf{e}}_0^3+\left(3{\mathbf{D}}_4+45{\mathbf{D}}_2^3-16{\mathbf{D}}_2{\mathbf{D}}_3-12{\mathbf{D}}_3{\mathbf{D}}_2\right){\mathbf{e}}_0^4+(-210{\mathbf{D}}_2^4+90{\mathbf{D}}_2^2{\mathbf{D}}_3\hfill \\ & -24{\mathbf{D}}_2{\mathbf{D}}_4-24{\mathbf{D}}_3^2-16{\mathbf{D}}_4{\mathbf{D}}_2+60{\mathbf{D}}_3{\mathbf{D}}_2^2+70{\mathbf{D}}_2{\mathbf{D}}_3{\mathbf{D}}_2+4{\mathbf{D}}_5){\mathbf{e}}_0^5+(-{\mathbf{D}}_6-329{\mathbf{D}}_2^2{\mathbf{D}}_3{\mathbf{D}}_2\hfill \\ & -299{\mathbf{D}}_2{\mathbf{D}}_3{\mathbf{D}}_2^2-260{\mathbf{D}}_3{\mathbf{D}}_2^3+95{\mathbf{D}}_2{\mathbf{D}}_4{\mathbf{D}}_2+90{\mathbf{D}}_3^2{\mathbf{D}}_2+80{\mathbf{D}}_4{\mathbf{D}}_2^2+140{\mathbf{D}}_2{\mathbf{D}}_3^2+120{\mathbf{D}}_3{\mathbf{D}}_2{\mathbf{D}}_3\hfill \\ & +135{\mathbf{D}}_2^2{\mathbf{D}}_4-420{\mathbf{D}}_2^3{\mathbf{D}}_3+748{\mathbf{D}}_2^5-32{\mathbf{D}}_2{\mathbf{D}}_5-36{\mathbf{D}}_3{\mathbf{D}}_4-32{\mathbf{D}}_4{\mathbf{D}}_3-20{\mathbf{D}}_5{\mathbf{D}}_2+88{\mathbf{D}}_2^3+32{\mathbf{D}}_2^6)\hfill \\ & {\mathbf{e}}_0^6+\mathrm{O}\left({\mathbf{e}}_0^7\right)\hfill \end{array}\end{array}$$

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\mathsf{\phi }}_5& =2{\mathbf{D}}_2^3{\mathbf{e}}_0^4+\left(4{\mathbf{D}}_2{\mathbf{D}}_3{\mathbf{D}}_2-20{\mathbf{D}}_2^4+4{\mathbf{D}}_2^2{\mathbf{D}}_3\right){\mathbf{e}}_0^5+(124{\mathbf{D}}_2^5+6{\mathbf{D}}_2{\mathbf{D}}_4{\mathbf{D}}_2-36{\mathbf{D}}_2^2{\mathbf{D}}_3{\mathbf{D}}_2-28{\mathbf{D}}_2{\mathbf{D}}_3{\mathbf{D}}_2^2\hfill \\ & +8{\mathbf{D}}_2{\mathbf{D}}_3^2-40{\mathbf{D}}_2^3{\mathbf{D}}_3+6{\mathbf{D}}_2^2{\mathbf{D}}_4){\mathbf{e}}_0^6+\mathrm{O}\left({\mathbf{e}}_0^7\right)\hfill \end{array}\end{array}$$

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\mathsf{\phi }}_6& =2{\mathbf{D}}_2^3{\mathbf{e}}_0^4+\left(4{\mathbf{D}}_2{\mathbf{D}}_3{\mathbf{D}}_2-24{\mathbf{D}}_2^4+4{\mathbf{D}}_2^2{\mathbf{D}}_3\right){\mathbf{e}}_0^5+(8{\mathbf{D}}_2{\mathbf{D}}_3^2-36{\mathbf{D}}_2{\mathbf{D}}_3{\mathbf{D}}_2^2-48{\mathbf{D}}_2^3{\mathbf{D}}_3+176{\mathbf{D}}_2^5\hfill \\ & -42{\mathbf{D}}_2^2{\mathbf{D}}_3{\mathbf{D}}_2+6{\mathbf{D}}_2{\mathbf{D}}_4{\mathbf{D}}_2+6{\mathbf{D}}_2^2{\mathbf{D}}_4){\mathbf{e}}_0^6+\mathrm{O}\left({\mathbf{e}}_0^7\right)\hfill \end{array}\end{array}$$

By using all of the values of ϕ’s, we get the following error equation.

$$\begin{array}{c}\hfill {\mathbf{e}}_2=\left(2{\mathbf{D}}_2^2{\mathbf{D}}_3{\mathbf{D}}_2-9{\mathbf{D}}_2{\mathbf{D}}_3{\mathbf{D}}_2^2+14{\mathbf{D}}_3{\mathbf{D}}_2^3+16{\mathbf{D}}_2^6-88{\mathbf{D}}_2^3+114{\mathbf{D}}_2^5\right){\mathbf{e}}_0^6+\mathrm{O}\left({\mathbf{e}}_0^7\right)\end{array}$$

(30)

The error Equation (30) directly implies that the local order of convergence of the iterative method Equation (23) is at least six. ☐

It is beneficial to construct the multi-step version of the iterative method Equation (27) because we get an increment of factor two in the convergence order attained in the previous step, at the cost of single function evaluation plus the solution of two lower and upper triangular systems. The multi-step version of Equation (27) can be written as:

$$\begin{array}{c}\hfill \left\{\begin{array}{cc}& {\mathbf{u}}_0=\text{Initial}\text{guess}\hfill \\ & \mathbf{B}=\mathbf{A}+\text{diag}\left({\mathrm{f}}^{\prime }\left({\mathbf{u}}_0\right)\right)\hfill \\ & \mathbf{B}{\mathsf{\phi }}_1=\mathbf{A}{\mathbf{u}}_0+\mathrm{f}\left({\mathbf{u}}_0\right)-\mathrm{w}\left(\mathbf{x}\right)\hfill \\ & {\mathbf{u}}_1={\mathbf{u}}_0-{\mathsf{\phi }}_1\hfill \\ & \mathbf{B}{\mathsf{\phi }}_2=\mathbf{A}{\mathbf{u}}_1+\mathrm{f}\left({\mathbf{u}}_1\right)-\mathrm{w}\left(\mathbf{x}\right)\hfill \\ & \mathbf{B}{\mathsf{\phi }}_3=\mathbf{A}{\mathsf{\phi }}_2+{\mathrm{f}}^{\prime }\left({\mathbf{u}}_1\right)\odot {\mathsf{\phi }}_2\hfill \\ & \mathbf{B}{\mathsf{\phi }}_4=\mathbf{A}{\mathsf{\phi }}_3+{\mathrm{f}}^{\prime }\left({\mathbf{u}}_1\right)\odot {\mathsf{\phi }}_3\hfill \\ & \mathbf{B}{\mathsf{\phi }}_5={\mathrm{f}}^{″}\left({\mathbf{u}}_1\right)\odot {\mathsf{\phi }}_2^2\hfill \\ & \mathbf{B}{\mathsf{\phi }}_6={\mathrm{f}}^{″}\left({\mathbf{u}}_1\right)\odot {\mathsf{\phi }}_2\odot {\mathsf{\phi }}_3\hfill \\ & {\mathbf{u}}_2={\mathbf{u}}_1+{\mathsf{\beta }}_1{\mathsf{\phi }}_2+{\mathsf{\beta }}_2{\mathsf{\phi }}_3+{\mathsf{\beta }}_3{\mathsf{\phi }}_4+{\mathsf{\beta }}_4{\mathsf{\phi }}_5+{\mathsf{\beta }}_5{\mathsf{\phi }}_6\hfill \\ & \text{for}\mathrm{i}=3,\mathrm{m}\hfill \\ & \mathbf{B}{\mathsf{\psi }}_1=\mathbf{A}{\mathbf{u}}_{\mathrm{i}-1}+\mathrm{f}\left({\mathbf{u}}_{\mathrm{i}-1}\right)-\mathrm{w}\left(\mathbf{x}\right)\hfill \\ & \mathbf{B}{\mathsf{\psi }}_2={\mathrm{f}}^{″}\left({\mathbf{u}}_1\right)\odot {\mathsf{\phi }}_1\odot {\mathsf{\psi }}_1\hfill \\ & {\mathbf{u}}_{\mathrm{i}}={\mathbf{u}}_{\mathrm{i}-1}-{\mathsf{\psi }}_1-{\mathsf{\psi }}_2\hfill \\ & \text{end}\hfill \end{array}\right.\end{array}$$

(31)

We claim that the convergence order of the m-step iterative method Equation (31) is $2(\mathrm{m}+1)$.

Theorem 2. 

Let $\mathbf{F}:\Gamma \subseteq {\mathbb{R}}^{\mathrm{n}}\to {\mathbb{R}}^{\mathrm{n}}$ be a sufficiently Fréchet differentiable function on an open convex neighborhood Γ of ${\mathbf{u}}^{*}\in {\mathbb{R}}^{\mathrm{n}}$ with $\mathbf{F}\left({\mathbf{u}}^{*}\right)=0$ and $det\left({\mathbf{F}}^{\prime }\left({\mathbf{u}}^{*}\right)\right)\ne 0$, where ${\mathbf{F}}^{\prime }\left(\mathbf{u}\right)$ denotes the Fréchet derivative of $\mathbf{F}\left(\mathbf{u}\right)$. Let ${\mathbf{D}}_1={\mathbf{F}}^{\prime }\left({\mathbf{u}}^{*}\right)$ and ${\mathbf{D}}_{\mathrm{s}}=\frac{1}{\mathrm{s}!}{{\mathbf{F}}^{\prime }\left({\mathbf{u}}^{*}\right)}^{-1}{\mathbf{F}}^{\left(\mathrm{s}\right)}\left({\mathbf{u}}^{*}\right)$, for $\mathrm{s}\ge 2$, where ${\mathbf{F}}^{\left(\mathrm{s}\right)}\left(\mathbf{u}\right)$ denotes the s-order Fréchet derivative of $\mathbf{F}\left(\mathbf{u}\right)$. Then, the iterative method Equation (31), with an initial guess in the neighborhood of ${\mathbf{u}}^{*}$, converges to ${\mathbf{u}}^{*}$ with local order of convergence at least $2(\mathrm{m}+1)$ and error:

$$\begin{array}{c}\hfill {\mathbf{e}}_{\mathrm{m}}=\mathbf{L}{{\mathbf{e}}_0}^{2(\mathrm{m}+1)}+\mathrm{O}\left({{\mathbf{e}}_0}^{2(\mathrm{m}+1)+1}\right)\end{array}$$

where ${\mathbf{e}}_0={\mathbf{u}}_0-{\mathbf{u}}^{*}$, ${{\mathbf{e}}_0}^{\mathrm{p}}=\stackrel{\mathrm{p}times}{\overbrace{({\mathbf{e}}_0,{\mathbf{e}}_0,\dots ,{\mathbf{e}}_0)}}$ and $\mathbf{L}={\left(3\left(2{\mathbf{D}}_2^2+{\mathbf{D}}_3\right)\right)}^{\mathrm{m}-2}\left(42{\mathbf{D}}_2^5+2{\mathbf{D}}_2^2{\mathbf{D}}_3{\mathbf{D}}_2-9{\mathbf{D}}_2{\mathbf{D}}_3{\mathbf{D}}_2^2+14{\mathbf{D}}_3{\mathbf{D}}_2^3\right)$ is a $2(\mathrm{m}+1)$-linear function, i.e., $\mathbf{L}\in \mathbb{L}\stackrel{2(\mathrm{m}+1)times}{\overbrace{\left({\mathbb{R}}^{\mathrm{n}},{\mathbb{R}}^{\mathrm{n}},{\mathbb{R}}^{\mathrm{n}},\cdots ,{\mathbb{R}}^{\mathrm{n}}\right)}}$ with $\mathbf{L}{{\mathbf{e}}_0}^6\in {\mathbb{R}}^{\mathrm{n}}$.

Proof. 

We will prove the convergence order via mathematical induction. Suppose that the claimed order of convergence is true for $\mathrm{m}=\mathrm{s}$, i.e., ${\mathbf{e}}_{\mathrm{s}}=\mathrm{O}\left({\mathbf{e}}^{2(\mathrm{s}+1)}\right)$. For the rest of the proof, we proceed as follows:

$$\begin{array}{c}\hfill \begin{array}{cc}& {\mathbf{u}}_{\mathrm{s}}-{\mathbf{u}}^{*}\sim {\mathbf{e}}_{\mathrm{s}}\hfill \\ & \mathbf{F}\left({\mathrm{u}}_{\mathrm{s}}\right)\sim {\mathbf{D}}_1{\mathbf{e}}_{\mathrm{s}}\hfill \\ & {\mathbf{B}}^{-1}\sim \left(\mathbf{I}-2{\mathbf{D}}_2{\mathbf{e}}_0\right){\mathbf{D}}_1^{-1}\hfill \\ & {\mathsf{\psi }}_1\sim \left(\mathbf{I}-2{\mathbf{D}}_2{\mathbf{e}}_0\right){\mathbf{e}}_{\mathrm{s}}={\mathbf{e}}_{\mathrm{s}}-2{\mathbf{D}}_2{\mathbf{e}}_0{\mathbf{e}}_{\mathrm{s}}\hfill \\ & {\mathbf{B}}^{-1}{\mathbf{F}}^{″}\left({\mathbf{u}}_1\right){\mathsf{\phi }}_1\sim 2{\mathbf{D}}_2{\mathbf{e}}_0-6{\mathbf{D}}_2^2{\mathbf{e}}^2\hfill \\ & {\mathsf{\psi }}_2\sim \left(2{\mathbf{D}}_2{\mathbf{e}}_0-6{\mathbf{D}}_2^2{\mathbf{e}}_0^2\right)\left(\mathbf{I}-2{\mathbf{D}}_2{\mathbf{e}}_0\right){\mathbf{e}}_{\mathrm{s}}\sim 2{\mathbf{D}}_2{\mathbf{e}}_0{\mathbf{e}}_{\mathrm{s}}-10{\mathbf{D}}_2^2{\mathbf{e}}_0^2{\mathbf{e}}_{\mathrm{s}}\hfill \\ & {\mathbf{u}}_{\mathrm{s}+1}-{\mathbf{u}}^{*}={\mathbf{u}}_{\mathrm{s}}-{\mathbf{u}}^{*}-{\mathsf{\psi }}_1-{\mathsf{\psi }}_2\sim {\mathbf{e}}_{\mathrm{s}}-{\mathbf{e}}_{\mathrm{s}}+2{\mathbf{D}}_2{\mathbf{e}}_0{\mathbf{e}}_{\mathrm{s}}-2{\mathbf{D}}_2{\mathbf{e}}_0{\mathbf{e}}_{\mathrm{s}}+10{\mathbf{D}}_2^2{\mathbf{e}}_0^2{\mathbf{e}}_{\mathrm{s}}\hfill \\ & {\mathbf{u}}_{\mathrm{s}+1}-{\mathbf{u}}^{*}\sim {\mathbf{e}}_0^2{\mathbf{e}}_{\mathrm{s}}\hfill \\ & {\mathbf{e}}_{\mathrm{s}+1}=\mathrm{O}\left({\mathbf{e}}_0^{2(\mathrm{s}+2)}\right)\hfill \end{array}\end{array}$$

(32)

☐

7. Numerical Testing

To show the effectiveness of our proposed multi-step iterative method, we solve a carefully-chosen set of ODEs and PDEs. It is also important to verify the claimed order of convergence numerically, so we adopt the following definition to compute the computational order of convergence (COC):

$$\begin{array}{c}\hfill \text{COC}=\frac{\mathrm{l}\mathrm{o}\mathrm{g}\left(\left|\right|\mathbf{F}\left({\mathbf{x}}_{\mathrm{k}+1}\right){\left|\right|}_{\infty }/\left|\right|\mathbf{F}\left({\mathbf{x}}_{\mathrm{k}}\right){\left|\right|}_{\infty }\right)}{\mathrm{l}\mathrm{o}\mathrm{g}\left(\left|\right|\mathbf{F}\left({\mathbf{x}}_{\mathrm{k}}\right){\left|\right|}_{\infty }/\left|\right|\mathbf{F}\left({\mathbf{x}}_{\mathrm{k}-1}\right){\left|\right|}_{\infty }\right)}\end{array}$$

Consider a system of four nonlinear equations:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\mathrm{F}}_1\left(\mathbf{x}\right)& ={\mathrm{x}}_2{\mathrm{x}}_3+{\mathrm{x}}_4\left({\mathrm{x}}_2+{\mathrm{x}}_3\right)=0\hfill \\ \hfill {\mathrm{F}}_2\left(\mathbf{x}\right)& ={\mathrm{x}}_1{\mathrm{x}}_3+{\mathrm{x}}_4\left({\mathrm{x}}_1+{\mathrm{x}}_3\right)=0\hfill \\ \hfill {\mathrm{F}}_3\left(\mathbf{x}\right)& ={\mathrm{x}}_1{\mathrm{x}}_2+{\mathrm{x}}_4\left({\mathrm{x}}_1+{\mathrm{x}}_2\right)=0\hfill \\ \hfill {\mathrm{F}}_4\left(\mathbf{x}\right)& ={\mathrm{x}}_1{\mathrm{x}}_2+{\mathrm{x}}_3\left({\mathrm{x}}_1+{\mathrm{x}}_2\right)-1=0\hfill \end{array}\end{array}$$

(33)

Let $\mathbf{w}={[{\mathrm{w}}_1,{\mathrm{w}}_2,{\mathrm{w}}_3,{\mathrm{w}}_4]}^{\mathrm{T}}$ be a constant vector, and ${\mathbf{F}}^{\prime }\left(\mathbf{x}\right)$ and ${\mathbf{F}}^{″}\left(\mathbf{x}\right)\mathbf{w}={\left({\mathbf{F}}^{\prime }\left(\mathbf{x}\right)\mathbf{w}\right)}^{\prime }$ can be written as:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\mathbf{F}}^{\prime }\left(\mathbf{x}\right)& =\left[\begin{array}{cccc}0& {\mathrm{x}}_3+{\mathrm{x}}_4& {\mathrm{x}}_2+{\mathrm{x}}_4& {\mathrm{x}}_2+{\mathrm{x}}_3\\ {\mathrm{x}}_3+{\mathrm{x}}_4& 0& {\mathrm{x}}_1+{\mathrm{x}}_4& {\mathrm{x}}_1+{\mathrm{x}}_3\\ {\mathrm{x}}_2+{\mathrm{x}}_4& {\mathrm{x}}_1+{\mathrm{x}}_4& 0& {\mathrm{x}}_1+{\mathrm{x}}_2\\ {\mathrm{x}}_2+{\mathrm{x}}_3& {\mathrm{x}}_1+{\mathrm{x}}_3& {\mathrm{x}}_1+{\mathrm{x}}_2& 0\end{array}\right]\hfill \\ \hfill {\mathbf{F}}^{″}\left(\mathbf{x}\right)\mathbf{w}& =\left[\begin{array}{cccc}0& {\mathrm{w}}_3+{\mathrm{w}}_4& {\mathrm{w}}_2+{\mathrm{w}}_4& {\mathrm{w}}_2+{\mathrm{w}}_3\\ {\mathrm{w}}_3+{\mathrm{w}}_4& 0& {\mathrm{w}}_1+{\mathrm{w}}_4& {\mathrm{w}}_1+{\mathrm{w}}_3\\ {\mathrm{w}}_2+{\mathrm{w}}_4& {\mathrm{w}}_1+{\mathrm{w}}_4& 0& {\mathrm{w}}_1+{\mathrm{w}}_2\\ {\mathrm{w}}_2+{\mathrm{w}}_3& {\mathrm{w}}_1+{\mathrm{w}}_3& {\mathrm{w}}_1+{\mathrm{w}}_2& 0\end{array}\right]\hfill \end{array}\end{array}$$

(34)

In

Table 1. Verification of the convergence order of different iterative methods for the problem Equation (33); the initial guess is ${\mathbf{x}}_{\mathrm{k}}=15/10$. COC, computational order of convergence.

Iterations \  Methods	(19)	(20)	(21)	(22)
1      $\left|\right|\mathbf{F}\left({\mathbf{x}}_{\mathrm{k}}\right){\left|\right|}_{\infty }$	8.88e-01	5.80e-01	4.12e-01	4.12e-01
2       -	3.57e-02	2.48e-03	9.94e-05	9.94e-05
3       -	1.33e-06	6.41e-14	5.51e-25	5.51e-25
4       -	7.985e-21	4.48e-58	4.63e-129	4.63e-129
5       -	1.91e-64	1.67e-236	3.09e-652	3.09e-652
6       -	2.90e-196	4.99e-952	6.59e-3271	6.59e-3271
7       -	1.13e-592	6.26e-3816	4.63e-16,367	4.63e-16,367
8       -	7.53e-1783	2.43e-15,273	1.27e-81,850	1.27e−81,850
COC	3.00	4.00	5.00	5.00

, we computed the COC’s of different iterative methods, and numerical results confirm the theoretically-claimed convergence orders. The last two columns of Table 1 are similar because the structure of iterative methods Equations (21) and (22) is similar, and iterative method Equation (22) has an extra parameter that does not enhance the convergence order.

Table 2. Multi-step iterative method Equation (31): verification of convergence order for the problem Equation (33).

Iterations\ Steps	$\mathrm{m}=2$	$\mathrm{m}=3$	$\mathrm{m}=4$	$\mathrm{m}=5$
1      $\left|\right|\mathbf{F}\left({\mathbf{x}}_{\mathrm{k}}\right){\left|\right|}_{\infty }$	6.65e-02	1.66e-02	4.69e-03	1.28e-03
2       -	6.10e-12	5.46e-21	1.24e-32	7.41e-47
3       -	1.42e-75	3.71e-175	2.15e-337	1.38e-577
4       -	2.85e-461	8.08e-1415	5.24e-3394	3.16e-6958
COC	6.06	8.04	10.0	12.0

shows the COC for different values of m for the iterative method Equation (31), and the COC satisfies the formula $2(\mathrm{m}+1)$ that was claimed. The numerical computations in Table 1 and Table 2 are performed in Mathematica${}^{\circledR }$ 10 by setting $MinPrecision = 82,000 and $MinPrecision = 7000, respectively.

7.1. 3D Nonlinear Poisson Problem

We consider the following nonlinear Poisson–Dirichlet boundary value problem:

$$\begin{array}{c}\hfill {\mathrm{u}}_{\mathrm{x}\mathrm{x}}+{\mathrm{u}}_{\mathrm{y}\mathrm{y}}+{\mathrm{u}}_{\mathrm{z}\mathrm{z}}+\mathrm{f}\left(\mathrm{u}\right)=\mathrm{p}(\mathrm{x},\mathrm{y},\mathrm{z}),(\mathrm{x},\mathrm{y},\mathrm{z})\in {(0,1)}^3\end{array}$$

(35)

where $\mathrm{f}(·)$ is a nonlinear function of $\mathrm{u}(\mathrm{x},\mathrm{y},\mathrm{z})$ and $\mathrm{p}(\mathrm{x},\mathrm{y},\mathrm{z})$ is the source term. The discretization of Equation (35) is performed by using the Chebyshev pseudo-spectral collocation method [23,24,25,26,27]. The discretized form of Equation (35) is:

$$\begin{array}{c}\hfill \begin{array}{cc}& \mathbf{F}\left(\mathbf{U}\right)=\left(\left({\mathbf{T}}_{\mathrm{x}\mathrm{x}}\otimes {\mathbf{I}}_{\mathrm{y}}\otimes {\mathbf{I}}_{\mathrm{z}}\right)+\left({\mathbf{I}}_{\mathrm{x}}\otimes {\mathbf{T}}_{\mathrm{y}\mathrm{y}}\otimes {\mathbf{I}}_{\mathrm{z}}\right)+\left({\mathbf{I}}_{\mathrm{x}}\otimes {\mathbf{I}}_{\mathrm{y}}\otimes {\mathbf{T}}_{\mathrm{z}\mathrm{z}}\right)\right)\mathbf{U}+\mathrm{f}\left(\mathbf{U}\right)-\mathbf{p}=\mathbf{0}\hfill \\ & {\mathbf{F}}^{\prime }\left(\mathbf{U}\right)=\left(\left({\mathbf{T}}_{\mathrm{x}\mathrm{x}}\otimes {\mathbf{I}}_{\mathrm{y}}\otimes {\mathbf{I}}_{\mathrm{z}}\right)+\left({\mathbf{I}}_{\mathrm{x}}\otimes {\mathbf{T}}_{\mathrm{y}\mathrm{y}}\otimes {\mathbf{I}}_{\mathrm{z}}\right)+\left({\mathbf{I}}_{\mathrm{x}}\otimes {\mathbf{I}}_{\mathrm{y}}\otimes {\mathbf{T}}_{\mathrm{z}\mathrm{z}}\right)\right)+\text{diag}\left({\mathrm{f}}^{\prime }\left(\mathbf{U}\right)\right)\hfill \\ & {\mathbf{F}}^{″}\left(\mathbf{U}\right)=\text{diag}\left({\mathrm{f}}^{″}\left(\mathbf{U}\right)\right)\hfill \end{array}\end{array}$$

(36)

where ${\mathbf{T}}_{··}$ is the Chebyshev pseudo-spectral collocation discretization of the second order derivative, $\mathbf{I}$ is the identity matrix and ⊗ is a Kronecker product. We assume the solution of Equation (35) is $\mathrm{u}=\mathrm{s}\mathrm{i}\mathrm{n}(\mathrm{x}+\mathrm{y}+\mathrm{z})$ and $\mathrm{f}\left(\mathrm{u}\right)={\mathrm{u}}^2$. In

Table 3. Multi-step iterative method Equation (31): absolute error in the solution of Equation (36) versus different grid sizes. Number of iterations = 1; initial guess $\mathbf{U}=\mathbf{0}$.

m \N	$8\times 8\times 8$	$10\times 10\times 10$	$12\times 12\times 12$
2      $\left|\right|\mathbf{F}\left({\mathbf{U}}_{\mathrm{k}}\right){\left|\right|}_{\infty }$	6.70e-07	7.51e-07	9.07e-07
3       -	1.54e-09	1.42e-09	1.47e-09
4       -	5.62e-10	1.09e-11	1.12e-11
5       -	5.62e-10	3.20e-13	7.45e-14

, we show the absolute error in the solution of Equation (36). In all of the cases, we perform a single iteration and get maximum reduction in the absolute error when $\mathrm{m}=5$ and the grid size is $12\times 12\times 12$. In one iteration, we perform the LU factorization only once, and then, we solve lower and upper triangular systems for linear subproblems, which make the procedure computationally inexpensive.

7.2. 2D Nonlinear Wave Equation

The following 2D nonlinear wave equation is studied.

$$\begin{array}{c}\hfill {\mathrm{u}}_{\mathrm{t}\mathrm{t}}-{\mathrm{c}}^2\left({\mathrm{u}}_{\mathrm{x}\mathrm{x}}+{\mathrm{u}}_{\mathrm{y}\mathrm{y}}\right)+\mathrm{f}\left(\mathrm{u}\right)=\mathrm{p}(\mathrm{x},\mathrm{y}),(\mathrm{x},\mathrm{y},\mathrm{t})\in {(-1,1)}^2\times (0,2]\end{array}$$

(37)

where $\mathrm{f}\left(\mathrm{u}\right)={\mathrm{u}}^{\mathrm{s}}$ is nonlinear function and c, s are constants. The source term $\mathrm{p}(\mathrm{x},\mathrm{y})$ can be computed by assuming solution $\mathrm{u}=\mathrm{e}\mathrm{x}\mathrm{p}(-\mathrm{t})\mathrm{s}\mathrm{i}\mathrm{n}(\mathrm{x}+\mathrm{y})$. Dirichlet boundary conditions are imposed. We use the Chebyshev pseudo-spectral collocation method to discretize Equation (37), and we obtain the following system of nonlinear equations.

$$\begin{array}{c}\hfill \begin{array}{cc}& \mathbf{F}\left(\mathbf{U}\right)=\left(\left({\mathbf{T}}_{\mathrm{t}\mathrm{t}}\otimes {\mathbf{I}}_{\mathrm{x}}\otimes {\mathbf{I}}_{\mathrm{y}}\right)-{\mathrm{c}}^2\left(\left({\mathbf{I}}_{\mathrm{t}}\otimes {\mathbf{T}}_{\mathrm{x}\mathrm{x}}\otimes {\mathbf{I}}_{\mathrm{y}}\right)+\left({\mathbf{I}}_{\mathrm{t}}\otimes {\mathbf{I}}_{\mathrm{x}}\otimes {\mathbf{T}}_{\mathrm{y}\mathrm{y}}\right)\right)\right)\mathbf{U}+\mathrm{f}\left(\mathbf{U}\right)-\mathbf{p}=\mathbf{0}\hfill \\ & {\mathbf{F}}^{\prime }\left(\mathbf{U}\right)=\left({\mathbf{T}}_{\mathrm{t}\mathrm{t}}\otimes {\mathbf{I}}_{\mathrm{x}}\otimes {\mathbf{I}}_{\mathrm{y}}\right)-{\mathrm{c}}^2\left(\left({\mathbf{I}}_{\mathrm{t}}\otimes {\mathbf{T}}_{\mathrm{x}\mathrm{x}}\otimes {\mathbf{I}}_{\mathrm{y}}\right)+\left({\mathbf{I}}_{\mathrm{t}}\otimes {\mathbf{I}}_{\mathrm{x}}\otimes {\mathbf{T}}_{\mathrm{y}\mathrm{y}}\right)\right)+\text{diag}\left({\mathrm{f}}^{\prime }\left(\mathbf{U}\right)\right)\hfill \\ & {\mathbf{F}}^{″}\left(\mathbf{U}\right)=\text{diag}\left({\mathrm{f}}^{″}\left(\mathbf{U}\right)\right)\hfill \end{array}\end{array}$$

(38)

Table 4. Multi-step iterative method Equation (31): absolute error in the solution of Equation (38) versus different grid sizes. Number of iterations = 1, $\mathrm{s}=3$, $\mathrm{c}=1$; initial guess $\mathbf{U}=\mathbf{0}$.

m \N	$10\times 10\times 10$	$20\times 20\times 20$	$20\times 20\times 30$
2      $\left|\right|\mathbf{F}\left({\mathbf{U}}_{\mathrm{k}}\right){\left|\right|}_{\infty }$	2.82e-04	8.34e-04	5.46e-05
3       -	2.81e-05	3.42e-05	2.99e-05
4       -	8.32e-06	5.99e-06	1.66e-05
5       -	8.23e-07	2.02e-06	9.22e-06
6       -		8.71e-07	5.12e-06
9       -		5.53e-08	8.77e-07
12      -		3.45e-09	1.50e-07
13      -		1.37e-09	8.35e-08
14      -		5.40e-10	4.64e-08
16      -		8.75e-11	1.43e-08
17      -			7.94e-09
21      -			7.55e-10
25      -			7.19e-11

shows that the number of steps, to get a reduction in the absolute error in the solution of Equation (38), depends on the step-size in the spectral approximation.

7.3. Verification of the Mesh-Independence Principle

In the paper [28], the authors proved the mesh-independence principle for the classical Newton method in order to solve systems of nonlinear equations associated with a general class of nonlinear boundary value problems. They also have stated some conditions on the discretization of the nonlinear boundary value problems. In simple words, suppose we have two discretizations (one is finer than the other) of a nonlinear boundary value problem, and suppose that we use the Newton method to solve two associated systems of nonlinear equations, one for each discretization. The magnitude of absolute errors for both discretizations against the number of Newton iterations will be almost equal. Table 4 and

Table 5. Multi-step Newton method: absolute error in the solution of Equation (38) versus different grid sizes. Number of iterations = 1, $\mathrm{s}=3$, $\mathrm{c}=1$; initial guess $\mathbf{U}=\mathbf{0}$.

m \N	$10\times 10\times 10$	$20\times 20\times 20$	$20\times 20\times 30$
2      $\left|\right|\mathbf{F}\left({\mathbf{U}}_{\mathrm{k}}\right){\left|\right|}_{\infty }$	6.14e-03	2.13e-02	3.69e-03
3       -	1.23e-03	1.79e-03	1.70-04
4       -	2.11e-04	7.85e-04	6.29e-05
5       -	4.68e-05	4.41e-04	1.02e-05
6       -	2.09e-06	2.15e-04	5.12e-06
9       -	2.60e-07	2.93e-05	7.04e-07
12      -		3.90e-06	4.91e-08
13      -		1.99e-06	2.03e-08
14      -		1.02e-06	8.36e-09
16      -		2.66e-07	1.42e-09
17      -		1.36e-07	5.88e-10
21      -		9.28e-09	2.14e-11
25      -		6.36e-10	2.13e-11
28      -		8.77e-11	2.12e-11

show that multi-step iterative methods do not obey the mesh-independence principle, when we perform a single evaluation of the Jacobian and its inversion.

Table 6. Newton method and iterative method Equation (27): absolute error in the solution of Equation (38) versus different grid sizes. Number of iterations = 4; number of steps = 0, $\mathrm{s}=3$, $\mathrm{c}=1$; initial guess $\mathbf{U}=\mathbf{0}$.

Iterations \N	$\mathbf{10}\mathbf{\times }\mathbf{10}\mathbf{\times }\mathbf{10}$	$\mathbf{15}\mathbf{\times }\mathbf{15}\mathbf{\times }\mathbf{15}$	$\mathbf{20}\mathbf{\times }\mathbf{20}\mathbf{\times }\mathbf{20}$
Newton method
1      $\left|\right|\mathbf{F}\left({\mathbf{U}}_{\mathrm{k}}\right){\left|\right|}_{\infty }$	5.78e-02	4.52e-02	5.13e-02
2       -	1.66e-03	2.70e-04	2.24e-04
3       -	7.84e-07	2.43e-08	2.16e-09
4       -	1.94e-07	4.74e-12	5.39e-11
Iterations \N	$30\times 20\times 20$	$20\times 30\times 20$	$20\times 20\times 30$
1      $\left|\right|\mathbf{F}\left({\mathbf{U}}_{\mathrm{k}}\right){\left|\right|}_{\infty }$	5.11e-02	5.11-02	5.19e-02
2       -	2.24e-04	2.24e-04	2.23e-04
3       -	2.55e-09	2.55e-09	2.65e-09
4       -	5.36e-11	5.37e-11	2.13e-11
Iterations \N	$10\times 10\times 10$	$15\times 15\times 15$	$20\times 20\times 20$
Iterative method Equation (27)
1      $\left|\right|\mathbf{F}\left({\mathbf{U}}_{\mathrm{k}}\right){\left|\right|}_{\infty }$	2.82e-04	8.42e-05	8.34e-04
2       -	1.94e-07	4.74e-12	5.39e-11
3       -	1.94e-07	4.74e-12	5.39e-11
4       -	1.94e-07	4.74e-12	5.39e-11
Iterations \N	$30\times 20\times 20$	$20\times 30\times 20$	$20\times 20\times 30$
1      $\left|\right|\mathbf{F}\left({\mathbf{U}}_{\mathrm{k}}\right){\left|\right|}_{\infty }$	5.23e-04	5.23e-04	5.46e-05
2       -	5.37e-11	5.37e-11	2.14e-11
3       -	5.37e-11	5.37e-11	2.14e-11
4       -	5.37e-11	5.37e-11	2.14e-11

displays that the Newton method and the iterative method Equation (27) almost verify the mesh-independence principle. Notice that the iterative method in Equation (31) is the multi-step version of the iterative method Equation (27). The multi-step iterative methods for the problem Equation (35) does not satisfy the mesh-independence property for the problem Equation (35), while Newton and the iterative method in Equation (27) have this very nice feature.

Remark. In a single iteration of multi-step iterative methods, we use the frozen Jacobian at the initial guess, and we employ its LU factors repeatedly in the multi-step part to solve the related lower and upper triangular systems. The idea of the frozen Jacobian is fruitful because the computation of a new Jacobian and its LU factorization are computationally expensive operations. Now, we ask whether the applicability of such methods is limited. What we have observed in numerical simulations is the following: when one has to solve a stiff initial or boundary value problem, in the start of the simulations, it is recommended to perform more than one iteration with a limited number of multi-steps. In fact, in the case of stiff problems, the computation of the Jacobian repeatedly is crucial for the convergence of the whole iterative process.

8. Conclusions

The homotopy methods provide a way to construct models for the development of higher order iterative methods, to solve systems of nonlinear equations. We have shown that the inclusion of parameters in the homotopy can help to enhance the convergence order of the proposed techniques. Usually, we get higher order Fréchet derivatives that are computationally unsuitable for a general system of nonlinear equations. However, for particular classes of ODEs and PDEs, we have shown that the computational cost of higher order Fréchet derivatives is not high. The numerical results show the validity and effectiveness of the proposed multi-step iterative methods.