Вычисление производных высших порядков полинома Ньютона / Хабр

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

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

$P_n(x) = c_1 + c_2(x - x_1) + c_3(x - x_1)(x - x_2) + \ldots + c_n(x - x_1)(x - x_2) \ldots (x - x_{n - 1})$ ,

где:

— полином Ньютона n-й степени
$c_i = \frac{\Delta^iy_0}{i!h^i}$ — разделенная разность
— шаг между узлами
— значение функции в узле
$i = \overline{1, n}$

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

Пусть $(x - x_1) \ldots (x - x_n) = \prod_{i = 1}^n (x - x_i) = g_n(x),$ тогда:

$g_n'(x) = \sum_{i=1}^{n} \prod_{j=1,(j \ne i)}^{n}(x - x_j)$ $g_n''(x) = 2! \cdot \sum_{i_1 = 1}^{n} \sum_{i_2 = i_1 + 1}^{n} \prod_{j = 1, j \notin \{i_1, i_2\}}^n (x - x_j)$ $g_n'''(x) = 3! \cdot \sum_{i_1 = 1}^{n} \sum_{i_2 = i_1 + 1}^{n} \sum_{i_3 = i_2 + 1}^{n}\prod_{j = 1, j \notin \{i_1, i_2, i_3\}}^n (x - x_j)$ $\ldots$ $g_n^{(m)}(x) = m! \cdot \sum_{i_1 = 1}^{n} ( \sum_{i_2 = i_1 + 1}^{n} \ldots ( \sum_{i_m = i_{m-1}+1}^{n} \prod_{j = 1, j \notin \{i_1, i_2, \ldots,i_m\}}^n (x - x_j)))$

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

$P_n'(x) = c_2 + c_3 g_2'(x) + \ldots + c_n g_{n-1}'(x)$ $P_n''(x) = 2c_3 + c_4 g_3''(x) + \ldots + c_n g_{n-1}''(x)$ $\ldots$ $P_n^{(m)}(x) = m! \cdot c_{m + 1} + c_{m + 2} g_{m + 1}^{(m)}(x) + \ldots + c_n g_{n-1}^{(m)}(x),$

где:

$P_n^{(m)}$ — производная m-й степени полинома Ньютона
$m = \overline{0, n}$ , иначе производная равна 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;
}