polynomial.h

Go to the documentation of this file.

 
 
#ifndef THEORETICA_INTERP_POLYNOMIAL_H
#define THEORETICA_INTERP_POLYNOMIAL_H
 
#include <vector>
#include "../core/real_analysis.h"
#include "../polynomial/polynomial.h"
#include "../algebra/algebra.h"
#include "../core/function.h"
 
 
namespace theoretica {
 
 
    template <
        typename Vector1, typename Vector2,
        enable_vector<Vector1> = true,
        enable_vector<Vector2> = true,
        typename Type = vector_element_t<Vector2>
    >
    inline polynomial<Type> lagrange(const Vector1& x, const Vector2& y) {
 
        if(!x.size()) {
            TH_MATH_ERROR("lagrange", x.size(), MathError::InvalidArgument);
            return polynomial<Type>(make_error<Type>());
        }
 
        if(x.size() != y.size()) {
            TH_MATH_ERROR("lagrange", y.size(), MathError::InvalidArgument);
            return polynomial<Type>(make_error<Type>());
        }
 
        // Check that all x_i are different to prevent division by zero
        for (unsigned int i = 0; i < x.size() - 1; ++i) {
 
            if(abs(x[i] - x[i + 1]) < MACH_EPSILON) {
                TH_MATH_ERROR("lagrange", x[i], MathError::InvalidArgument);
                return polynomial<Type>(make_error<Type>());
            }
        }
 
        // Lagrange polynomial to construct
        polynomial<Type> L = {0};
 
        for (unsigned int j = 0; j < x.size(); ++j) {
            
            // The Lagrange polynomial is a linear combination of all l_j
            polynomial<real> l_j = {1.0};
 
            for (unsigned int m = 0; m < x.size(); ++m) {
                
                if(m == j)
                    continue;
 
                // l_j = product(x - x_m / x_j - x_m)
                l_j *= polynomial<real>({-x[m], 1.0}) / (x[j] - x[m]);
            }
 
            // L = sum(y_j * l_j)
            polynomial<Type> term;
            term.resize(l_j.size());
 
            for (size_t i = 0; i < l_j.size(); ++i)
                term[i] = y[j] * l_j[i];
 
            L += term;
        }
 
        return L;
    }
 
 
    template<typename Vector = vec<real>>
    inline Vector chebyshev_nodes(real a, real b, unsigned int n) {
 
        Vector nodes;
        nodes.resize(n);
 
        real m = (b + a) / 2.0;
        real c = (b - a) / 2.0;
 
        for (unsigned int i = 1; i < n + 1; ++i)
            nodes[i - 1] = m + c * cos(real(2 * i - 1) / real(2 * n) * PI);
 
        return nodes;
    }
 
 
    template <
        typename RealFunction,
        typename Type = return_type_t<RealFunction>
    >
    inline polynomial<Type> interpolate_grid(RealFunction f, real a, real b, unsigned int order) {
 
        vec<real> x (order + 1);
        vec<Type> y (order + 1);
 
        // Sample <order + 1> equidistant points
        for (unsigned int i = 0; i < order + 1; ++i) {
            x[i] = (b - a) / real(order) * real(i) + a;
            y[i] = f(x[i]);
        }
 
        return lagrange(x, y);
    }
 
 
    template <
        typename RealFunction,
        typename Type = return_type_t<RealFunction>
    >
    inline polynomial<Type> interpolate_chebyshev(RealFunction f, real a, real b, unsigned int order) {
 
        vec<real> x = chebyshev_nodes(a, b, order + 1);
        vec<Type> y (order + 1);
 
        // Sample <order + 1> equidistant points
        for (unsigned int i = 0; i < order + 1; ++i)
            y[i] = f(x[i]);
 
        return lagrange(x, y);
    }
 
 
    template <
        typename RealFunction,
        typename Type = return_type_t<RealFunction>
    >
    inline polynomial<Type> interpolate(RealFunction f, real a, real b, unsigned int order) {
        return interpolate_chebyshev(f, a, b, order);
    }
 
}
 
#endif