COST (GBP)

1.49

0.00

Maths › Rootfinding ›

Bernoulli

viewed 2879 times and licensed 114 times

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

Controller: CodeCogs

Contents

Interface
Bernoulli
Page Comments

Interface

C++

Bernoulli

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$

(2)

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$

(3)

Given any starting values $\inline 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$

(4)

If we suppose that the zeros $\inline 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$

(5)

The $\inline 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}$

(6)

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}$

(7)

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

(8)

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 }$

(9)

Furthermore,

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

(10)

and it follows that

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

(11)

Having obtained this zero, the polynomial is deflated and the procedure is repeated. Thus this method can only furnishone 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 $\inline z_2$ has nearly the same magnitude as $\inline 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

Authors

Lucian Bentea (August 2005)

Source Code

Source code is available when you agree to a GP Licence or buy a Commercial Licence.

Not a member, then Register with CodeCogs. Already a Member, then Login.