Roots of Polynomials

O2scl

Contents

Introduction

Classes are provided for solving quadratic, cubic, and quartic equations as well as general polynomials. There is a standard nomenclature: classes which handle polynomials with real coefficients and real roots end with the suffix _real (quadratic_real, cubic_real and quartic_real), classes which handle real coefficients and complex roots end with the suffix _real_coeff ( quadratic_real_coeff, cubic_real_coeff, quartic_real_coeff, and poly_real_coeff), and classes which handle complex polynomials with complex coefficients end with the suffix _complex (quadratic_complex, cubic_complex, quartic_complex, and poly_complex). As a reminder, complex roots may not occur in conjugate pairs if the coefficients are not real. Most of these routines will return an error if the leading coefficient is zero.

In the public interfaces to the polynomial solvers, the complex type std::complex<double> is used.

For quadratics, quadratic_real_coeff_gsl2 is the best if the coefficients are real, while if the coefficients are complex, use quadratic_complex_std. Both of these classes work with multiprecision floating point types.

For cubics with real coefficients, cubic_real_coeff_gsl, cubic_real_coeff_gsl2, and cubic_real_coeff_cern are nearly equivalent. the GSL-based classes are a bit faster and more accurate, but the CERNLIB-based class is a bit better at identifying the nature of the roots. If the coefficients are complex, use cubic_complex_std, but keep in mind this class sometimes fails to properly identify the nature of the roots. The classes cubic_real_coeff_gsl2, cubic_real_coeff_cern, and cubic_complex_std all work with multiprecision floating point types.

Currently, quartics with real coefficients are best solved with poly_real_coeff_gsl, albeit a bit more slowly than direct solvers. The direct solvers are faster, but fail unpredictably, and are thus still experimental.

Polynomial solver example

This example shows how to find the roots of the second-, third-, fourth-, and fifth-order Chebyshev polynomials

\[\begin{split}\begin{eqnarray} &2x^2-1& \nonumber \\ &4x^3-3 x& \nonumber \\ &8x^4-8x^2+1& \nonumber \\ &16x^5-20x^3+5x& \nonumber \end{eqnarray}\end{split}\]

For the Chebyshev polynomial of order \(n\), the roots are given by

\[\cos \left[ \frac{\pi(k-1/2)}{n}\right]\]

for \(k = 1,\ldots,n\) These roots are used in cheb_approx_tl to approximate functions using Chebyshev polynomials .

/* Example: ex_poly.cpp
   -------------------------------------------------------------------
   Demonstrate the solution of the Chebyshev polynomials. See "License 
   Information" section of the documentation for license information.
*/

#include <boost/numeric/ublas/vector.hpp>
#include <o2scl/poly.h>
// For pi
#include <o2scl/constants.h>
#include <o2scl/vector.h>
#include <o2scl/test_mgr.h>

using namespace std;
using namespace o2scl;
using namespace o2scl_const;

int main(void) {

  cout.setf(ios::scientific);
  cout.setf(ios::showpos);

  test_mgr t;
  t.set_output_level(1);

  typedef boost::numeric::ublas::vector<double> ubvector;

  // Quadratic solver
  quadratic_real_coeff_gsl2<> quad;
  // Cubic solver
  cubic_real_coeff_cern<> cubic;
  // Quartic solver
  quartic_real_coeff_cern<> quart;
  // Generic polynomial solver
  poly_real_coeff_gsl<> gen;

  // Storage for the roots
  ubvector v(5);
  double d;
  std::vector<std::complex<double> > ca(5);

  // The second order polynomial
  cout << "Second order roots: " << endl;
  quad.solve_r(2.0,0.0,-1.0,v[0],v[1]);

  // Sort the roots and compare with the exact results
  vector_sort<ubvector,double>(2,v);
  for(size_t i=0;i<2;i++) {
    double exact=cos(pi*(((double)(2-i))-0.5)/2.0);
    cout << v[i] << " " << exact << endl;
    t.test_abs(v[i],exact,1.0e-14,"2nd order");
  }
  cout << endl;

  // The third order polynomial
  cout << "Third order roots: " << endl;
  cubic.solve_rc(4.0,0.0,-3.0,0.0,v[0],ca[0],ca[1]);

  // Sort the roots and compare with the exact results
  v[1]=ca[0].real();
  v[2]=ca[1].real();
  vector_sort<ubvector,double>(3,v);
  for(size_t i=0;i<3;i++) {
    double exact=cos(pi*(((double)(3-i))-0.5)/3.0);
    cout << v[i] << " " << exact << endl;
    if (i==1) {
      t.test_abs(v[i],exact,1.0e-14,"3rd order");
    } else {
      t.test_abs(v[i],exact,1.0e-14,"3rd order");
    }
  }
  cout << endl;

  // The fourth order polynomial
  cout << "Fourth order roots: " << endl;
  quart.solve_rc(8.0,0.0,-8.0,0.0,1.0,ca[0],ca[1],ca[2],ca[3]);

  // Sort the roots and compare with the exact results
  for(size_t i=0;i<4;i++) v[i]=ca[i].real();
  vector_sort<ubvector,double>(4,v);
  for(size_t i=0;i<4;i++) {
    double exact=cos(pi*(((double)(4-i))-0.5)/4.0);
    cout << v[i] << " " << exact << endl;
    t.test_abs(v[i],exact,1.0e-14,"4th order");
  }
  cout << endl;

  // The fifth order polynomial
  cout << "Fifth order roots: " << endl;
  vector<double> co={16.0,0.0,-20.0,0.0,5.0,0.0};
  gen.solve_rc_arr(5,co,ca);

  // Sort the roots and compare with the exact results
  for(size_t i=0;i<5;i++) v[i]=ca[i].real();
  vector_sort<ubvector,double>(5,v);
  for(size_t i=0;i<5;i++) {
    double exact=cos(pi*(((double)(5-i))-0.5)/5.0);
    cout << v[i] << " " << exact << endl;
    t.test_abs(v[i],exact,1.0e-14,"5th order");
  }
  cout << endl;

  t.report();
  return 0;
}