Bernoulli — Rootfinding

Calculates the zeros of a polynomial using Bernoulli's algorithm.

Interface

#include <codecogs/maths/rootfinding/bernoulli.h>

using namespace Maths::Rootfinding;

std::vector<double> bernoulli(const std::vector<double> &polynomial)

Bernoulli's method exploits the connection between a linear difference equation and the zeros of its characteristic polynomial in order to find the zeros of a polynomial without knowing crude first approximations. Given the polynomial

$$p(z) = a_0 z^k + a_1 z^{k - 1} + \ldots + a_k \qquad k \geq 1 \qquad a_0 a_k \neq 0$$

(1)

the difference equation which has p as its characteristic polynomial is given as

$$a_0 x_n + a_1 x_{n - 1} + \ldots + a_k x_{n - k} = 0$$

(2)

Given any starting values $x_0, x_1, \ldots , x_{k - 1}$, the corresponding solution of the previous equation can be found numerically using the recurrence relation

$$x_n = - \frac {1} {a_0} (a_1 x_{n - 1} + a_2 x_{n - 2} + \ldots + a_k x_{n - k}) \qquad n \geq k$$

(3)

If we suppose that the zeros $z_1, z_2, \ldots , z_k$ of p all have multiplicity 1 (i.e. distinct zeros) then the solution can be expressed in the form

$$x_n = c_1 z_1^n + c_2 z_2^n + \ldots + c_k z_k^n$$

(4)

The $c_j$ coefficients can be computed if the zeros are known. But since these are not known, the coefficients are unknown. However, the quotients

$$q_n = \frac {x_{n + 1}} {x_n}$$

(5)

using the previous equation, are analytically represented by

$$q_n = \frac {c_1 z_1^{n + 1} + c_2 z_2^{n + 1} + \ldots + c_k z_k^{n + 1}} {c_1 z_1^n + c_2 z_2^n + \ldots c_k z_k^n}$$

(6)

If a polynomial p of degree K has zeros $z_1, z_2, \ldots , z_k$, not necessarily distinct, then the zero $z_j$ is called dominant if its modulus is strictly greater than the moduli of the other zeros, i.e. $|z_j| > |z_i|$ for all $i \neq j$. Now suppose the above polynomial p has a single dominant zero, and let this be $z_1$ for our case (this zero will be real if the coefficients of p are real). Also, assume that the starting values of $x_0, x_1, \ldots , x_{k - 1}$ of the solution ${x_n}$ are chosen so that $c_1 \neq 0$. It can be shown that this condition is always satisfied if the starting values are chosen so that

$$x_{-k + 1} = x_{-k + 2} = \ldots = x_{-1} = 0 \qquad x_0 = 1$$

(7)

Under these assumptions, we may write the quotients as

$$q_n = z_1 \frac {1 + \frac{c_2}{c_1} \left( \frac {z_2} {z_1} \right) ^ {n + 1} + \ldots + \frac {c_k} {c_1} \left( \frac {z_k} {z_1} \right) ^ {n + 1} } { 1 + \frac {c_2} {c_1} \left( \frac {z_2} {z_1} \right) ^ n + \ldots + \frac {c_k} {c_1} \left( \frac {z_k} {z_1} \right) ^ n }$$

(8)

Furthermore,

$$\lim_{n \rightarrow \infty} \left( \frac {z_j} {z_1} \right) ^ n = 0 \qquad j = 2, 3, \ldots , k$$

(9)

and it follows that

$$\lim_{n \rightarrow \infty} q_n = z_1$$

(10)

Having obtained this zero, the polynomial is deflated and the procedure is repeated. Thus this method can only furnish one or two zeros of a given polynomial at a time, and these zeros are those of largest or smallest magnitude. So, if a zero of intermediate modulus is desired, it is necessary to compute all larger (or all smaller) zeros and then remove them from the polynomial by deflation. Also if $z_2$ has nearly the same magnitude as $z_1$ the convergence process is very slow. However, if the zero of largest or smallest magnitude is the zero that is desired and is distinct, Bernoulli's method can be useful.

This algorithm finds the roots of the polynomial given as an array of coefficients and returns them wrapped inside a C++ vector object. The method described above has been used.

References

Jean-Pierre Moreau's Home Page, http://perso.wanadoo.fr/jean-pierre.moreau/
Claude Nowakowski, "Methodes de calcul numerique", PSI Editions, France, 1981
Wankere R. Mekwi, "Iterative Methods for Roots of Polynomials", Exeter College, University of Oxford

Example 1

#include <codecogs/maths/rootfinding/bernoulli.h>

#include <iostream>
#include <iomanip>

int main()  
{
    double A[5] = { 1, 2, -13, -14, 24 };

    std::vector<double> poly(A, A + 5), roots = Maths::Rootfinding::bernoulli(poly);

    std::cout << "The roots of the polynomial are : " << std::endl << std::endl;
    int i;
    for (i = 0; i < roots.size(); i++)
    {
      std::cout << std::setw(5) << "x_" << i << " = ";
      std::cout << std::setprecision(12) << roots[i] << std::endl;
    }
    return 0;
}

Output:

The roots of the polynomial are :

   x_0 = -3.99998329472
   x_1 = 3.00000133652
   x_2 = -2.00001264664
   x_3 = 0.999994604828

Parameters

polynomial

the coefficients of the polynomial given as a C++ vector object