Интерполяция  — задача восстановления функции по заданному дискретному набору её значений. При этом предполагается, что исходная функция является непрерывной на рассматриваемом отрезке. Одним из методов её решения является интерполяционный полином Ньютона, применяемый в случае равномерной сетки узлов. Существует две формы записи данного полинома: прямая и обратная.

В данной статье будет рассмотрен прямой полином Ньютона:

,

где:

• — полином Ньютона n-й степени

• — разделенная разность

• — шаг между узлами

• — значение функции в узле

• 

Поскольку интерполяционный полином Ньютона аппроксимирует исходную функцию, его дифференцирование позволяет получить приближённые значения производных интерполируемой функции. Далее выведем общую аналитическую формулу m-й производной полинома Ньютона, позволяющей вычислять производные произвольного порядка.

Пусть тогда:

Следовательно:

где:

• — производная m-й степени полинома Ньютона

• , иначе производная равна 0

Далее рассмотрим программную реализацию полинома Ньютона

Для начала рассмотрим реализацию класса таблицы узлов, которая используется в классе полинома Ньютона:

struct Point {
    double x;
    double y;
};

class PointsTable {
    std::vector<Point> data;

public:
    explicit PointsTable(void) : data({}) {}
    explicit PointsTable(const std::initializer_list<Point>& points) : data(points) {}

    std::size_t size(void) const;
  
    Point operator[](std::size_t index) const;

    std::vector<double> get_ys(void) const;
    std::vector<double> get_xs(void) const;
};


std::vector<double> PointsTable::get_ys(void) const {
    std::vector<double> y_data;
    for (std::size_t i = 0; i < (this->data).size(); i++)
        y_data.push_back((this->data)[i].y);
    return y_data;
}

std::vector<double> PointsTable::get_xs(void) const {
    std::vector<double> x_data;
    for (std::size_t i = 0; i < data.size(); i++)
        x_data.push_back(data[i].x);
    return x_data;
}

std::size_t PointsTable::size(void) const { return data.size(); }
Point PointsTable::operator[](std::size_t index) const { return data[index]; }

Реализация класса полинома Ньютона:

// Класс полинома Ньютона
class NewtonPolinom {
private:
    std::vector<double> _x; // Узлы интерполяции (x_1, x_2, ..., x_n)
    std::vector<double> coefs; // Коэффициенты полинома (c_1, c_2, ..., c_n)

    std::size_t d_degree = 0; // Порядок производной (0 - сам полином)
    // Функция для вычисления члена полинома g^{(m)}(x)
    // Параметры: x - точка вычисления, n - количество множителей, ignore_index - индексы пропускаемых узлов
    std::function<double(double, std::size_t, std::vector<std::size_t>)> term_calc;
public:
    explicit NewtonPolinom(const PointsTable table); 

    NewtonPolinom diff(void) const; // Вычисление производной от полинома
    std::size_t diff_degree(void) const; // Получение порядка текущей производной

    double operator()(double x) const;  // Вычисление значения полинома (или его производной) в точке x

private:
    void divided_differences(const PointsTable &table); // Вычисление коэффициентов c
    std::function<double(double, std::size_t, std::vector<std::size_t>)> term_calc_init(void); // Инициализация term_calc
};

int factorial(int n) { return !abs(n) ? 1 : n * factorial(n - 1); }

std::function<double(double, std::size_t, std::vector<std::size_t>)> NewtonPolinom::term_calc_init(void) {
    auto __x = _x;
    auto in_vec = [](std::vector<std::size_t> vec, std::size_t i) -> bool {
        for (auto _i : vec)
            if (_i == i)
                return true;
        return false;
    };

    return [__x, in_vec](double x, std::size_t n, std::vector<std::size_t> ignore_index) -> double {
        double term = 1;
        for (std::size_t i = 0; i < n; ++i)
            if (!ignore_index.size() || !in_vec(ignore_index, i))
                term *= (x - __x[i]);
        return term;
    };
}

NewtonPolinom::NewtonPolinom(const PointsTable table) {
    divided_differences(table);
    term_calc = term_calc_init();
}

NewtonPolinom NewtonPolinom::diff(void) const {
    NewtonPolinom _new = *this;
    ++_new.d_degree;

    auto __x = _x;
    auto prev_term = term_calc;
    std::size_t _prev_d = this->d_degree;

    auto in_vec = [](std::vector<std::size_t> vec, std::size_t i) -> bool {
        for (auto _i : vec)
            if (_i == i)
                return true;
        return false;
    };

    _new.term_calc = [__x, prev_term, _prev_d, in_vec](double x, std::size_t n, std::vector<std::size_t> ignore_index) -> double {
        double sum = 0;
        auto new_ignore = ignore_index;
        new_ignore.resize(ignore_index.size() + 1);
        for (std::size_t i = 0; i < n; ++i) {
            if (!ignore_index.size() || !in_vec(ignore_index, i)) {
                new_ignore[new_ignore.size() - 1] = i;
                sum += prev_term(x, n, new_ignore);
            }
        }
        return sum;
    };
    return _new;
};
std::size_t NewtonPolinom::diff_degree(void) const { return d_degree; };

double NewtonPolinom::operator()(double x) const {
    double result = 0;
    std::size_t n = coefs.size();
    double term = factorial(d_degree);
    for (std::size_t i = d_degree; i < n; i++) {
        result += coefs[i] * term;
        term = term_calc(x, i + 1, std::vector<std::size_t>());
    }

    return result;
}

void NewtonPolinom::divided_differences(const PointsTable &table) {
    auto y_data = table.get_ys();
    std::vector<double> result = { table[0].y };

    std::size_t k = 1, n = table.size() - 1;
    for (; y_data.size() != 1 && k <= n; k++) {
        std::vector<double> data;
        for (std::size_t i = 0; i < y_data.size() - 1; i++) {
            data.push_back((y_data[i] - y_data[i + 1]) / (table[i].x - table[i + k].x));
        }
        result.push_back((y_data = data)[0]);
    }
    _x = table.get_xs();
    coefs = result;
}