Bernoulli's method

This article is about the root-finding algorithm. For an unrelated topic in ordinary differential equation, see Bernoulli differential equation. For an unrelated topic in fluid dynamics, see Bernoulli's principle.

Excerpt from Observations Concerning Series by Daniel Bernoulli published in 1729.

In numerical analysis, Bernoulli's method, named after Daniel Bernoulli, is a root-finding algorithm which calculates the largest root (or zero) in magnitude of a polynomial. The method uses a quotient ratio of successive terms in a sequence calculated using only the coefficients of the given polynomial and does not require an initial starting root approximation. Since it converges slowly with a linear order, it can be used for finding an initial guess for another root-finding algorithm, such as Newton's method.

History

Lagrange on Bernoulli's method.

Bernoulli's method was first introduced by Swiss-French mathematician and physicist Daniel Bernoulli (1700-1782) in 1729. He noticed a trend from recurrent series created using polynomial coefficients growing by a ratio related to a root of the polynomial but did not prove why it worked. In 1725, Bernoulli moved to Saint Petersburg with his brother Nicolaus II Bernoulli, who unfortunately died of fever in 1726. While there, he worked closely with Leonhard Euler, a student of Johann Bernoulli, and made many advancements in harmonics, mathematical economics (see St. Petersburg paradox), and hydrodynamics. Euler called Bernoulli's method "frequently very useful" and gave a justification for why it works in 1748. The mathematician Joseph-Louis Lagrange expanded on this for the case of multiple roots in 1798. Bernoulli's method predates other root-finding algorithms like Graeffe's method (1826 to Dandelin) and is contemporary to Halley's method (1694). Since then, it has influenced the development of more modern algorithms such as the QD method.

The method

Given a polynomial p written as:

$p(z)=a_{0}z^{d}+a_{1}z^{d-1}+\dots +a_{d}$

where d is the degree of p. Pick d arbitrary starting values for a sequence x written as $x=\{x_{-d+1},x_{-d},\dots ,x_{-1},x_{0}\}$, where usually $\{0,0,0,\dots ,0,1\}$ is used for initial starting values. Then the subsequent terms in the sequence are computed with:

$x_{n}=-{\frac {a_{1}x_{n-1}+a_{2}x_{n-2}+\dots +a_{d}x_{n-d}}{a_{0}}}$

The sucessive ratios of two successive numbers in the sequence x, written as $q_{n}={\frac {x_{n+1}}{x_{n}}}$ in the limit tends to the dominant root, meaning the one with the largest magnitude, of p. For complex roots, this refers to the root with the largest complex modulus or distance from the origin in the complex plane.

Derivation of the method

Consider the n-th order difference equation:

$a_{n}x_{m}+a_{n-1}x_{m-1}+\dots +a_{0}x_{m-n}=0\quad (a_{n}\neq 0;\ m=n,n+1,\dots )$

which has a solution:

$x_{m}=c_{1}r_{1}^{m}+c_{2}r_{2}^{m}+\dots +c_{n}r_{n}^{m}$

where $r_{i}$ is a root of p and $c_{i}$ are suitable constants. By division of successive terms in this sequence we get:

${\frac {x_{m+1}}{x_{m}}}={\frac {c_{1}r_{1}^{m+1}+c_{2}r_{2}^{m+1}+\dots +c_{n}r_{n}^{m+1}}{c_{1}r_{1}^{m}+c_{2}r_{2}^{m}+\dots +c_{n}r_{n}^{m}}}$

Assuming $r_{1}$ is the dominant root such that $|r_{1}|>|r_{k}|,\ k=2,3,\dots ,n$ then factoring out $r_{1}$ gives:

${\frac {x_{m+1}}{x_{m}}}=r_{1}\cdot {\frac {c_{1}+c_{2}\left({\frac {r_{2}}{r_{1}}}\right)^{m+1}+\dots +c_{n}\left({\frac {r_{n}}{r_{1}}}\right)^{m+1}}{c_{1}+c_{2}\left({\frac {r_{2}}{r_{1}}}\right)^{m}+\dots +c_{n}\left({\frac {r_{n}}{r_{1}}}\right)^{m}}}$

Since each ${\frac {r_{k}}{r_{1}}}$ term with $k\geq 2$ has absolute value less than 1, then $\lim _{m\to \infty }\left({\frac {r_{k}}{r_{1}}}\right)^{m}=0$. Hence $\lim _{m\to \infty }{\frac {x_{m+1}}{x_{m}}}=r_{1}$. This assumes $c_{1}\neq 0$, which Blum shows is satisfied by using initial values of all zeros followed by a final 1.

Convergence

Bernoulli's method will converge on the largest root in magnitude of a given polynomial with a linear order of convergence. This can lead to problems when there is multiplicity of the dominant root, a dominant conjugate pair, or otherwise multiple equidistant dominant roots, although variations on Bernoulli's method exist to combat these issues. To speed convergence, Alexander Aitken developed his Aitken delta-squared process as part of an improvement on his extension to Bernoulli's method, which also found all of the roots simultaneously.

In the case of converging on the smallest root of a polynomial (instead of the largest), the coefficients of the original polynomial are reversed, the method applied on the new polynomial, and the result inverted, which would give the smallest root of the original. When using root deflation with something like Horner's method, deflating from the smallest root is more stable. Another extension of Bernoulli's method is the Quotient-Difference method, which also finds all roots simultaneously, even though it can be unstable. Given the slow convergence of Bernoulli's method, and the instability of QD method, they can instead be used as reliable ways to find initial starting values for other root-finding algorithms, rather than iterated until tolerance.

Example

Plot of the example polynomial on a Cartesian plane showing the dominant root at $\left({\frac {1+{\sqrt {5}}}{2}},0\right)$.

Bernoulli's method convergence on the root of $p(z)=z^{2}-z-1$.

Let $p(z)=z^{2}-z-1=0$. Then $a_{0}=1$, $a_{1}=-1$, and $a_{2}=-1$, so the recurrence becomes:

${\begin{aligned}x_{n}&=-{\frac {(-1)x_{n-1}+(-1)x_{n-2}}{1}}\\&=x_{n-1}+x_{n-2}\end{aligned}}$

Using the recommended initial values $x_{-1}=0$, $x_{0}=1$ generates the following table:

n	xn	qn
-1	0	−
0	1	1
1	1	2
2	2	1.5
3	3	1.666
4	5	1.6
5	8	1.625
6	13	1.61538461538
7	21	1.61904761905
8	34	1.61764705882
9	55	1.61818181818

This eventually converges on ${\frac {1+{\sqrt {5}}}{2}}\approx 1.618034$, also known as the Golden ratio, which is the largest root of the example polynomial. The sequence ${x_{n}}$ is also the well-known Fibonacci sequence. Bernoulli's method works even if the sequence used different starting values instead of 0 and 1; the limit of the quotient $q_{n}$ remains the same.

Comparison with other methods

The example above illustrates how Bernoulli's method generates a sequence converging on the dominant root. Compared to other root-finding algorithms, Bernoulli's method offers distinct advantages and limitations.

Advantages

• No initial guess: Newton's method, Secant method, Halley's method, and other similar approaches, all require one or more starting values. Bernoulli's method requires only the polynomial coefficients, eliminating the need for an initial guess.
• No derivatives: Although derivatives of polynomials are straightforward with the power rule, this is a computation that is not required in Bernoulli's method.
• Naturally finds a dominant root: Normally, finding large roots is considered less stable, but substituting z in p with $\left({\frac {1}{z}}\right)$, which reverses the order of coefficients, and then inverting the result of Bernoulli's method gives the smallest root of p, which is more stable.

Limitations

• Slow convergence: Fröberg writes "As a rule, Bernoulli's method converges slowly, so instead, one ought to use, for example, the Newton-Raphson method." This is in contrast to Jennings who writes "The approximate zeros obtained by the Bernoulli method can be further improved by applying, say, the Newton-Raphson method". One author argues for instead-of while the other promotes in-conjunction-with, due to the linear order of convergence. It is important to note that the method's slow convergence can be improved with Aitken's delta-squared process.
• Finds one root at a time: The standard version of Bernoulli's method finds a single root requiring deflation to find another. When compared to algorithms such as Durand–Kerner method, Aberth method, Bairstow's method, and the "RPOLY" version of Jenkins–Traub algorithm they find multiple roots by default. One can overcome this limitation by applying an extension of Bernoulli's method such as the method by Aitken or QD method.
• Issues with multiples: Multiplicity and multiple dominant roots, such as conjugate pairs, can exacerbate the slowness of Bernoulli's method, yet improvements can be made to counter this.

Modern applications

Polynomial root-finding algorithms each have their own niches to which they are suited and traditional Bernoulli's method can provide initial approximations for other algorithms. It can also be used to find complex roots yet the more sophisticated extensions of Bernoulli's method, such as the one by Aitken and QD method, are able to find complex roots while solving for all of the roots simultaneously. There are also variations on Bernoulli's method that improve stability and handle multiple roots.

In related applications, Bernoulli's method has been shown to be equivalent Power method on a companion matrix for finding eigenvalues. Advancements in systolic arrays have led to a parallelized version of Bernoulli's method. The method has also been generalized to find poles of rational functions, extending to the field of complex analysis. An implementation of Bernoulli's method is included with the CodeCogs open source numerical methods library.

See also

• Aitken's delta-squared process
• Graeffe's method
• Horner's method
• Lehmer-Schur algorithm
• List of things named after members of the Bernoulli family
• Polynomial root-finding