A multi-step Ulm-Chebyshev-like method for solving nonlinear operator equations

1.   Introduction

Let $X$ and $Y$ be Banach spaces, $D\in X$ be an open subset, and $F$: $D\in X\rightarrow Y$ be a nonlinear operator with the continuous Fr$\acute{{\rm e}}$chet derivative denoted by $F'$. We consider the problem of approximating a solution ${\bf{x}}^{*}$ of a nonlinear equation

which applies to the inverse eigenvalue problems ^[1,2,3,4], the generalized inverse eigenvalue problems ^[5], the inverse singular value problems ^[6,7,8], the Chandrasekhar integral equation ^[9], the neutral differential-algebraic equations ^[10,11], and so on. Without any doubt, Newton's method is the most-used iterative process to solve this problem. It is given by the algorithm:

for ${\bf{x}}_{0}$ given. This iterative process has quadratic R-order of convergence under some mild conditions ^[12,13,14]. For recent progress on Newton's method, one may refer to ^[15,16,17].

Other methods, such as higher-order methods, also include in their expression the inverse of the operator $F'$. To avoid this problem, Newton-type methods:

where $H_{n}$ is an approximation of $F'({\bf{x}}_{n})^{-1}$ are considered. One of these methods was proposed by Moser in ^[18]. Given ${\bf{x}}_{0}\in D$ and $B_{0}\in \mathcal{L}(Y, X)$, the following sequences are defined:

The first equation is similar to Newton's method, but replaces the operator $F'({\bf{x}}_{n})^{-1}$ with a linear operator $B_{n}$. The second equation is Newton's method applied to the equation

where $g_{n}$: $\mathcal{L}(Y, X)\rightarrow \mathcal{L}(X, Y)$ is defined by

So $\{B_{n}\}$ gives us an approximation of $F'({\bf{x}}_{n})^{-1}$. It can be shown that the rate of convergence for the above scheme is $(1+\sqrt{5})/2$, provided the root of (1.1) is simple ^[18]. However, from a numerical perspective, this is unsatisfactory because the scheme uses the same amount of information in each step as Newton's method, but its convergence speed is not faster than the secant method. For that, in ^[19], Ulm proposed the following iterative method to solve nonlinear equations. Given ${\bf{x}}_{0}\in D$ and $B_{0}\in \mathcal{L}(Y, X)$, Ulm defines

Notice that, here $F'({\bf{x}}_{n+1})$ appears instead of $F'({\bf{x}}_{n})$ in (1.2). This is crucial for obtaining fast convergence. Under the classical assumption that the derivative $F'$ is Lipschitz continuous around the solution, Ulm showed that the method generates successive approximations that converge to a solution of (1.1) asymptotically as fast as Newton's method. For recent progress on Newton-Moser type method ^[20,21], one may refer to Moser's method ^[22], Ulm's method ^{[23,24,25,26]}, and Ulm-like method ^[27,28].

Considering the previous antecedent, in order to extend the above ideas to high-order convergent iterative methods (cubic convergence), in ^[29], Ezquerro and Hernández considered Chebyshev's method and proposed the following iterative method to solve nonlinear equations. Given ${\bf{x}}_{0}\in D$ and $B_{0}\in \mathcal{L}(Y, X)$, Ulm-Chebyshev defines

which does not use any inverse operator in its application. Ezquerro and Hernández showed that the method generates successive approximations that converge to a solution of (1.1), and has cubical convergence. Recently, some authors have employed Ulm-Chebyshev's method to solve inverse eigenvalue problems and inverse singular value problems ^[1,2,5,7]. There, they found that computing exactly the derivative $F'({\bf{x}}_{n})\ (n = 0, 1, 2, \ldots)$ at each iteration is costly, especially in the case when the system is large.

In order to reduce the cost of accurately calculating the derivative $F'({\bf{x}}_{n})\ (n = 0, 1, 2, \ldots)$ in each iteration, motivated by the Ulm and Ulm-Chebyshev methods, we propose a multi-step Ulm-Chebyshev-like method for solving the nonlinear operator equation

Given ${\bf{x}}_{0}\in D$ and $B_{0}\in \mathcal{L}(Y, X)$, the multi-step Ulm-Chebyshev-like method is defined by

where $A_{n+1}$ is an approximation of the derivative $F'({\bf{x}}_{n+1})$. This method exhibits several attractive features. First, it is inverse free: we do not need to solve a linear equation at each iteration. Second, it is derivative free: we do not need to computer the Fr$\acute{{\rm e}}$chet derivative at each iteration. Third, in addition to solving the nonlinear Eq (1.1), the method produces successive approximations $\{B_{n}\}$ to the value of $F'({\bf{x}}^{*})^{-1}$, having ${\bf{x}}^{*}$ as a solution of (1.1). This property is very helpful, especially when one investigates the sensitivity of the solution to small perturbations. Fourth, the method converges to the solution with $R$-convergence rate 4.

Further more, in Section 2, we analyze the local convergence of the new iterative method. Under certain assumptions, the radius of the convergence ball for the multi-step Ulm-Chebyshev-like method is estimated, and the $R$-convergence rate 4 of the multi-step Ulm-Chebyshev-like method is proved. Section 3 is devoted to showing some of the most important features of the new iterative method by means of five examples. Compared with other existing methods, the proposed method has higher convergence order and/or requires less operations.

2.   Convergence analysis

Let $\mathbf{B}({\bf{x}}, r)$ stand for the open ball in $X$ with center ${\bf{x}}$ and radius $r > 0$. Let ${\bf{x}}^{*}\in D$ be a solution of the nonlinear Eq (1.1) such that $F'({\bf{x}}^{*})$ is invertible and that $F'$ satisfies the Lipschitz condition on $\mathbf{B}({\bf{x}}^{*}, r)$ with the Lipschitz constant $L$:

Let $\{ {\bf{x}}_{n}\}$ be generated by the multi-step Ulm-Chebyshev-like method. Let $A_{n}$ be an approximation of the derivative $F'({\bf{x}}_{n})$ such that

where $\{\eta_{n}\}$ is a nonnegative-valued sequence satisfying $\sup_{n\geq 0}\eta_{n}\leq\eta$ where $\eta$ is a nonnegative constant. Let

The following lemma is crucial for the proof of the main theorem.

Lemma 2.1. (^[27]) If ${\bf{x}}_{n}\in \mathbf{B}({\bf{x}}, r_{L})$, then the following assertions hold:

and $A_{n}$ is invertible and $\|A_{n}^{-1}\|\leq\mu$.

Note that in the multi-step Ulm-Chebyshev-like method, sequence $\{B_{n}\}$ is generated by the algorithm except for $B_{0}$. Below, we prove that if $B_{0}$ approximates $A_{0}^{-1}$, then the sequence $\{ {\bf{x}}_{n}\}$ generated by the multi-step Ulm-Chebyshev-like method converges locally to ${\bf{x}}^{*}$ with $R$-convergence rate 4. For this end, let $B_{0}$ satisfy that

where $\xi$ is defined in (2.4).

Theorem 2.1. Suppose that the Jacobian matrix $F'({\bf{x}}^{*})$ is invertible and that $F'$ satisfies the Lipschitz condition (2.1) on $\mathbf{B}({\bf{x}}^{*}, r_{L})$. Then there exist positive numbers $\beta$ and $\xi$ such that for any ${\bf{x}}_{0}\in\mathbf{B}({\bf{x}}^{*}, \beta)$ and $B_{0}$ satisfying (2.5), the sequence $\{ {\bf{x}}_{n}\}$ generated by the multi-step Ulm-Chebyshev-like method with initial point ${\bf{x}}_{0}$ converges to ${\bf{x}}^{*}$. Moreover, the following estimates hold for each $n = 0, 1, \ldots.$

and

where $\alpha$ and $\beta$ are defined in (2.4).

Proof. We proceed by mathematical induction. Clearly, (2.6) is trivial for $n = 0$ by the assumption. By (2.4) and (2.5), we obtain

That is, (2.7) holds for $n = 0$. Now we assume that (2.6) and (2.7) hold for $n = m$. Then one has

and

By (2.4), we get

It follows from (2.4), (2.8), and Lemma 2.1 that

and

Then

and

By (1.5), we have

where

for $0\leq\theta\leq1$. Since

and

it follows from (2.4), (2.8), (2.11), (2.12), and the Lipschitz condition that

which together with (2.4), (2.8), and $\alpha\leq1$ gives

It follows from (2.11), (2.12), and the Lipschitz condition that

Similar to (2.13)–(2.15) and by (2.4) and $\alpha\leq1$, we have

and

By (2.4), (2.16), (2.18), and $\alpha\leq1$, we get

Consequently, (2.6) holds for $n = m+1$, and by (2.4), (2.8), and (2.19), we get

By (2.4), we obtain

and it follows from (2.4), (2.8), and Lemma 2.1 that

Together with (2.1), (2.4), (2.10), and (2.20), we have

which follows from (2.9) and (2.11), and we get

From the fourth equation in (1.5), we obtain

which together with (2.23) gives

Notice that

and we have

It follows from (2.4), (2.24), and $\alpha\leq1$ that

This confirms that (2.7) holds for $n = m+1$ and the proof is complete.□

Remark 2.1. Under the conditions as in Theorem 2.2, the sequence ${\bf{x}}^{k}$ converges to the limit ${\bf{x}}^{*}$ with $R$-convergence rate 4.

3.   Numerical experiments

In this section, we report the numerical performance of the multi-step Ulm-Chebyshev-like method for solving the nonlinear operator Eq (1.1). We compare the multi-step Ulm-Chebyshev-like method (Algorithm MSUCL) with the Newton-type method (Algorithm NT) in ^[16], the Chebyshev-like method (Algorithm CL) in ^[30], the Ulm-like method (Algorithm UL) in ^[27], and the Ulm-Chebyshev method (Algorithm UC) in ^[29]. All tests were carried out in MALAB 7.10 running on a PC Intel Pentium IV with a 3.0 GHz CPU.

Example 3.1. (^[30]) We consider the following system of 3 nonlinear equations:

The Jacobian is given by

and

correct to $14$ decimal places in this case. We choose the starting vector

for $n = 0, 1, \ldots, \eta_{n} = \frac{1}{10}$. For all algorithms, the stopping tolerance for the iterations is $10^{-12}$.

From Table 1, we observe that the Algorithm MSMCL converges to the solution with $R$-convergence rate 4, the Algorithm UL converges quadratically, and the Algorithm NT, the Algorithm UL, and the Algorithm UC converge cubically in the root sense.

Example 3.2. (^[27]) We next consider the two-point boundary value problem

We divide the interval $[0, 1]$ into $m+1$ subintervals and we get

Let $d_{0}, d_{1}, \ldots, d_{m+1}$ be the points of subdivision with

An approximation for the second derivative may be chosen as

Let the operator $\phi$: $\mathbf{R}^{m}\rightarrow \mathbf{R}^{m}$ be defined by

To get an approximation to the solution of (3.1), we need to solve the following nonlinear equation:

where

Obviously, ${\bf{x}}^{*} = 0$ is a solution of (3.2) and

Hence

Furthermore, it is easy to verify that

where $\|\cdot\|$ denotes the F-norm. For different choices of $m$ and ${\bf{x}}_{0}$, the convergence performance of the algorithm is illustrated in the following tables. Here we consider the following three problem sizes:

(a) $m = 10$ and ${\bf{x}}_{0} = \sigma(1, 1, \ldots, 1)^{T}$;

(b) $m = 100$ and ${\bf{x}}_{0} = \sigma(1, 1, \ldots, 1)^{T}$;

(c) $m = 1000$ and ${\bf{x}}_{0} = \sigma(1, 1, \ldots, 1)^{T}$, where $\sigma = 0.2\ {\rm or}\ 0.02$.

For all algorithms, the stopping tolerance for the iterations is $10^{-12}$.

Tables 2–4 show that the new method has higher convergence order and from Table 5, we see that the CPU time by the Algorithm MSUCL is less than the other existing methods.

Remark 3.1. Various indices can be employed to compare the efficiency of iterative methods. One of the most prevalent efficiency indices is introduced as follows ^[31,32]:

where OC, IT, FE, JE, and CPU are the order of convergence, number of iterations, number of total function evaluations, total number of Jacobian evaluations, and CPU time, respectively.

We can see from Table 6 that Algorithm MSUCL is better than the other methods for solving all test problems in criteria CR2, indicating its superior performance.

3.1.   Chandrasekhar H-equation

The Chandrasekhar integral equation ^[9] which arises from radiative transfer theory is a nonlinear integral equation which gives a full nonlinear system of equations if discretized. The Chandrasekhar integral equation is given by

with parameter $c$ and the operator $F$ as

If we discretize the integral Eq (3.3) using the midpoint integration rule with $n$ grid points:

we obtain the resulting system of nonlinear equations:

When starting with a $(1, 1, \cdots, 1)^{T}$ vector, the system (3.4) has a solution for all $c\in(0, \ 1)$. The $c$ were equally spaced with $\Delta c = 0.01$ in the interval $c\in(0, \ 1)$ and we choose $n = 100$. We note that in this case the Jacobian is a full matrix for $n = 0, 1, \ldots, \eta_{n} = \frac{1}{10}$. For all algorithms, the stopping tolerance for the iterations is $10^{-12}$.

Let $Iter_{Total}$ denote the total number of iterations for all $c$ considered and $Iter$ be its mean iteration number. From Table 7, we find that the new method has the least total and mean number of iterations, which has the lowest computational cost and therefore is the most efficient of the other existing methods in terms of CPU time.

3.2.   Inverse eigenvalue problem

We consider the following inverse eigenvalue problem (IEP): given $n+1$ real symmetric $n\times n$ matrices $\{A_{i}\}_{i = 0}^{n}$ and $n$ real numbers

find a vector

such that

where

and $\{\lambda_{j}(A({\bf c}))\}_{j = 1}^n$ are the eigenvalues of $A({\bf c})$ with

The above IEP can be represented mathematically through a set of non-linear equations:

To further illustrate the effectiveness of the new algorithm, we present a practical engineering application in vibrations ^[2,33,34]. We consider the vibration of a taut string with $n$ beads. Figure 1 shows such a model for the case $n = 4$.

Here, we assume that the $n$ beads are placed along the string, where the ends of the string are clamped. The mass of the $j$th bead is denoted by $m_j$. The horizontal lengths between masses $m_j$ and $m_{j+1}$ (and between each beads at each end and the clamped support) are set to be a constant $L$. The horizontal tension is set to be a constant $T$. Then the equation of motion is governed by

where

That is, the ends of the string are fixed. The matrix form of (3.6) is given by:

where

with

and $J$ is the discrete Laplacian matrix

The general solution of (3.7) is given in terms of the eigenvalue problem

where $\lambda$ is the square of the natural frequency of the vibration system and the nonzero vector $y$ accounts for the interplay between the masses. The inverse problem for the beaded string is to compute the masses $\{m_j\}_{j = 1}^n$ so that the resulting system has a prescribed set of natural frequencies. It is easy to check that the eigenvalues of $J$ are given by:

Thus $J$ is symmetric and positive definite and $CJ$ is similar to $L^TCL$, where $L$ is the Cholesky factor of

Then the inverse problem is converted into the form of the IEP where

with

for $j = 1, 2, \ldots, n$. The beaded string data in Example 3.3 comes from the website: http://www.caam.rice.edu/beads.

Example 3.3. ^[2]) This is an inverse problem for the beaded string with $n = 6$ beads, where

We report our numerical results for starting point

We solve Example 3.3 using all algorithms, and the stopping tolerance is set to be

for $n = 0, 1, \ldots, \eta_{n} = \frac{1}{10}$. The numerical results are listed in Tables 8 and 9. We observe from Table 8 that the proposed method has higher convergence order and/or requires less operations than the other existing methods. Table 9 displays the computed masses for the beaded string. As expected, the desired masses are recovered: http://www.caam.rice.edu/beads.

Example 3.4. (^[1]) This is an inverse problem with $n = 10$. Define

where

Now

Then

We report our numerical results for different starting points:

● (a) ${\bf c}^{0} = (-77.95824, -62.08697, 96.54128, 40.10535, -44.33137, 20.79310)^T$,

● (b) ${\bf c}^{0} = (-76.86213, -63.46336, 95.28928, 41.39452, -42.24157, 17.37889)^T$,

● (c) ${\bf c}^{0} = (-78.58345, -65.97678, 97.83621, 43.47844, -49.26789, 23.67335)^T$,

● (d) ${\bf c}^{0} = (-85.47863, -67.28566, 80.28746, 35.38552, -45.45096, 23.47528)^T$.

Here we take

For all algorithms, the stopping tolerance is set to be

Table 10 displays the error of $\| {\bf c}^{k}- {\bf c}^{*}\|$ for the above four initial points ${\bf c}^{0}$, where "It." represents the number of outer iterations, and "*" denotes that the corresponding algorithm fails to converge, respectively.

We see from Table 10 that for these choices of the initial points, Algorithm UL, Algorithm UC, and Algorithm MSUCL converge but Algorithm NT and Algorithm CL do not, and Algorithm MSUCL needs less iterations than Algorithm UL and Algorithm UC.

4.   Conclusions

In this paper, we have proposed a multi-step Ulm-Chebyshev-like method for solving nonlinear operator equations, which does not contain inverse operators in its expression, and does not require computing Jacobian matrices for solving Jacobian equations. We prove that the multi-step Ulm-Chebyshev-like method converges locally to the solution with $R$-convergence rate 4 under appropriate conditions. As an application, it is demonstrated how this result can be used to analyze the Chandrasekhar integral equation and to solve the inverse eigenvalue problems. The proposed method has higher convergence order and/or requires less operations than the other existing methods, indicating its superior performance. The study of the stability analysis of our new method is also one for our future work.

Author contributions

Wei Ma: algorithms, software, numerical examples, writing–original draft, writing–review and editing; Ming Zhao: algorithms, software, numerical examples, writing–original draft, writing–review and editing; Jiaxin Li: algorithms, software, numerical examples, writing–original draft, writing–review and editing. All authors of this article have been contributed equally. All authors have read and approved the final version of the manuscript for publication.

Acknowledgments

This work was supported by the Henan Province College Student Innovation and Entrepreneurship Training Program Project (No. 202410481003).

Conflict of interest

The authors declare no conflicts of interest.