多项式插值法(C++)

多项式插值法是一种通过已知的若干点找到一个多项式, 使得该多项式的图像能够经过这些点的方法. 这个方法在数学和科学领域中被广泛使用, 可以用于拟合实验数据, 近似复杂函数等.

具体来说, 给定一组插值点, 多项式插值法就是要找到一个多项式, 使得这个多项式的零点或极值与给定的插值点相符. 这个多项式的次数通常不会超过给定点的数目. 在数学上, 这就是一个求解插值多项式的问题.

常用的多项式插值法有直接法, 拉格朗日插值法和牛顿插值法等. 直接法就是根据插值点直接求解线性方程组来构造插值多项式. 拉格朗日插值法是通过求解拉格朗日插值多项式来找到插值点的一致多项式. 牛顿插值法则是利用牛顿插值公式来求解插值多项式.

在理论上可以证明插值多项式的存在唯一性, 即设$f(x)$在$[a,b]$上有定义, 且已知$f(x)$在$[a,b]$上$n+1$个互异节点$x_0,x_1,\cdots ,x_n$的函数值$f(x_0),f(x_1),\cdots ,f(x_n)$, 则存在一个唯一的不超过$n$次的多项式$P_n(x)$, 满足

$$P_n(x_i)=f(x_i),i=0,1,\cdots ,n$$

总的来说, 多项式插值法是一种通过已知的若干点找到一个多项式, 使得该多项式的图像能够经过这些点的方法. 它被广泛应用于各种科学和工程领域, 如数据拟合, 图像处理, 数值分析和控制系统设计等.

本文代码要包含的头文件如下:

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

Lagrange插值法

介绍

Lagrange插值法是一种通过已知的若干点找到一个多项式, 使得该多项式的图像能够经过这些点的方法. 这个方法最早由英国数学家爱德华·华林于1779年发现, 不久后由莱昂哈德·欧拉在1783年再次发现, 并在1795年由法国数学家约瑟夫·路易斯·拉格朗日发表在他的著作《师范学校数学基础教程》中, 因此被称为拉格朗日插值法.

Lagrange插值法的核心思想是利用已知的插值节点构造一个多项式函数, 使得该函数在这些节点上的值等于已知的值, 而在其他节点上的值等于0, 即Lagrange插值基函数

$$l_k(x)=\prod _{\begin{array}{c}i=0\\ i\ne k\end{array}}^n\frac{x-x_i}{x_k-x_i}=\frac{{\omega }_{n+1}(x)}{(x-x_k){\omega }_{n+1}^{\prime }(x_k)}$$

其中,

$${\omega }_{n+1}(x)=\prod _{i=0}^n(x-x_i)$$

而Lagrange插值多项式即为
$$L_n(x)=\sum _{k=0}^n{y}_k{l}_k(x)$$

其截断误差为
$$R_n(x)=\frac{f^{(n+1)}(\xi )}{(n+1)!}{\omega }_{n+1}(x)$$

在具体实现中, Lagrange插值法需要计算每个节点的权重系数, 然后通过线性组合这些节点的多项式函数来构造Lagrange插值多项式. 在计算权重系数时, 需要使用到每个节点的相邻节点的信息, 因此需要计算每个节点的前驱节点和后继节点的数量和位置.

Lagrange插值法的优点是它具有很好的数值稳定性, 即当插值节点的数量增加时, 插值误差通常会迅速减小. 此外, Lagrange插值法还可以处理多个自变量的情况, 即可以用于多元插值问题. 但是, 当插值节点数量增加时, 计算量也会相应增加, 因此需要考虑计算效率和精度之间的平衡.

代码实现

首先编写一个检查向量是否包含重复元素的函数:

namespace _MYFUNCTION
{
    /*
     * 检查向量是否包含重复元素
     *  包含  :true
     *  不包含:false
     */
    bool check_unique(const vec &x)
    {
        unordered_set<double> s;
        for (const auto &i : x)
            if (!s.insert(i).second)
                return true;
        return false;
    }
}

接着可以实现Lagrange插值法:

/*
 * Lagrange插值
 * x0:观测点
 * y0:观测数据
 * x1:插值点
 * y1:预测值
 */
void Lagrange(const vec &x0, const vec &y0, const vec &x1, vec &y1)
{
    if (x1.empty())
    {
        y1.clear();
        return;
    }
    if (x0.empty() || y0.empty())
        throw "观测集为空！";
    if (x0.n_elem != y0.n_elem)
        throw "观测点与观测数据向量大小不一致！";
    if (x0.n_elem < 2)
        throw "观测点数量小于2！";
    if (_MYFUNCTION::check_unique(x0))
        throw "观测集包含重复元素！";
    y1.zeros(x1.n_elem);
    const unsigned start(-1);
    for (unsigned i(0); i != x1.n_elem; ++i)
        for (unsigned j(0); j != x0.n_elem; ++j)
        {
            double t(y0.at(j));
            unsigned k(start);
            while (++k != j)
                (t *= x1.at(i) - x0.at(k)) /= x0.at(j) - x0.at(k);
            while (++k != x0.n_elem)
                (t *= x1.at(i) - x0.at(k)) /= x0.at(j) - x0.at(k);
            y1.at(i) += t;
        }
}

Newton插值法

介绍

Newton插值法是一种代数插值方法, 可以用于在给定的一组插值节点处构造一个多项式函数, 使得该函数在这些节点上的值等于已知的值. 这个方法最早由英国数学家艾萨克·牛顿在17世纪提出, 因此被称为Newton插值法.

Newton插值法引入了差商的概念, 用于在插值节点增加时便于计算. 差商可以理解为两个函数值的比值, 它能够描述给定节点处函数的局部行为. 通过差商, 可以将一个复杂的多项式插值问题转化为求解一组线性方程组的问题, 大大简化了计算. 其算法可表示为如下公式:
$$\begin{array}{rl}{N}_n\left(x\right)=& f(x_0)+f\left[x,x_0\right]\left(x-x_0\right)+f\left[x_0,x_1,x_2\right]\left(x-x_0\right)\left(x-x_1\right)+\\ & +\cdots +f\left[x_0,x_1,\cdots ,x_n\right]\left(x-x_0\right)\left(x-x_1\right)\cdots \left(x-x_{n-1}\right)\end{array}$$

其中,

$$\{\begin{array}{l}f\left[x_0,x_k\right]=\frac{f(x_k)-f(x_0)}{x_k-x_0}\\ f[x_0,x_1,\dots ,x_k]=\frac{f[x_1,x_2,\cdots ,x_k]-f[x_0,x_1,\cdots ,x_{k-1}]}{x_k-x_0}\end{array}$$

其截断误差为

$$R_n(x)=f(x)-N_n(x)=f\left[x,x_0,x_1,\cdots ,x_n\right]{\omega }_{n+1}(x)$$

在具体实现中, 首先需要计算差商表和余项表. 差商表的计算是从第一个节点开始, 逐步向后计算, 直到最后一个节点. 每个节点的差商都是根据其前驱节点的差商和当前节点的函数值来计算的. 而余项表的计算则是根据差商表和已知的函数值来逐步向后计算的.

Newton插值法的优点是它具有很好的数值稳定性和高效性. 当插值节点的数量增加时, 插值误差通常会迅速减小. 此外, Newton插值法还可以处理多个自变量的情况, 即可以用于多元插值问题.

代码实现

/*
 * Newton插值
 * x0:观测点
 * y0:观测数据
 * x1:插值点
 * y1:预测值
 */
void Newton(const vec &x0, const vec &y0, const vec &x1, vec &y1)
{
    if (x1.empty())
    {
        y1.clear();
        return;
    }
    if (x0.empty() || y0.empty())
        throw "观测集为空！";
    if (x0.n_elem != y0.n_elem)
        throw "观测点与观测数据向量大小不一致！";
    if (x0.n_elem < 2)
        throw "观测点数量小于2！";
    if (_MYFUNCTION::check_unique(x0))
        throw "观测集包含重复元素！";
    unsigned k(x0.n_elem), l(k), g(0);
    double *a(new double[k]), *p(a + k);
    auto q(y0.cend());
    do
        *--p = *--q;
    while (--l);
    y1.set_size(x1.n_elem).fill(*p);
    vec xt(ones(x1.n_elem));
    while (l = --k)
    {
        xt %= x1 - x0.at(g);
        p = a;
        q = x0.cbegin();
        y1 += ((*p -= *(p + 1)) /= *q - *(q + (++g))) * xt;
        while (--l)
        {
            ++q, ++p;
            (*p -= *(p + 1)) /= *q - *(q + g);
        }
    }
    delete[] a;
}

测试案例

设函数$f(x)=\sqrt{x}$, 并已知$f(x)$在$2,2.1,2.2$处的函数值, 试用2次插值多项式$P_2(x)$计算$f(2.15)$的近似值.

测试代码如下:

#include <iostream>
#include <math.h>
// 插值函数代码略
int main()
{
    double (*p)(double) = sqrt;
    vec x0{2.0, 2.1, 2.2}, y0(x0), x1(1), y1;
    y0.transform(p);
    x1.at(0) = 2.15;
    try
    {
        Newton(x0, y0, x1, y1);
        std::cout << y1.at(0) << std::endl;
        Lagrange(x0, y0, x1, y1);
        std::cout << y1.at(0) << std::endl;
    }
    catch (const char *s)
    {
        std::cout << s;
    }
    return 0;
}

求解得$P_2(2.15)=1.46629$.