Fractional Newton-Raphson Method and Some Variants for the Solution of Non-linear Systems

TL;DR

This paper introduces novel fractional derivative-based numerical methods for solving non-linear systems in complex space, improving convergence properties and addressing discontinuities present in traditional fractional Newton-Raphson methods.

Contribution

It proposes new fractional Newton-Raphson variants with variable derivative orders, enhancing convergence and stability for complex non-linear systems.

Findings

Methods achieve at least quadratic convergence.

New approaches avoid discontinuities in fractional derivatives.

Applicable to systems with real initial conditions.

Abstract

The following document presents some novel numerical methods valid for one and several variables, which using the fractional derivative, allow to find solutions for some non-linear systems in the complex space using real initial conditions. The origin of these methods is the fractional Newton-Raphson method but unlike the latter, the orders of fractional derivatives proposed here are functions. In the first method, a function is used to guarantee an order of convergence (at least) quadratic, and in the others, a function is used to avoid the discontinuity that is generated when the fractional derivative of the constants is used, and with this, it is possible that the methods have at most an order of convergence (at least) linear.

Equations

$$\left\{\xi \in \Omega :f(\xi )=0\right\},$$

$$f(x)=\sin (x)-\frac{1}{x}.$$

$$x_{i+1}=\Phi (x_i),$$

$$\underset{i\to \infty }{\lim }{x}_i\to \xi .$$

$$x_0\in B(\xi ;\delta ):=\left\{y\in {\mathbb{R}}^n:\Vert y-\xi \Vert <\delta \right\},$$

$$\underset{i\to \infty }{\lim }\Vert x_i-\xi \Vert \to 0$$

$$\xi =\underset{i\to \infty }{\lim }{x}_{i+1}=\underset{i\to \infty }{\lim }\Phi (x_i)=\Phi \left(\underset{i\to \infty }{\lim }{x}_i\right)=\Phi (\xi ),$$

$$\Vert \Phi (x)-\Phi (y)\Vert \le \beta \Vert x-y\Vert ,$$

$$\Vert x_{k+1}-\xi \Vert \le \frac{\beta }{1-\beta }\Vert x_{k+1}-x_k\Vert ,$$

$$\Vert x_{k+1}-\xi \Vert \le C{\Vert x_k-\xi \Vert }^p,$$

$$\displaystyle\left\{\begin{array}[]{cc}\left\lvert\Phi^{(k)}(\xi)\right\rvert=0,\ \forall k\leq p-1,&\mbox{if }p\geq 2,\\ \left\lvert\Phi^{(1)}(\xi)\right\rvert<1,&\mbox{if }p=1,\end{array}\right.$$

$$\Phi (x_i)=\Phi (\xi )+\sum _{s=1}^p\frac{{\Phi }^{(s)}(\xi )}{s!}(x_i-\xi {)}^s+o\left((x_i-\xi {)}^p\right),$$

$$\left|\Phi (x_i)-\Phi (\xi )\right|\le \sum _{s=1}^p\left|\frac{{\Phi }^{(s)}(\xi )}{s!}\right|{\left|x_i-\xi \right|}^s+o\left({\left|x_i-\xi \right|}^p\right),$$

$$\frac{\left|x_{i+1}-\xi \right|}{{\left|x_i-\xi \right|}^p}=\frac{\left|\Phi (x_i)-\Phi (\xi )\right|}{{\left|x_i-\xi \right|}^p}\le \left|\frac{{\Phi }^{(p)}(\xi )}{p!}\right|+\frac{o\left({\left|x_i-\xi \right|}^p\right)}{{\left|x_i-\xi \right|}^p}\underset{i\to \infty }{⟶}\left|\frac{{\Phi }^{(p)}(\xi )}{p!}\right|,$$

$$\left|x_{i+1}-\xi \right|\le \left|\frac{{\Phi }^{(p)}(\xi )}{p!}\right|{\left|x_i-\xi \right|}^p,$$

$$\Phi (x)=x-A(x)f(x),$$

$$A(x):=\left(a_{jk}(x)\right)=\left(\begin{array}{cccc}{a}_{11}(x)& a_{12}(x)& \cdots & a_{1n}(x)\\ a_{12}(x)& a_{22}(x)& \cdots & a_{2n}(x)\\ ⋮& ⋮& \ddots & ⋮\\ a_{n1}(x)& a_{n2}(x)& ⋮& a_{nn}(x)\end{array}\right),$$

$${\Phi }^{(1)}(x):=\left({\partial }_j{\Phi }_k(x)\right),$$

$${\partial }_j{\Phi }_k(x):=\frac{\partial }{\partial x_j}{\Phi }_k(x),$$

$$\displaystyle\begin{array}[]{cccc}\displaystyle\lim_{x\to\xi}\left\lVert\Phi^{(1)}(x)\right\rVert=0&\Rightarrow&\displaystyle\lim_{x\to\xi}\partial_{j}\Phi_{k}(x)=0,&\forall j,k\leq n,\end{array}$$

$${\Phi }_k(x)=x_k-\sum _{j=1}^n{a}_{kj}(x)f_j(x),$$

$${\partial }_l{\Phi }_k(x)={\delta }_{lk}-\sum _{j=1}^n\left[a_{kj}(x){\partial }_l{f}_j(x)+\left({\partial }_l{a}_{kj}(x)\right)f_k(x)\right],$$

$$\displaystyle\delta_{lk}=\left\{\begin{array}[]{cc}1,&\mbox{si }l=k,\\ 0,&\mbox{si }l\neq k.\end{array}\right.$$

$${\partial }_l{\Phi }_k(\xi )={\delta }_{lk}-\sum _{j=1}^n{a}_{kj}(\xi ){\partial }_l{f}_j(\xi )=0$$

$$A(\xi )f^{(1)}(\xi )=I_n$$

$$\underset{x\to \xi }{\lim }A(x)={\left(f^{(1)}(\xi )\right)}^{-1},$$

$$x_{i+1}:=\Phi (x_i)=x_i-{\left(f^{(1)}(x_i)\right)}^{-1}f(x_i),$$

$$f(x)=(x-\xi {)}^{m}g(x),$$

$$f(x)=$$

$$f^{(1)}(x)=$$

$$\displaystyle\begin{array}[]{c}\Phi^{(1)}(x)=1-\dfrac{mg^{2}(x)+(x-\xi)^{2}\left[\left(g^{(1)}(x)\right)^{2}-g(x)g^{(2)}(x)\right]}{\left(mg(x)+(x-\xi)g^{(1)}(x)\right)^{2}},\end{array}$$

Full text

Fractional Newton-Raphson Method and Some Variants for the Solution of Non-linear Systems

A. Torres-Hernandez

Department of Physics - UNAM

anthony.torres@ciencias.unam.mx

  

F. Brambila-Paz

Department of Mathematics - UNAM

fernandobrambila@gmail.com

1.   Introduction

When starting to study the fractional calculus, the first difficulty is that, when wanting to solve a problem related to physical units, such as determining the velocity of a particle, the fractional derivative seems to make no sense, this is due to the presence of physical units such as meter and second raised to non-integer exponents, opposite to what happens with operators of integer order. The second difficulty, which is a recurring topic of debate in the study of fractional calculus, is to know what is the order “optimal” $\alpha$ that should be used when the goal is to solve a problem related to fractional operators. To face these difficulties, in the first case, it is common to dimensionless any equation in which non-integer operators are involved, while for the second case different orders $\alpha$ are used in fractional operators to solve some problem, and then it is chosen the order $\alpha$ that provides the “best solution” based about an established criteria.

Based on the two previous difficulties, arises the idea of looking for applications with dimensionless nature such that the need to use multiple orders $\alpha$ can be exploited in some way. The aforementioned led to the study of Newton-Raphson method and a particular problem related to the search for roots in the complex space for polynomials: if it is needed to find a complex root of a polynomial using Newton-Raphson method, it is necessary to provide a complex initial condition $x_{0}$ and, if the right conditions are selected, this leads to a complex solution, but there is also the possibility that this leads to a real solution. If the root obtained is real, it is necessary to change the initial condition and expect that this leads to a complex solution, otherwise, it is necessary to change the value of the initial condition again.

The process described above, it is very similar to what happens when using different values $\alpha$ in fractional operators until we find a solution that meets some desired condition. Seeing Newton-Raphson method from the perspective of fractional calculus, one can consider that an order $\alpha$ remains fixed , in this case $\alpha=1$, and the initial conditions $x_{0}$ are varied until obtaining one solution that satisfies an established criteria. Then reversing the behavior of $\alpha$ and $x_{0}$, i.e., leave the initial condition fixed and varying the order $\alpha$, the fractional Newton-Raphson method is obtained, which is nothing other thing than Newton-Raphson method using any definition of fractional derivative that fits the function with which one is working. This change, although in essence simple, allows to find roots in the complex space using real initial conditions, because fractional operators generally do not carry polynomials to polynomials.

1.1.   Fixed Point Method

A classic problem of common interest in both Physics and Mathematics is to find the zeros of a function $f:\Omega\subset\mathbb{R}^{n}\to\mathbb{R}^{n}$, i.e.,

[TABLE]

this problem often arises as a consequence of wanting to solve other problems, for instance, if we want to determine the eigenvalues of a matrix or want to build a box with a given volume but with minimal surface; in the first example, we need to find the zeros (or roots) of the characteristic polynomial of the matrix, while in the second one we need to find the zeros of the gradient of a function that relates the surface of the box with its volume.

Although finding the zeros of a function may seem like a simple problem, in general, it involves solving nonlinear equations, which in many cases does not have an analytical solution, an example of this is present when we are trying to determine the zeros of the following function

[TABLE]

Because in many cases there is no analytical solution, numerical methods are needed to try to determine the solutions to these problems; it should be noted that when using numerical methods, the word “determine” should be interpreted as approach a solution with a degree of precision desired. The numerical methods mentioned above are usually of the iterative type and work as follows: suppose we have a function $f:\Omega\subset\mathbb{R}^{n}\to\mathbb{R}^{n}$ and we search a value $\xi\in\mathbb{R}^{n}$ such that $f(\xi)=0$, then we can start by giving an initial value $x_{0}\in\mathbb{R}^{n}$ and then calculate a value $x_{i}$ close to the searched value $\xi$ using an iteration function $\Phi:\mathbb{R}^{n}\to\mathbb{R}^{n}$ as follows [3]:

[TABLE]

this generates a sequence $\left\{x_{i}\right\}_{i=0}^{\infty}$, which under certain conditions satisfies that

[TABLE]

To understand the convergence of the iteration function $\Phi$ it is necessary to have the following definition [4]:

Definition 1.1**.**

Let $\Phi:\mathbb{R}^{n}\to\mathbb{R}^{n}$ be an iteration function. The method given in (1) to determine $\xi\in\mathbb{R}^{n}$, it is called (locally) convergent, if exists $\delta>0$ such that for all initial value

[TABLE]

it holds that

[TABLE]

where $\left\lVert\ \cdot\ \right\rVert:\mathbb{R}^{n}\to\mathbb{R}$ denotes any vector norm.

When it is assumed that the iteration function $\Phi$ is continuous at $\xi$ and that the sequence $\left\{x_{i}\right\}_{i=0}^{\infty}$ converges to $\xi$ under the condition given in (2), it is true that

[TABLE]

the previous result is the reason why the method given in (1) is called Fixed Point Method.

1.1.1 Convergence and Order of Convergence

The (local) convergence of the iteration function $\Phi$ established in (2), it is useful for demonstrating certain intrinsic properties of the fixed point method. Before continuing it is necessary to have the following definition [3]:

Definition 1.2**.**

Let $\Phi:\Omega\subset\mathbb{R}^{n}\to\mathbb{R}^{n}$. The function $\Phi$ is a contraction on a set $\Omega_{0}\subset\Omega$, if exists a non-negative constant $\beta<1$ such that

[TABLE]

where $\beta$ is called a contraction constant.

The previous definition guarantee that if the iteration function $\Phi$ is a contraction on a set $\Omega_{0}$, then it is Lipschitz continuous and has at least one fixed point. The existence of a fixed point is guaranteed by the following theorem [5]:

** Theorem 1.3****.**

Contraction Mapping Theorem: Let $\Phi:\Omega\subset\mathbb{R}^{n}\to\mathbb{R}^{n}$. Assuming that $\Phi$ is a contraction on a closed set $\Omega_{0}\subset\Omega$, and that $\Phi(x)\in\Omega_{0}\ \forall x\in\Omega_{0}$. Then $\Phi$ has a unique fixed point $\xi\in\Omega_{0}$ and for any initial value $x_{0}\in\Omega_{0}$, the sequence $\left\{x_{i}\right\}_{i=0}^{\infty}$ generated by (1) converges to $\xi$. Moreover

[TABLE]

where $\beta$ is the contraction constant given in (4).

When the fixed point method given by (1) is used, in addition to convergence, exists a special interest in the order of convergence, which is defined as follows [4]:

Definition 1.4**.**

Let $\Phi:\Omega\subset\mathbb{R}^{n}\to\mathbb{R}^{n}$ be an iteration function with a fixed point $\xi\in\Omega$. Then the method (1) is called (locally) convergent of (at least) order $\boldsymbol{p}$ ($p\geq 1$), if exists $\delta>0$ and exists a non-negative constant $C$ (with $C<1$ if $p=1$) such that for any initial value $x_{0}\in B(\xi;\delta)$ it is true that

[TABLE]

where $C$ is called convergence factor.

The order of convergence is usually related to the speed at which the sequence generated by (1) converges. For the particular cases $p=1$ or $p=2$ it is said that the method has (at least) linear or quadratic convergence, respectively. The following theorem for the one-dimensional case [4], allows characterizing the order of convergence of an iteration function $\Phi$ with its derivatives

** Theorem 1.5****.**

Let $\Phi:\Omega\subset\mathbb{R}\to\mathbb{R}$ be an iteration function with a fixed point $\xi\in\Omega$. Assuming that $\Phi$ is $p$-times differentiable in $\xi$ for some $p\in\mathbb{N}$, and furthermore

[TABLE]

then $\Phi$ is (locally) convergent of (at least) order $p$.

Proof.

Assuming that $\Phi:\Omega\subset\mathbb{R}\to\mathbb{R}$ is an iteration function $p$-times differentiable with a fixed point $\xi\in\Omega$, then we can expand in Taylor series the function $\Phi(x_{i})$ around $\xi$ and order $p$

[TABLE]

then we obtain

[TABLE]

assuming that the sequence $\left\{x_{i}\right\}_{i=0}^{\infty}$ generated by (1) converges to $\xi$ and also that $\left\lvert\Phi^{(s)}(\xi)\right\rvert=0\ \forall s<p$, the previous expression implies that

[TABLE]

as consequence, exists a value $k>0$ such that

[TABLE]

∎

A version of the previous theorem for the case $n$-dimensional may be found in the reference [3].

1.2.   Newton-Raphson Method

The previous theorem in its n-dimensional form is usually very useful to generate a fixed point method with an order of convergence desired, an order of convergence that is usually appreciated in iterative methods is the (at least) quadratic order. If we have a function $f:\Omega\subset\mathbb{R}^{n}\to\mathbb{R}^{n}$ and we search a value $\xi\in\Omega$ such that $f(\xi)=0$, we may build an iteration function $\Phi$ in general form as [6]:

[TABLE]

with $A(x)$ a matrix like

[TABLE]

where $a_{jk}:\mathbb{R}^{n}\to\mathbb{R}$ ($1\leq j,k\leq n$). Notice that the matrix $A(x)$ is determined according to the order of convergence desired. Before continuing, it is necessary to mention that the following conditions are needed:

Suppose we can generalize the Theorem 1.5 to the case $n$-dimensional, although for this it is necessary to guarantee that the iteration function $\Phi$ given by (10) near the value $\xi$ can be expressed in terms of its Taylor series in several variables. 2. 2.

It is necessary to guarantee that the norm of the equivalent of the first derivative in several variables of the iteration function $\Phi$ tends to zero near the value $\xi$.

Then, we will assume that the first condition is satisfied; for the second condition we have that the equivalent to the first derivative in several variables is the Jacobian matrix of the function $\Phi$, which is defined as follows [5]:

[TABLE]

where

[TABLE]

with $\Phi_{k}:\mathbb{R}^{n}\to\mathbb{R}$, the competent $k$-th of the iteration function $\Phi$. Now considering that

[TABLE]

Assume that we have a function $f(x):\Omega\subset\mathbb{R}^{n}\to\mathbb{R}^{n}$ with a zero $\xi\in\Omega$, such that all of its first partial derivatives are defined in $\xi$. Then taking the iteration function $\Phi$ given by (10), the $k$-th component of the iteration function may be written as

[TABLE]

then

[TABLE]

where $\delta_{lk}$ is the Kronecker delta, which is defined as

[TABLE]

Assuming that (14) is fulfilled

[TABLE]

the previous expression may be written in matrix form as

[TABLE]

where $f^{(1)}$ and $I_{n}$ are the Jacobian matrix of the function $f$ and the identity matrix of $n\times n$, respectively. Denoting by $\mbox{det}\left(A\right)$ the determinant of the matrix $A$, then any matrix $A(x)$ that fulfill the following condition

[TABLE]

guarantees that exists $\delta>0$ such that the iteration function $\Phi$ given by (10), converges (locally) with an order of convergence (at least) quadratic in $B(\xi;\delta)$. The following fixed point method can be obtained from the previous result

[TABLE]

which is known as Newton-Raphson method ($n$-dimensional), also known as Newton’s method [7].

Although the condition given in (16) could seem that Newton-Raphson method always has an order of convergence (at least) quadratic, unfortunately, this is not true; the order of convergence is now conditioned by the way in which the function f is constituted, the mentioned above may be appreciated in the following proposition

Proposition 1.6**.**

Let $f:\Omega\subset\mathbb{R}\to\mathbb{R}$ be a function that is at least twice differentiable in $\xi\in\Omega$. So if $\xi$ is a zero of $f$ with algebraic multiplicity $m$ ($m\geq 2$), i.e.,

[TABLE]

the Newton-Raphson method (one-dimensional) has an order of convergence (at least) linear.

Proof.

Suppose we have $f:\Omega\subset\mathbb{R}\to\mathbb{R}$ a function with a zero $\xi\in\Omega$ of algebraic multiplicity $m\geq 2$, and that $f$ is at least twice differentiable in $\xi$, then

[TABLE]

as consequence the derivative of the iteration function $\Phi$ of Newton-Raphson method may be expressed as

[TABLE]

therefore

[TABLE]

and by the Theorem 1.5, then Newton-Raphson method under the hypothesis of the proposition converges (locally) with an order of convergence (at least) linear.

∎

Notice that the order of convergence (at least) quadratic, does not necessarily imply a higher speed of convergence when using Newton-Raphson method; to verify this, just is necessary to take a polynomial of the form

[TABLE]

now assuming that $\xi\in\mathbb{R}$ is a zero of the polynomial and that the sequence $\left\{x_{i}\right\}_{i=0}^{\infty}$ generated by (17) has an initial condition $x_{0}\in\mathbb{R}$ away from $\xi$, then

[TABLE]

so a small change is obtained between the values $x_{i}$ and $x_{i+1}$; therefore, the sequence generated by Newton method will have a slow convergence to the value of $\xi$ [3].

2.   Basic Definitions of the Fractional Derivative

2.1.   Introduction to the Definition of Riemann-Liouville

One of the key pieces in the study of fractional calculus is the iterated integral, which is defined as follows [8]:

Definition 2.1**.**

Let $L_{loc}^{1}(a,b)$, the space of locally integrable functions in the interval $(a,b)$. If $f$ is a function such that $f\in L_{loc}^{1}(a,\infty)$, then the $n$-th iterated integral of the function $f$ is given by [8]:

[TABLE]

where

[TABLE]

Considerate that $(n-1)!=\Gamma\left(n\right)$ , a generalization of (20) may be obtained for an arbitrary order $\alpha>0$

[TABLE]

similarly, if $f\in L_{loc}^{1}(-\infty,b)$, we may define

[TABLE]

the equations (21) and (22) correspond to the definitions of right and left fractional integral of Riemann-Liouville, respectively.

Below, there are some valid properties for right fractional integrals. However, the same properties for left fractional integrals have some changes. Fractional integrals satisfy the semigroup property, which is given in the following proposition [8]:

Proposition 2.2**.**

Let $f$ be a function. If $f\in L_{loc}^{1}(a,\infty)$, then the fractional integrals of $f$ satisfy that

[TABLE]

From the previous result, it is obtained that in particular

[TABLE]

since the $d/dx$ operator is the inverse operator to the left of the operator ${}_{a}{I}_{x}$, any integral $\alpha$-th of a function $f\in L_{loc}^{1}(a,\infty)$ may be written as

[TABLE]

It should be mentioned that the above results are also valid for $\alpha$ in complexes, with $Re\left(\alpha\right)>0$. Under appropriate conditions in (24) its follows

[TABLE]

then, from (24) and (25) the following operator is built, which corresponds to the ** (right) fractional derivative of Riemann-Liouville**

[TABLE]

where $n=\lfloor Re\left(\alpha\right)\rfloor+1$, with $\lfloor Re\left(\alpha\right)\rfloor$ the largest integer less than or equal to $Re\left(\alpha\right)$. We can unify the definitions of integral and fractional derivative of Riemann-Liouville, given by (21) and (26), as follows

[TABLE]

where $n=\lfloor Re\left(\alpha\right)\rfloor+1$. Applying the operator (30) with $a=0$ and $\alpha\in\mathbb{R}\setminus\mathbb{Z}$ to the function $x^{\mu}$, we obtain

[TABLE]

2.2.   Introduction to the Definition of Caputo

Michele Caputo (1969) published a book and introduced a new definition of fractional derivative, he created this definition with the objective of modeling anomalous diffusion phenomena. The definition of Caputo had already been discovered independently by Gerasimov (1948). This fractional derivative is of the utmost importance since it allows us to give a physical interpretation of the initial value problems, moreover to being used to model fractional time. In some texts, it is known as the fractional derivative of Gerasimov-Caputo [9].

Let $f$ be a function, such that $f$ is $n$-times differentiable with $f^{(n)}\in L_{loc}^{1}(a,b)$, then the (right) fractional derivative of Caputo is defined as:

[TABLE]

where $n=\lfloor Re\left(\alpha\right)\rfloor+1$. It should be mentioned that the fractional derivative of Caputo behaves as the inverse operator to the left of fractional integral of Riemann-Liouville , that is,

[TABLE]

On the other hand, the relation between the fractional derivatives of Caputo and Riemann-Liouville is given by the following expression [9]

[TABLE]

then if $f^{(k)}(a)=0\ \ \forall k<n$, we obtain

[TABLE]

considering the previous particular case, it is possible to unify the definitions of fractional integral of Riemann-Liouville and fractional derivative of Caputo as follows

[TABLE]

3.   Fractional Newton Method

Let $\mathbb{P}_{n}(\mathbb{R})$, the space of polynomials of degree less than or equal to $n$ with real coefficients. The zeros $\xi$ of a function $f\in\mathbb{P}_{n}(\mathbb{R})$ are usually named as roots. The Newton-Raphson method is useful for finding the roots of a function $f$. However, this method is limited because it cannot find roots $\xi\in\mathbb{C}\setminus\mathbb{R}$, if the sequence $\left\{x_{i}\right\}_{i=0}^{\infty}$ generated by (17) has an initial condition $x_{0}\in\mathbb{R}$. To solve this problem and develop a method that has the ability to find roots, both real and complex, of a polynomial if the initial condition $x_{0}$ is real, we propose a new method called fractional Newton-Raphson method, which consists of Newton-Raphson method with the implementation of the fractional derivative. Before continuing, it is necessary to define the fractional Jacobian matrix of a function $f:\Omega\subset\mathbb{R}^{n}\to\mathbb{R}^{n}$ as follows

[TABLE]

where

[TABLE]

with $f_{k}:\mathbb{R}^{n}\to\mathbb{R}$. The operator $\partial^{\alpha}/\partial x_{j}^{\alpha}$ denotes any fractional derivative, applied only to the variable $x_{j}$, that satisfies the following condition of continuity respect to the order of the derivative

[TABLE]

then, the matrix (37) satisfies that

[TABLE]

where $f^{(1)}(x)$ denotes the Jacobian matrix of the function $f$.

Taking into account that a polynomial of degree $n$ it is composed of $n+1$ monomials of the form $x^{m}$, with $m\geq 0$, we can take the equation (31) with (17), to define the following iteration function that results in the Fractional Newton-Raphson Method [1, 2]

[TABLE]

To try to guarantee that the sequence $\left\{x_{i}\right\}_{i=0}^{\infty}$ generated by (40) has a order of convergence (at least) quadratic, the condition (16) is combined with (38) to define the following function

[TABLE]

then, in (40) we make $\alpha\to\alpha_{f}(x)$ obtaining the following iteration function that results in the Fractional Newton Method

[TABLE]

the difference between the methods (40) and (45), is that the just for the second can exists $\delta>0$ such that if the sequence $\left\{x_{i}\right\}_{i=0}^{\infty}$ generated by (45) converges to a root $\xi$ of $f$, exists $k>0$ such that $\forall i\geq k$, the sequence has an order of convergence (at least) quadratic in $B(\xi;\delta)$.

The value of $\alpha$ in (43) is assigned with the following reasoning: when we use the definitions of fractional derivatives given by (30) and (36) in a function $f$, it is necessary that the function be $n$-times integrable and $n$-times differentiable, where $n=\lfloor Re\left(\alpha\right)\rfloor+1$, therefore $\left\lvert\alpha\right\rvert<n$ and, on the other hand, for using Newton method it is just necessary that the function be one-time differentiable, as a consequence of $\eqref{eq:c2.21}$ it is obtained that

[TABLE]

Without loss of generality, to understand why the sequence $\left\{x_{i}\right\}_{i=0}^{\infty}$ generated by the fractional Newton method when we use a function $f\in P_{n}(\mathbb{R})$, has the ability to enter the complex space starting from an initial condition $x_{0}\in\mathbb{R}$, it is only necessary to observe the fractional derivative of Riemann-Liouville of order $\alpha=1/2$ of the monomial $x^{m}$

[TABLE]

whose result is a function with a rational exponent, contrary to what would happen when using the conventional derivative. When the iteration function given by (45) is used, we must taken an initial condition $x_{0}\neq 0$, as a consequence of the fact that the fractional derivative of order $\alpha>0$ of a constant is usually proportional to the function $x^{-\alpha}$.

When we take $\alpha\neq 1$, the sequence $\left\{x_{i}\right\}_{i=0}^{\infty}$ generated by the iteration function (45), presents among its behaviors, the following particular cases depending on the initial condition $x_{0}$:

If we take an initial condition $x_{0}>0$, the sequence $\left\{x_{i}\right\}_{i=0}^{\infty}$ may be divided into three parts, this occurs because it may exists a value $M\in\mathbb{N}$ for which $\left\{x_{i}\right\}_{i=0}^{M-1}\subset\mathbb{R}^{+}$ with $\left\{x_{M}\right\}\subset\mathbb{R}^{-}$, in consequence $\left\{x_{i}\right\}_{i\geq M+1}\subset\mathbb{C}$. 2. 2.

On the other hand, if we take an initial $x_{0}<0$ condition, the sequence $\left\{x_{i}\right\}_{i=0}^{\infty}$ may be divided into two parts, $\left\{x_{0}\right\}\subset\mathbb{R}^{-}$ and $\left\{x_{i}\right\}_{i\geq 1}\subset\mathbb{C}$.

Unlike classical Newton-Raphson method; which uses tangent lines to generate a sequence $\left\{x_{i}\right\}_{i=0}^{\infty}$, the fractional Newton method uses lines more similar to secants (see Figure 1). A consequence of the fact that the lines are not tangent when using (45), is that different trajectories can be obtained for the same initial condition $x_{0}$ just by changing the order $\alpha$ of the derivative (see Figure 2).

A typical inconvenience that arises in problems related to fractional calculus, it is the fact that it is not always known what is the appropriate order $\alpha$ to solve these problems. As a consequence, different values of $\alpha$ are generally tested and we choose the value that allows to find the best solution considering an criterion of precision established. Based on the aforementioned, it is necessary to follow the instructions below when using the fractional Newton method to find the zeros $\xi$ of a function $f$:

Without considering the integers $-1,\ 0$ and $1$, a partition of the interval $[-2,2]$ is created as follows

[TABLE]

and using the partition, the following sequence $\left\{\alpha_{m}\right\}_{m=1}^{s}$ is created. 2. 2.

We choose a non-negative tolerance $TOL<1$, a limit of iterations $L_{IT}>1$ for all $\alpha_{m}$, an initial condition $x_{0}\neq 0$ and a value $M>L_{IT}$. 3. 3.

We choose a value $\delta>0$ to use $\alpha_{f}(x)$ given by (43), such that $TOL<\delta<1$. In addition, it is taken a fractional derivative that satisfies the condition of continuity (38), and it is unified with the fractional integral in the same way as in the equations (30) and (36). 4. 4.

The iteration function (45) is used with all the values of the partition $\left\{\alpha_{m}\right\}_{m=1}^{s}$, and for each value $\alpha_{m}$ is generated a sequence $\left\{{}^{m}x_{i}\right\}_{i=0}^{R_{m}}$, where

[TABLE]

then is generated a sequence $\left\{x_{R_{m_{k}}}\right\}_{k=1}^{r}$ , with $r\leq s$, such that

[TABLE] 5. 5.

We choose a value $\varepsilon>0$ and we take the values $x_{R_{m_{1}}}$ and $x_{R_{m_{2}}}$, then is defined $X_{1}=x_{R_{m_{1}}}$. If the following condition is fulfilled

[TABLE]

is taken $X_{1}=x_{R_{m_{2}}}$. On the other hand if

[TABLE]

is defined $X_{2}=x_{R_{m_{2}}}$. Without loss of generality, it may be assumed that the second condition is fulfilled, then is taken $X_{3}=X_{R_{m_{3}}}$ and are checked the conditions (49) and (50) with the values $X_{1}$ and $X_{2}$. The process described before is repeated for all values $X_{k}=x_{R_{m_{k}}}$, with $k\geq 4$, and that generates a sequence $\left\{X_{m}\right\}_{m=1}^{t}$, with $t\leq r$, such that

[TABLE]

By following the steps described before to implement the fractional Newton method, a subset of the solution set of zeros, both real and complex, may be obtained from the function $f$. We will proceed to give an example where is found the solution set of zeros of a function $f\in\mathbb{P}_{n}(\mathbb{R})$.

Example 3.1**.**

Let the function:

[TABLE]

then the following values are chosen to use the iteration function given by (45)

[TABLE]

and using the fractional derivative given by (30), we obtain the results of the Table 1

Although the fractional Newton method was originally defined for polynomials, the method can be extended to a broader class of functions, as shown in the following examples:

Example 3.2**.**

Let the function:

[TABLE]

then the following values are chosen to use the iteration function given by (45)

[TABLE]

and using the fractional derivative given by (30), we obtain the results of Table 2

In the previous example, a subset of the solution set of zeros of the function (56) was obtained, because this function has an infinite amount of zeros. Using the fractional Newton method does not guarantee that all zeros of a function $f$ can be found, leaving an initial condition $x_{0}$ fixed and varying the orders $\alpha_{m}$ of the derivative. As in the classical Newton’s method, finding most of the zeros of the function will depend on giving a proper initial condition $x_{0}$. If we want to find a larger subset of zeros of the function (56), there are some strategies that are usually useful, for example:

To change the initial condition $x_{0}$. 2. 2.

To use a larger amount of $\alpha_{m}$ values. 3. 3.

To increase the value of $L_{IT}$.

In general, the last strategy is usually the most useful, but this causes the method (45) to become more costly, because it requires a longer runtime for all values $\alpha_{m}$.

Example 3.3**.**

Let the function:

[TABLE]

then the following values are chosen to use the iteration function given by (45)

[TABLE]

and using the fractional derivative given by (30), we obtain the results of Table 3

Example 3.4**.**

Let the function:

[TABLE]

then the following values are chosen to use the iteration function given by (45)

[TABLE]

and using the fractional derivative given by (30), we obtain the results of Table 4

3.1.    Fractional Quasi-Newton Method

Although the fractional Newton method is useful for finding multiple zeros of a function $f$, it has the disadvantage that in many cases calculating the fractional derivative of a function is not a simple task. To try to minimize this problem, we use that commonly for many definitions of the fractional derivative, the arbitrary order derivative of a constant is not always zero, that is,

[TABLE]

Then, we may define the function

[TABLE]

it should be noted that the previous function is almost a linear approximation of the function $f$ in the initial condition $x_{0}$. Then, for any fractional derivative that satisfies the condition (65), and using (37) the Fractional Quasi-Newton Method may be defined as

[TABLE]

with $Q_{g,\beta}(x_{i})$ given by the following matrix

[TABLE]

where $(x_{i})_{j}$ is the component $j$-th of the value $x_{i}$, and $\beta(\alpha,(x_{i})_{j})$ is defined as follows

[TABLE]

Since the iteration function (68) does not satisfy the condition (16), then any sequence $\left\{x_{i}\right\}_{i=0}^{\infty}$ generated by this iteration function has at most one order of convergence (at least) linear. As a consequence, the speed of convergence is slower compared to what would be obtained when using (45), and then it is necessary to use a larger value of $L_{IT}$. It should be mentioned that the value $\alpha=1$ in (72), is not taken to try to guarantee an order of convergence as in (43), but to avoid the discontinuity that is generated when using the fractional derivative of constants in the value $x=0$. An example is given using the fractional quasi-Newton method, where is found a subset of the solution set of zeros of the function $f$.

Example 3.5**.**

Let the function:

[TABLE]

then the following values are chosen to use the iteration function given by (68)

[TABLE]

and using the fractional derivative given by (30), we obtain the results of Table 5

It is worth mentioning that the function (74) has infinite zeros, so it will be used with the following method.

3.2.    Fractional Pseudo-Newton Method

A conventional way in which Newton’s method is usually deduced is based on a linear approximation of a function $f$ in a value $x_{i}$. To perform the aforementioned, we may consider the Taylor polynomial of the function $f$ around a value $x_{i}$, and disregarding higher order terms we have that

[TABLE]

then, assuming that $\xi$ is a zero of $f$, from the previous expression we have that

[TABLE]

as consequence, a sequence $\left\{x_{i}\right\}_{i=0}^{\infty}$ that approximates the value $\xi$ may be generated using the iteration function

[TABLE]

However, the equation (76) is not the only way to generate a linear approximation to the function $f$ in the point $x_{i}$, in general it may be taken as

[TABLE]

where $m$ is any constant value of a slope, that allows the approximation (77) to the function $f$ to be valid. The previous equation allows to obtain the following iteration function

[TABLE]

which originates the Parallel Chord Method [7]. The iteration function (78) can be generalized to larger dimensions as follows

[TABLE]

where $I_{n}$, is the identity matrix of $n\times n$. The previous equation implies that it is enough to apply the method (78) component to component for the case in several variables.

Using as a basis the idea of (79), and considering any fractional derivative that satisfies the condition (65), we can define the Fractional Pseudo-Newton Method as follows

[TABLE]

with $P_{\epsilon,\beta}(x_{i})$ given by the following matrix evaluated at the value $x_{i}$

[TABLE]

where $(x_{i})_{j}$ is the component $j$-th of the value $x_{i}$, with $\delta_{kj}$ the Kronecker delta, $\epsilon$ a positive constant $\ll 1$, and $\beta(\alpha,(x_{i})_{j})$ is given by (72).

Two examples are given using the fractional pseudo-Newton method, where is found a subset of the solution set of zeros of some functions.

Example 3.6**.**

Let the function:

[TABLE]

then the following values are chosen to use the iteration function given by (81)

[TABLE]

and using the fractional derivative given by (30), we obtain the results of Table 6

Example 3.7**.**

Let the function:

[TABLE]

then the following values are chosen to use the iteration function given by (81)

[TABLE]

and using the fractional derivative given by (30), we obtain the results of Table 7

4.   Conclusions

The fractional Newton-Raphson method and its variants are useful for finding multiple solutions of nonlinear systems, in the complex space using real initial conditions. However, it should be clarified that they present some advantages and disadvantages between each of them, for example, although the fractional Newton method generally has an order of convergence (at least) quadratic, this method has the disadvantage that it is not an easy task to find the fractional Jacobian matrix for many functions, and also the need to reverse this matrix must be added for each new iteration. But it has an advantage over the other methods, and this is because it can be used with few iterations, which allows occupying a greater number of $\alpha_{m}$ values belonging to the partition of the interval $(-2,2)$.

The quasi-Newton fractional method has the advantage that the fractional Jacobian matrix with which it works, compared to the fractional Newton method, is easy to obtain. But a disadvantage is that the method may have at most an order of convergence (at least) linear, so the speed of convergence is lower and it is necessary to use a greater number of iterations to ensure success in the search for solutions. As a consequence, the method is more costly because it requires a longer runtime to use all values $\alpha_{m}$. An additional advantage of the method is that if the initial condition is close enough to a solution, its behavior is very similar to fractional Newton method and may converge with a relatively small number of iterations, but it still has the disadvantage that we need to reverse a matrix in each iteration.

Finally, the fractional pseudo-Newton method seems to solve the problem of the need to invert a matrix in each iteration that are present in the previous methods. However, this method, as the fractional quasi-Newton method, may have at most an order of convergence (at least) linear, and it has a lower speed of convergence compared to the latter. As a consequence, it is necessary to use a larger number of values $\alpha_{m}$ and a greater number of iterations, then we require a longer runtime to use all values to find solutions, so it can be considered as a slow and costly method. But its advantage over the other methods is that it is not necessary to calculate and invert any fractional Jacobian matrix.

The methods may solve some nonlinear systems and are really efficient to find multiple solutions, both real and complex, using real initial conditions. It should be mentioned that these methods are extremely recommended in systems that have infinite solutions or a large number of them.