Hermite插值法(C++)

Hermite插值法是一种常用的数值方法, 主要用于在给定的数据点集上构造一个可微的函数来近似描述这些数据点的趋势和特征. 它可以通过已知的函数值和导数值来推断出未知的函数值和导数值, 并且可以用于多元插值问题. Hermite插值法构造的多项式是可微的, 这意味着它可以用来近似描述数据点的导数和曲率等特征, 对于需要对数据进行微积分分析的问题非常有用, 例如优化、控制和工程设计等领域. 此外, Hermite插值法可以通过增加数据点的函数值和导数来扩展到高次多项式, 从而提高插值的精度和可微性.

问题描述

给定函数$f(x)$在区间$[a,b]$上$n+1$个互异节点$x_0,x_1,\cdots ,x_n$处的函数值, 以及直到$m_i$阶的导数值$f(x_i),f^{\prime }(x_i),\cdots ,f^{(m_i)}(x_i)$. 令$m=\sum _{i=0}^n(m_i+1)-1$, 求一个次数不超过$m$的多项式$H_m(x)$, 使得
$$H_m^{(j)}(x_i)=f^{(j)}(x_i),i=0,1,\cdots ,n,j=0,\cdots ,m_i$$

求解方法

Lagrange型Hermite插值法

$m$次Lagrange插值多项式可以表示为:
$$H_m(x)=\sum _{i=0}^n\sum _{k=0}^{m_i}{f}^{(k)}(x_i)l_{ik}(x)$$

其中,
$$l_{ik}(x)=\frac{1}{k!}{\left\{\frac{1}{\prod _{\genfrac{}{}{0ex}{}{j=0}{j\ne i}}^n(x-x_j{)}^{m_j+1}}\right\}}_{x_i}^{m_i-k}(x-x_i{)}^k\prod _{\genfrac{}{}{0ex}{}{j=0}{j\ne i}}^n(x-x_j{)}^{m_j+1}$$
为Lagrange插值基函数, { φ ( x ) } x 0 k \{\varphi(x)\}_{x_0}^k {φ(x)}x0​k​表示函数$\phi (x)$在$x_0$处的$k$次泰勒展开.

由于Lagrange型Hermite插值法的表达式中出现了泰勒展开, 无法化简为一般的初等函数形式, 因此不便于实现一般的Lagrange型Hermite插值法.

Newton型Hermite插值法

Newton型Hermite插值法是一种更常用和易于实现的Hermite插值方法. 它使用Newton插值多项式的形式, 并考虑已知函数值和导数值来构造插值多项式. 以下是Newton型Hermite插值法的描述:

给定函数$f(x)$在区间$[a,b]$上$n+1$个互异节点$x_0,x_1,\dots ,x_n$处的函数值, 以及直到$m_i$阶的导数值$f(x_i),f^{\prime }(x_i),\dots ,f^{(m_i)}(x_i)$. 令$m=\sum _{i=0}^n(m_i+1)-1$, 我们希望找到一个次数不超过$m$的插值多项式$H_m(x)$, 满足以下条件:
$$H_m^{(j)}(x_i)=f^{(j)}(x_i),i=0,1,\dots ,n,j=0,1,\dots ,m_i$$

Newton型Hermite插值法的多项式表示如下:
$$\begin{array}{rl}{H}_m(x)=& f(x_0)+f[x_0,x_0](x-x_0)+\cdots +f[\underset{m_0+1}{\underbrace{x_0,x_0,\cdots ,x_0}}](x-x_0{)}^{m_0}+\\ & f[\underset{m_0+1}{\underbrace{x_0,x_0,\cdots ,x_0}},x_1](x-x_0{)}^{m_0+1}+\cdots +\\ & f[\underset{m_0+1}{\underbrace{x_0,x_0,\cdots ,x_0}},\underset{m_1+1}{\underbrace{x_1,x_1,\cdots ,x_1}}](x-x_0{)}^{m_0+1}(x-x_0{)}^{m_1}+\cdots +\\ & f[\underset{m_0+1}{\underbrace{x_0,x_0,\cdots ,x_0}},\cdots ,\underset{m_{n-1}+1}{\underbrace{x_{n-1},x_{n-1},\cdots ,x_{n-1}}},x_n]\prod _{i=0}^{n-1}(x-x_i{)}^{m_i+1}+\cdots +\\ & f[\underset{m_0+1}{\underbrace{x_0,x_0,\cdots ,x_0}},\cdots ,\underset{m_{n-1}+1}{\underbrace{x_{n-1},x_{n-1},\cdots ,x_{n-1}}},\underset{m_n+1}{\underbrace{x_n,x_n,\cdots ,x_n}}](x-x_n{)}^{m_n}\prod _{i=0}^{n-1}(x-x_i{)}^{m_i+1}\end{array}$$
其中, 上式中出现了重节点的差商, 对于这类特殊差商, 可以用如下极限的思想得出其值:
$$f\left[x_0,x_0\right]=\underset{x\to x_0}{\lim }f\left[x,x_0\right]=\underset{x\to x_0}{\lim }\frac{f(x_0)-f(x)}{x_0-x}=f^{\prime }(x_0)$$
$$f[\underset{k+1}{\underbrace{x_0,x_0,\cdots ,x_0}}]=\underset{x_i\to x_0}{\lim }f[x_0,x_1,\cdots ,x_k]=\underset{x_i\to x_0}{\lim }\frac{f^{(k)}(\eta )}{k!}=\frac{f^{(k)}(x_0)}{k!}$$
$$f\left[x_0,x_0,x_1\right]=\underset{x\to x_0}{\lim }f\left[x,x_0,x_1\right]=\underset{x\to x_0}{\lim }\frac{f\left[x_0,x_1\right]-f\left[x,x_0\right]}{x_1-x}=\frac{f\left[x_0,x_1\right]-f\left[x_0,x_0\right]}{x_1-x_0}$$

代码实现

本文使用Newton型Hermite插值多项式实现Hermite插值, 首先需要包含头文件:

#include <armadillo>
#include <vector>
using namespace arma;
using namespace std;

另外需要引入检查向量是否包含重复元素的函数, 见多项式插值法(基于c++,armadillo库):

namespace _MYFUNCTION {
    extern bool check_unique(const vec&);
}

首先定义差商表类:

// 差商表
class Difference_quotient_table
{
    double *x;                 // 观测值
    double *y;                 // 差商值
    unsigned *m;               // 当前导数阶数
    const mat *Y;              // 观测数据
    const vector<unsigned> *M; // 导数阶数
    unsigned k;                // 差商数
    unsigned step;             // 观测点的步长
    unsigned *pre_m;           // 当前导数阶数的前一个位置
    unsigned fact;             // 阶乘
    double *x_out;             // 输出当前x
    //    static constexpr double eps=1e-100;
public:
    Difference_quotient_table(const vec *x0, const mat *y0, const vector<unsigned> *M) : Y(y0), step(0), fact(1)
    {
        unsigned i(k = (this->M = M)->size()), *m_p(m = new unsigned[k] + k);
        auto M_p(M->cend());
        do
        {
            if ((*--m = *--M_p) > y0->n_cols)
            {
                delete[] (m -= --i);
                throw "导数数据缺失！";
            }
            k += *m;
        } while (--i);
        y = new double[k] + k;
        x = new double[k] + k;
        const double *y_p(&y0->at(M->size())), *x_p(&x0->at(M->size())); // armadillo库的矩阵是按列存储的
        do
        {
            *--x = *--x_p, *--y = *--y_p;
            if (i = *--m_p)
                do
                    *--x = *x_p, *--y = *y_p;
                while (--i);
        } while (m_p != m);
        pre_m = m - 1;
        --(x_out = x);
    }
    ~Difference_quotient_table() noexcept
    {
        delete[] y;
        delete[] x;
        delete[] m;
    }
    // 进行差商
    bool diff()
    {
        if (!--k)
            return false; // 只有1个元素, 不需要计算
        unsigned d_k(++step);
        fact *= step;
        unsigned *m_p(pre_m), i(0);
        double *p(y), *q(x), *q1(q + step);
        do
            if (*q != *q1) // abs(*p-=*(p+1))>eps
                (*p -= *(p + 1)) /= *q - *q1;
            else
            {
                while (!*++m_p)
                    ; // 找到第1个未完全计算的重节点差商
                unsigned j(*m_p);
                i += --*m_p;
                q += *m_p;
                q1 += *m_p;
                const double t(Y->at(m_p - m, d_k) / fact);
                *p = t;
                while (--j)
                    *++p = t;
            }
        while (++q, ++q1, ++p, ++i != k);
        return true;
    }
    double &get_diff() noexcept
    {
        return *y;
    }
    double &get_x() noexcept
    {
        return *++x_out;
    }
};

利用差商表可以很容易实现Newton型Hermite插值法如下:

/*
 * Hermite插值法
 * x0:观测点
 * y0:观测数据(第i列为i阶导数)
 * m :每个节点有效的导数阶数
 * x1:插值点
 * y1:预测值
 */
void Hermite(const vec &x0, const mat &y0, const vector<unsigned> &m, const vec &x1, vec &y1)
{
    if (x1.empty())
    {
        y1.clear();
        return;
    }
    if (x0.empty() || y0.empty())
        throw "观测集为空！";
    if (m.size() != x0.n_elem)
        throw "导数阶数向量与观测点向量大小不一致！";
    if (x0.n_elem != y0.n_rows)
        throw "观测点与观测数据向量大小不一致！";
    if (_MYFUNCTION::check_unique(x0))
        throw "观测集包含重复元素！";
    Difference_quotient_table D(&x0, &y0, &m);
    y1.set_size(x1.n_elem).fill(D.get_diff());
    vec x(ones(x1.n_elem));
    while (D.diff())
        y1 += D.get_diff() * (x %= x1 - D.get_x());
}

注意: 在此程序中观测点矩阵的第$i$列代表第$i-1$阶导数, 若输入的导数阶数大于等于观测点矩阵将弹出错误提示.

测试案例

求4次Newton型Hermite插值多项式$H_4(x)$, 使得$H_4(0)=1,H_4^{\prime }(0)=2,H_4^{\prime \prime }(0)=3,H_4(1)=4,H_4^{\prime }(1)=5$.

使用Newton型Hermite插值计算, 将观测点取为$(0,1{)}^T$, 观测数据取为
$$\left(\begin{array}{ccc}1& 2& 3\\ 4& 5& 0\end{array}\right)$$

其中, 观测数据的第2行第3列的数据是任取的.

#include <stdio.h>

// Hermite插值代码略

int main()
{
    // 正确答案
    auto ans = [](double x) -> double
    {
        return (((1.5 * x - 2) * x + 1.5) * x + 2) * x + 1;
    };
    vec x0{0., 1.}, x1(linspace(0, 1)), f_x1(x1), y1;
    f_x1.transform(ans);
    mat y0{{1, 2, 3}, {4, 5, 0}};
    vector<unsigned> m;
    m.push_back(2);
    m.push_back(1);
    try
    {
        Hermite(x0, y0, m, x1, y1);
        FILE *file;
        if (fopen_s(&file, "data.csv", "w"))
            throw "文件保存失败！";
        for (unsigned i(0); i != 101; ++i)
            fprintf(file, "%f,%f,%f\n", x1.at(i), y1.at(i), f_x1.at(i));
        fclose(file);
    }
    catch (const char *s)
    {
        printf(s);
    }
    return 0;
}

最终, 我们得到了运行结果, 并利用Origin画图如下: