Background: I needed some efficient but simple interpolation-methods aside of Newton's iteration because I want to have it in contexts, where the derivative of a function is not always known. So an obvious choice after bisection is regula-falsi and/or secant method. Now after a short course in wikipedia there it is mentioned, that the secant-method has convergence rate of about 1.6 (when it begins to converge), and a parabolic method named after Müller (using a quadratic interpolation-formula instead of the linear secant-formula) is a slight improvement.
Well, I thought that it would be interesting to extend this over the quadratic to cubic and finally dynamically increasing order, where the newly interpolated values and all old ones serve as a basis for the new estimate.
Rationale:
Assume some function $\small f(x)$; for a simple example I used $\small f(x) = exp(x) $. I want to find the coordinate $\small f(x_z)=z=5 $ Meaningful bounds for the initial interval are for instance $\small X_2 = [1,3] $
I compute the vandermonde-matrix $\small Y_2$ as $\small Y_2=\underset{r=1..2, c=1..2}{\operatorname{matrix}} \left[ f(X_2[r])^c \right] $ ; then $\ [z,z^2] \cdot Y_2^{-1} \cdot X_2 = x_3 $ gives an estimate for $\small x_z$ in the value of $\small x_3 $ .
Now I append the new value $\small x_3$ at the vector $\small X_2 $ to make it $\small X_3 = \operatorname{concat}(X_2,x_3) $. I create the quadratic vandermondematrix $\small Y_3 $ in the obvious way and iterate with increased dimensions.(The secant-method analogue were simply not to increase the dimension but to put the new estimate as new first or second component into the X-vector, depending on the sign of $\small z - f(x_k) $ )
Empirically, if I'm near enough the true root, I get nearly exactly quadratic convergence, which is at least an improvement over the Müller's-method. Since the quadratic convergence with the Newton-method was in the wiki-article put as somehow a limit (for the cost of need of the derivative), this is then also interesting.
However, for not-so-well chosen initial points it seems at least as possibly erratic as again the Newton-method (as a simple example I just tried the function $\small f(x)=exp(x)-4*x^2 $ having some quadratic disturbance over the exponential). And, since it requires inversion of increasing sizes of vandermonde-matrices with the additional problem, that with the increase of iterations the rows become massively linear depending, it is surely not a method for the practical use.
But for the principal understanding of that matter my questions are:
1) how could I attack the problem to show/argue/prove, that this process has in fact such a rate of convergence?
2) Is that quadratic convergence really a limit here? I hoped, that the dynamic of the method would as well allow a dynamic increasing of the rate of convergence.
Example: (iterations downwards; new estimates just appended to matrix, rightmost column log of absolute error)
$\qquad \small \begin{matrix} k & x_k=estim(x_z) & f(x_k) & \ln(|z-f(x_k)|) \\ 1 & 0 & 1.00000000000 & 1.38629436112 & \text{initial low x} \\ 2 & 5 & 148.413159103 & 4.96572968897 & \text{initial high x} \\ \hline \\ 3 & 0.135673098126 & 1.14530742925 & 1.34929125565 \\ 4 & 3.64058189042 & 38.1140084505 & 3.49995640888 \\ 5 & 3.30992901642 & 27.3831816466 & 3.10830985751 \\ 6 & 2.93800666159 & 18.8781781845 & 2.63031769174 \\ 7 & 2.53192098262 & 12.5776443810 & 2.02520238367 \\ 8 & 2.12577394716 & 8.37938017666 & 1.21769231299 \\ 9 & 1.79911612022 & 6.04430266580 & 0.0433493572149 \\ 10 & 1.63882927856 & 5.14913777620 & -1.90288472781 \\ 11 & 1.61020090108 & 5.00381639897 & -5.56844797984 & \text{convergence-rate quadratic}\\ 12 & 1.60943843944 & 5.00000263503 & -12.8466176558 \\ 13 & 1.60943791243 & 5.00000000000 & -27.3915350585 \\ 14 & 1.60943791243 & 5.00000000000 & -56.4731042545 \\ 15 & 1.60943791243 & 5.00000000000 & -114.629475172 \end{matrix} $
The error term $e_n=|x_\star-x_n|$ for the root estimate is given by
$$e_{n+1}=\frac1{n!}\left|(f^{-1})^{(n)}(\xi_n)f'(\xi_n)^n\right|e_ne_{n-1}\cdots e_1$$
for some $\xi_n$ in the convex combination of $\{x_\star\}\cup\{x_k\}_{k\le n}$. If $(f^{-1})^{(n)}(\xi_n)f'(\xi_n)^n/n!$ is bounded in $n$, it can be expected that the convergence is quadratic, as you have experimentally found.
Unfortunately this type of method is quite unreliable, as you have pointed out. Often it is better to sacrifice this theoretical asymptotic speed for simplicity and faster initial speed (especially since most applications only go to double precision).
As far as I know, the only ways to surpass quadratic convergence in general are to perform several function evaluations per iteration. In the sense of Ostrowski's index, your method is then optimal.