Sidi's generalized secant method

Sidi's generalized secant method is a root-finding algorithm, that is, a numerical method for solving equations of the form [math]\displaystyle{ f(x)=0 }[/math] . The method was published by Avram Sidi.

The method is a generalization of the secant method. Like the secant method, it is an iterative method which requires one evaluation of [math]\displaystyle{ f }[/math] in each iteration and no derivatives of [math]\displaystyle{ f }[/math]. The method can converge much faster though, with an order which approaches 2 provided that [math]\displaystyle{ f }[/math] satisfies the regularity conditions described below.

Algorithm

We call [math]\displaystyle{ \alpha }[/math] the root of [math]\displaystyle{ f }[/math], that is, [math]\displaystyle{ f(\alpha)=0 }[/math]. Sidi's method is an iterative method which generates a sequence [math]\displaystyle{ \{ x_i \} }[/math] of approximations of [math]\displaystyle{ \alpha }[/math]. Starting with k + 1 initial approximations [math]\displaystyle{ x_1 , \dots , x_{k+1} }[/math], the approximation [math]\displaystyle{ x_{k+2} }[/math] is calculated in the first iteration, the approximation [math]\displaystyle{ x_{k+3} }[/math] is calculated in the second iteration, etc. Each iteration takes as input the last k + 1 approximations and the value of [math]\displaystyle{ f }[/math] at those approximations. Hence the nth iteration takes as input the approximations [math]\displaystyle{ x_n , \dots , x_{n+k} }[/math] and the values [math]\displaystyle{ f(x_n) , \dots , f(x_{n+k}) }[/math].

The number k must be 1 or larger: k = 1, 2, 3, .... It remains fixed during the execution of the algorithm. In order to obtain the starting approximations [math]\displaystyle{ x_1 , \dots , x_{k+1} }[/math] one could carry out a few initializing iterations with a lower value of k.

The approximation [math]\displaystyle{ x_{n+k+1} }[/math] is calculated as follows in the nth iteration. A polynomial of interpolation [math]\displaystyle{ p_{n,k} (x) }[/math] of degree k is fitted to the k + 1 points [math]\displaystyle{ (x_n, f(x_n)), \dots (x_{n+k}, f(x_{n+k})) }[/math]. With this polynomial, the next approximation [math]\displaystyle{ x_{n+k+1} }[/math] of [math]\displaystyle{ \alpha }[/math] is calculated as

[math]\displaystyle{ x_{n+k+1} = x_{n+k} - \frac{f(x_{n+k})}{p_{n,k}'(x_{n+k})} }[/math]

(1)

with [math]\displaystyle{ p_{n,k}'(x_{n+k}) }[/math] the derivative of [math]\displaystyle{ p_{n,k} }[/math] at [math]\displaystyle{ x_{n+k} }[/math]. Having calculated [math]\displaystyle{ x_{n+k+1} }[/math] one calculates [math]\displaystyle{ f(x_{n+k+1}) }[/math] and the algorithm can continue with the (n + 1)th iteration. Clearly, this method requires the function [math]\displaystyle{ f }[/math] to be evaluated only once per iteration; it requires no derivatives of [math]\displaystyle{ f }[/math].

The iterative cycle is stopped if an appropriate stopping criterion is met. Typically the criterion is that the last calculated approximation is close enough to the sought-after root [math]\displaystyle{ \alpha }[/math].

To execute the algorithm effectively, Sidi's method calculates the interpolating polynomial [math]\displaystyle{ p_{n,k} (x) }[/math] in its Newton form.

Convergence

Sidi showed that if the function [math]\displaystyle{ f }[/math] is (k + 1)-times continuously differentiable in an open interval [math]\displaystyle{ I }[/math] containing [math]\displaystyle{ \alpha }[/math] (that is, [math]\displaystyle{ f \in C^{k+1} (I) }[/math]), [math]\displaystyle{ \alpha }[/math] is a simple root of [math]\displaystyle{ f }[/math] (that is, [math]\displaystyle{ f'(\alpha) \neq 0 }[/math]) and the initial approximations [math]\displaystyle{ x_1 , \dots , x_{k+1} }[/math] are chosen close enough to [math]\displaystyle{ \alpha }[/math], then the sequence [math]\displaystyle{ \{ x_i \} }[/math] converges to [math]\displaystyle{ \alpha }[/math], meaning that the following limit holds: [math]\displaystyle{ \lim\limits_{n \to \infty} x_n = \alpha }[/math].

Sidi furthermore showed that

[math]\displaystyle{ \lim_{n\to\infty} \frac{x_{n +1}-\alpha}{\prod^k_{i=0}(x_{n-i}-\alpha)} = L = \frac{(-1)^{k+1}} {(k+1)!}\frac{f^{(k+1)}(\alpha)}{f'(\alpha)}, }[/math]

and that the sequence converges to [math]\displaystyle{ \alpha }[/math] of order [math]\displaystyle{ \psi_k }[/math], i.e.

[math]\displaystyle{ \lim\limits_{n \to \infty} \frac{|x_{n+1}-\alpha|}{|x_n-\alpha|^{\psi_k}} = |L|^{(\psi_k-1)/k} }[/math]

The order of convergence [math]\displaystyle{ \psi_k }[/math] is the only positive root of the polynomial

[math]\displaystyle{ s^{k+1} - s^k - s^{k-1} - \dots - s - 1 }[/math]

We have e.g. [math]\displaystyle{ \psi_1 = (1+\sqrt{5})/2 }[/math] ≈ 1.6180, [math]\displaystyle{ \psi_2 }[/math] ≈ 1.8393 and [math]\displaystyle{ \psi_3 }[/math] ≈ 1.9276. The order approaches 2 from below if k becomes large: [math]\displaystyle{ \lim\limits_{k \to \infty} \psi_k = 2 }[/math]

Related algorithms

Sidi's method reduces to the secant method if we take k = 1. In this case the polynomial [math]\displaystyle{ p_{n,1} (x) }[/math] is the linear approximation of [math]\displaystyle{ f }[/math] around [math]\displaystyle{ \alpha }[/math] which is used in the nth iteration of the secant method.

We can expect that the larger we choose k, the better [math]\displaystyle{ p_{n,k} (x) }[/math] is an approximation of [math]\displaystyle{ f(x) }[/math] around [math]\displaystyle{ x=\alpha }[/math]. Also, the better [math]\displaystyle{ p_{n,k}' (x) }[/math] is an approximation of [math]\displaystyle{ f'(x) }[/math] around [math]\displaystyle{ x=\alpha }[/math]. If we replace [math]\displaystyle{ p_{n,k}' }[/math] with [math]\displaystyle{ f' }[/math] in (1) we obtain that the next approximation in each iteration is calculated as

[math]\displaystyle{ x_{n+k+1} = x_{n+k} - \frac{f(x_{n+k})}{f'(x_{n+k})} }[/math]

(2)

This is the Newton–Raphson method. It starts off with a single approximation [math]\displaystyle{ x_1 }[/math] so we can take k = 0 in (2). It does not require an interpolating polynomial but instead one has to evaluate the derivative [math]\displaystyle{ f' }[/math] in each iteration. Depending on the nature of [math]\displaystyle{ f }[/math] this may not be possible or practical.

Once the interpolating polynomial [math]\displaystyle{ p_{n,k} (x) }[/math] has been calculated, one can also calculate the next approximation [math]\displaystyle{ x_{n+k+1} }[/math] as a solution of [math]\displaystyle{ p_{n,k} (x)=0 }[/math] instead of using (1). For k = 1 these two methods are identical: it is the secant method. For k = 2 this method is known as Muller's method. For k = 3 this approach involves finding the roots of a cubic function, which is unattractively complicated. This problem becomes worse for even larger values of k. An additional complication is that the equation [math]\displaystyle{ p_{n,k} (x)=0 }[/math] will in general have multiple solutions and a prescription has to be given which of these solutions is the next approximation [math]\displaystyle{ x_{n+k+1} }[/math]. Muller does this for the case k = 2 but no such prescriptions appear to exist for k > 2.