C++：Lagrange插值方法

提示：文章写完后，目录可以自动生成，如何生成可参考右边的帮助文档

---

前言

本篇文章是在学习计算方法后，尝试用Lagrange插值方法进行插值，其中遇到了一些问题，询问了一下别人得到了解决方法，因此写下了这篇文章。

---

提示：以下是本篇文章正文内容，下面案例可供参考

一、Lagrange形式

为了构造出Lagrange形式的插值公式，可以先作插值数据点
{($x_0$,0),($x_1$,0),…,($x_{i-1}$,0),($x_i$,0),($x_{i+1}$,0),…,($x_n$,0)}
的n次插值多项式为$l_i$(x)，由插值条件可知：
$$l_i(x_k)={\delta }_{ik}=\{\begin{array}{cc}0,& k\ne i\\ 1,& k=i\end{array}$$
$l_i(x)$是有$n$个根$x_0,x_1,...,x_{i-1},x_{i+1},...,x_n$的$n$次多项式，因此可以写成下面这种形式：
$$l_i(x)=a(x-x_0)...(x-x_{i-1})(x-x_{i+1})...(x-x_n)$$
其中，$a$是一个待定的常数，令$x=x_i$，代入上式并可以有插值条件可得：
$$1=l_i(x)=a(x-x_0)...(x-x_{i-1})(x-x_{i+1})...(x-x_n)$$
因此：
$$a={\left[\prod _{\begin{array}{c}k\ne i\end{array}}(x_i-x_k)\right]}^{-1}$$
从而可以得到：
$$l_i(x)=\frac{{\prod }_{\begin{array}{c}k\ne i\end{array}}(x-x_k)}{{\prod }_{\begin{array}{c}k\ne i\end{array}}(x_i-x_k)},i=0,1,2,...,n\text{\hspace{1em}(1)}$$
若采用记号：
$$\omega (x)=\prod _{k=1}^{n}x-x_k\text{\hspace{1em}(2)}$$
则
$$l_i(x)=\frac{\omega (x)}{(x-x_i){\omega }^{\prime }(x_i)}\text{\hspace{1em}(3)}$$
有式$(1)$或式$(3)$定义的$n$次多项式$l_i(x)$称为Lagrange基本插值多项式，由它们可以构成插值多项式的Lagrange形式：
$$L(x)=\sum _{i=0}^n{y}_i{l}_i(x)=\sum _{i=0}^n{y}_i\frac{\omega (x)}{(x-x_i){\omega }^{\prime }(x_i)}\text{\hspace{1em}(4)}$$
式$(4)$为满足下列插值条件
$$P(x_i)=p_i=y_i,i=0,1,2,\dots ,m\text{\hspace{1em}(5)}$$
的次数不超过$n$的多项式。
但是直接利用式$(4)$来计算给定的$\bar{x}$处的值是不经济的，大约需要$n^2$flop。若先计算出与求值点$\bar{x}$无关的量
$$d_i={\omega }^{\prime }(x_i),i=0,1,2,\dots ,n$$
则式$(4)$可以写成
$$L(x)=\omega (x)\sum _{i=0}^n\frac{y_i}{d_i(x-x_i)}\text{\hspace{1em}(6)}$$
如果所有的$d_i$都已经计算出，则可以按照式$(6)$求值$l(\bar{x})$约需要$3n$flop。因为$d_i$不依赖于求值点$\bar{x}$，它们可以预先计算出来，以备作插值时使用。

二、程序算法

代码如下（示例）：

double Lagrange(double X[], double Y[], int n, double x)//采用Lagrange插值方法，分别表示，长度，深度，数组数目和所需要计算点的值
{
    double result = 0.0;//最后得到的函数解
    double li[200];//lagrange基本插值多项式
    double temp1, temp2;//分别用来存储Lagrange多项式的分母和分子
    int i, j;//循环变量
    for (i = 0; i < n; i++)
    {
        temp1 = 1.0; temp2 = 1.0;
        for (j = 0; j < n; j++)
        {
            if (i == j)
                continue;
            temp1 = temp1 * (x - X[j]);
            temp2 = temp2 * (X[i] - X[j]);
        }
        li[i] = temp1 / temp2;
    }
    for (i = 0; i < n; i++)
    {
        result = result + Y[i] * li[i];
    }
    return result;
}

三、病态问题

因为Lagrange插值本身存在的问题，可能会出现插值得到的结果不正常，产生一个特别大的误差。此时如果想在这得到准确的解可以进行分段插值的方法，在此就不在赘述了，只是在这里表明存在着这么一个问题，如果是不想这样的话，建议直接使用分段样条插值的方法。

总结

此篇文章对Lagrange插值方法进行了总结，指出其存在着病态问题，最后在计算的时候应该十分的注意。