polynomial.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2000 - 2018 by the deal.II authors
//
// This file is part of the deal.II library.
//
// The deal.II library is free software; you can use it, redistribute
// it, and/or modify it under the terms of the GNU Lesser General
// Public License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
// The full text of the license can be found in the file LICENSE at
// the top level of the deal.II distribution.
//
// ---------------------------------------------------------------------

#ifndef dealii_polynomial_h
#define dealii_polynomial_h



#include <deal.II/base/config.h>
#include <deal.II/base/exceptions.h>
#include <deal.II/base/subscriptor.h>
#include <deal.II/base/point.h>

#include <memory>
#include <vector>

DEAL_II_NAMESPACE_OPEN

namespace Polynomials
{

  template <typename number>
  class Polynomial : public Subscriptor
  {
  public:
    Polynomial (const std::vector<number> &coefficients);

    Polynomial (const unsigned int n);

    Polynomial (const std::vector<Point<1> > &lagrange_support_points,
                const unsigned int            evaluation_point);

    Polynomial ();

    number value (const number x) const;

    void value (const number         x,
                std::vector<number> &values) const;

    void value (const number         x,
                const unsigned int n_derivatives,
                number *values) const;

    unsigned int degree () const;

    void scale (const number factor);

    template <typename number2>
    void shift (const number2 offset);

    Polynomial<number> derivative () const;

    Polynomial<number> primitive () const;

    Polynomial<number> &operator *= (const double s);

    Polynomial<number> &operator *= (const Polynomial<number> &p);

    Polynomial<number> &operator += (const Polynomial<number> &p);

    Polynomial<number> &operator -= (const Polynomial<number> &p);

    bool operator == (const Polynomial<number> &p)  const;

    void print(std::ostream &out) const;

    template <class Archive>
    void serialize (Archive &ar, const unsigned int version);

  protected:

    static void scale (std::vector<number> &coefficients,
                       const number         factor);

    template <typename number2>
    static void shift (std::vector<number> &coefficients,
                       const number2        shift);

    static void multiply (std::vector<number> &coefficients,
                          const number factor);

    void transform_into_standard_form ();

    std::vector<number> coefficients;

    bool in_lagrange_product_form;

    std::vector<number> lagrange_support_points;

    number lagrange_weight;
  };


  template <typename number>
  class Monomial : public Polynomial<number>
  {
  public:
    Monomial(const unsigned int n,
             const double coefficient = 1.);

    static
    std::vector<Polynomial<number> >
    generate_complete_basis (const unsigned int degree);

  private:
    static std::vector<number> make_vector(unsigned int n,
                                           const double coefficient);
  };


  class LagrangeEquidistant: public Polynomial<double>
  {
  public:
    LagrangeEquidistant (const unsigned int n,
                         const unsigned int support_point);

    static
    std::vector<Polynomial<double> >
    generate_complete_basis (const unsigned int degree);

  private:

    static
    void
    compute_coefficients (const unsigned int n,
                          const unsigned int support_point,
                          std::vector<double> &a);
  };



  std::vector<Polynomial<double> >
  generate_complete_Lagrange_basis (const std::vector<Point<1> > &points);



  class Legendre : public Polynomial<double>
  {
  public:
    Legendre (const unsigned int p);

    static
    std::vector<Polynomial<double> >
    generate_complete_basis (const unsigned int degree);
  };

  class Lobatto : public Polynomial<double>
  {
  public:
    Lobatto (const unsigned int p = 0);

    static std::vector<Polynomial<double> >
    generate_complete_basis (const unsigned int p);

  private:
    std::vector<double> compute_coefficients (const unsigned int p);
  };



  class Hierarchical : public Polynomial<double>
  {
  public:
    Hierarchical (const unsigned int p);

    static
    std::vector<Polynomial<double> >
    generate_complete_basis (const unsigned int degree);

  private:
    static void compute_coefficients (const unsigned int p);

    static const std::vector<double> &
    get_coefficients (const unsigned int p);

    static std::vector<std::unique_ptr<const std::vector<double> > > recursive_coefficients;
  };



  class HermiteInterpolation : public Polynomial<double>
  {
  public:
    HermiteInterpolation (const unsigned int p);

    static std::vector<Polynomial<double> >
    generate_complete_basis (const unsigned int p);
  };



  class HermiteLikeInterpolation : public Polynomial<double>
  {
  public:
    HermiteLikeInterpolation (const unsigned int degree,
                              const unsigned int index);

    static std::vector<Polynomial<double> >
    generate_complete_basis (const unsigned int degree);
  };



  /*
   * Evaluate a Jacobi polynomial @f$ P_n^{\alpha, \beta}(x) @f$ specified by the
   * parameters @p alpha, @p beta, @p n, where @p n is the degree of the
   * Jacobi polynomial.
   *
   * @note The Jacobi polynomials are not orthonormal and are defined on the
   * unit interval @f$[0, 1]@f$ as usual for deal.II, rather than @f$[-1, +1]@f$ often
   * used in literature. @p x is the point of evaluation.
   */
  template <typename Number>
  Number
  jacobi_polynomial_value(const unsigned int degree,
                          const int alpha,
                          const int beta,
                          const Number x);


  template <typename Number>
  std::vector<Number>
  jacobi_polynomial_roots(const unsigned int degree,
                          const int alpha,
                          const int beta);
}


/* -------------------------- inline functions --------------------- */

namespace Polynomials
{
  template <typename number>
  inline
  Polynomial<number>::Polynomial ()
    :
    in_lagrange_product_form (false),
    lagrange_weight          (1.)
  {}



  template <typename number>
  inline
  unsigned int
  Polynomial<number>::degree () const
  {
    if (in_lagrange_product_form == true)
      {
        return lagrange_support_points.size();
      }
    else
      {
        Assert (coefficients.size()>0, ExcEmptyObject());
        return coefficients.size() - 1;
      }
  }



  template <typename number>
  inline
  number
  Polynomial<number>::value (const number x) const
  {
    if (in_lagrange_product_form == false)
      {
        Assert (coefficients.size() > 0, ExcEmptyObject());

        // Horner scheme
        const unsigned int m=coefficients.size();
        number value = coefficients.back();
        for (int k=m-2; k>=0; --k)
          value = value*x + coefficients[k];
        return value;
      }
    else
      {
        // direct evaluation of Lagrange polynomial
        const unsigned int m = lagrange_support_points.size();
        number value = 1.;
        for (unsigned int j=0; j<m; ++j)
          value *= x-lagrange_support_points[j];
        value *= lagrange_weight;
        return value;
      }
  }



  template <typename number>
  template <class Archive>
  inline
  void
  Polynomial<number>::serialize (Archive &ar, const unsigned int)
  {
    // forward to serialization function in the base class.
    ar &static_cast<Subscriptor &>(*this);
    ar &coefficients;
    ar &in_lagrange_product_form;
    ar &lagrange_support_points;
    ar &lagrange_weight;
  }



  template <typename Number>
  Number
  jacobi_polynomial_value(const unsigned int degree,
                          const int alpha,
                          const int beta,
                          const Number x)
  {
    Assert(alpha >=0 && beta >= 0,
           ExcNotImplemented("Negative alpha/beta coefficients not supported"));
    // the Jacobi polynomial is evaluated using a recursion formula.
    Number p0, p1;

    // The recursion formula is defined for the interval [-1, 1], so rescale
    // to that interval here
    const Number xeval = Number(-1) + 2. * x;

    // initial values P_0(x), P_1(x):
    p0 = 1.0;
    if (degree==0)
      return p0;
    p1 = ((alpha+beta+2)*xeval + (alpha-beta))/2;
    if (degree==1)
      return p1;

    for (unsigned int i=1; i<degree; ++i)
      {
        const Number v  = 2*i + (alpha + beta);
        const Number a1 = 2*(i+1)*(i + (alpha + beta + 1))*v;
        const Number a2 = (v + 1)*(alpha*alpha - beta*beta);
        const Number a3 = v*(v + 1)*(v + 2);
        const Number a4 = 2*(i+alpha)*(i+beta)*(v + 2);

        const Number pn = ((a2 + a3*xeval)*p1 - a4*p0)/a1;
        p0 = p1;
        p1 = pn;
      }
    return p1;
  }



  template <typename Number>
  std::vector<Number>
  jacobi_polynomial_roots(const unsigned int degree,
                          const int alpha,
                          const int beta)
  {
    std::vector<Number> x(degree, 0.5);

    // compute zeros with a Newton algorithm.

    // Set tolerance. For long double we might not always get the additional
    // precision in a run time environment (e.g. with valgrind), so we must
    // limit the tolerance to double. Since we do a Newton iteration, doing
    // one more iteration after the residual has indicated convergence will be
    // enough for all number types due to the quadratic convergence of
    // Newton's method

    const Number tolerance
      = 4 * std::max(static_cast<Number>(std::numeric_limits<double>::epsilon()),
                     std::numeric_limits<Number>::epsilon());

    // The following implementation follows closely the one given in the
    // appendix of the book by Karniadakis and Sherwin: Spectral/hp element
    // methods for computational fluid dynamics (Oxford University Press,
    // 2005)

    // If symmetric, we only need to compute the half of points
    const unsigned int n_points = (alpha == beta ? degree/2 : degree);
    for (unsigned int k=0; k<n_points; ++k)
      {
        // we take the zeros of the Chebyshev polynomial (alpha=beta=-0.5) as
        // initial values, corrected by the initial value
        Number r = 0.5-0.5*std::cos(static_cast<Number> (2*k+1)/(2*degree) *
                                    numbers::PI );
        if (k>0)
          r = (r + x[k-1])/2;

        unsigned int converged = numbers::invalid_unsigned_int;
        for (unsigned int it=1; it<1000; ++it)
          {
            Number s = 0.;
            for (unsigned int i=0; i<k; ++i)
              s += 1./(r - x[i]);

            // derivative of P_n^{alpha,beta}, rescaled to [0, 1]
            const Number J_x = (alpha+beta+degree+1)*
                               jacobi_polynomial_value(degree-1, alpha+1,
                                                       beta+1, r);

            // value of P_n^{alpha,beta}
            const Number f = jacobi_polynomial_value(degree, alpha, beta, r);
            const Number delta = f/(f*s - J_x);
            r += delta;
            if (converged == numbers::invalid_unsigned_int &&
                std::abs(delta) < tolerance)
              converged = it;

            // do one more iteration to ensure accuracy also for tighter
            // types than double (e.g. long double)
            if (it == converged + 1)
              break;
          }

        Assert(converged != numbers::invalid_unsigned_int,
               ExcMessage("Newton iteration for zero of Jacobi polynomial "
                          "did not converge."));

        x[k] = r;
      }

    // in case we assumed symmetry, fill up the missing values
    for (unsigned int k=n_points; k<degree; ++k)
      x[k] = 1.0-x[degree-k-1];

    return x;
  }

}
DEAL_II_NAMESPACE_CLOSE

#endif