numerics/ortho_poly.h

// -*- c++ -*-

// +--------------------------------------------------------------------+
// | This file is part of MathTL - the Mathematical Template Library    |
// |                                                                    |
// | Copyright (c) 2002-2009                                            |
// | Thorsten Raasch, Manuel Werner                                     |
// +--------------------------------------------------------------------+

#ifndef _MATHTL_ORTHO_POLY_H
#define _MATHTL_ORTHO_POLY_H

#include <cmath>
#include <utils/array1d.h>
#include <algebra/vector.h>
#include <algebra/polynomial.h>

namespace MathTL
{
  // some tools for polynomials that fulfill a three-term recurrence
  // relation (such as orthogonal polynomials, but also the monomials)
  //
  // reference: Deuflhard/Bornemann, Numerische Mathematik 1

  class OrthogonalPolynomial
  {
  public:
    virtual ~OrthogonalPolynomial() {}

    virtual double a(const unsigned int k) const = 0;

    virtual double b(const unsigned int k) const = 0;

    double operator () (const unsigned int n, const double x) const;

    Polynomial<double> assemble(const unsigned int n) const;

    double forwardSummation(const Vector<double>& coeffs, const double x) const;

    double adjointSummation(const Vector<double>& coeffs, const double x) const;
  };

  class Monomial : public OrthogonalPolynomial
  {
  public:
    double a(const unsigned int k) const;
    double b(const unsigned int k) const;
  };

  class ChebyshevPolynomial : public OrthogonalPolynomial
  {
  public:
    double a(const unsigned int k) const;
    double b(const unsigned int k) const;
  };

  class LegendrePolynomial : public OrthogonalPolynomial
  {
  public:
    double a(const unsigned int k) const;
    double b(const unsigned int k) const;
  };

  class GenMomentsPolynomial : public OrthogonalPolynomial
  {
  public:
    GenMomentsPolynomial(const Array1D<double>& moments,
                         const double a, const double b,
                         const unsigned int N);
    
    GenMomentsPolynomial(const Array1D<double>& moments,
                         const OrthogonalPolynomial& T,
                         const double a, const double b,
                         const unsigned int N);

    virtual ~GenMomentsPolynomial() {}

    double a(const unsigned int k) const;
    double b(const unsigned int k) const;
    
  private:
    Array1D<double> ak_, bk_;
  };
}

// include implementation of inline functions
#include <numerics/ortho_poly.cpp>

#endif