拉格朗日插值法

任意一个有穷数列都有通项公式。

给出 $n$ 个点 $(x_1,y_1),(x_2,y_2),\cdots,(x_n,y_n)$，怎么找到一个函数 $f(x)$ 同时过这 $n$ 个点？

拉格朗日的想法是这样的：构造一个开关函数

$$L_k(t)=\begin{cases} 1,&t=x_k\\ \\ 0,&t=\text{other}\ x \end{cases}$$

这个开关函数在学术上称为「插值基函数」。只有代入 $x_k$ 时开关才打开，代入其它 $x$ 时开关关闭。

进而构造

$$f(x)=y_1L_1(x)+y_2L_2(x)+\cdots+y_nL_n(x)$$

代入 $x_k$ 时，仅 $y_k$ 后面的开关是开着的，其余皆关闭，所以函数值为 $y_k$，即 $\forall k,f(x_k)=y_k$。这就是我们想要的函数。

问题就变成了如何构造开关函数。

容易发现 $\forall k,t=x_k,$

$$\prod (t-x_i)=(t-x_1)(t-x_2)\cdots(t-x_n)=0$$

把 $(t-x_k)$ 这一项拿走，得到

$$\prod_{i\not=k}(t-x_i)=\begin{cases} \not=0,&t=x_k\\ \\ 0,&t=\text{other}\ x \end{cases}$$

为使 $t=x_k$ 时得到 $1$，再作以下改变

$$\prod_{i\not=k}\frac{t-x_i}{x_k-x_i}=\begin{cases} 1,&t=x_k\\ \\ 0,&t=\text{other}\ x \end{cases}$$

于是我们成功地构造出了开关函数

$$L_k(t)=\prod_{i\not=k}\frac{t-x_i}{x_k-x_i}$$

模板 #

#include <iostream>
#include <vector>

using namespace std;

double f(vector<double>& x, vector<double>& y, double t) {
    double result = 0.0;
    int n = x.size();
    
    for (int i = 0; i < n; i ++) {
        double term = y[i];
        for (int j = 0; j < n; j ++) {
            if (j != i) {
                term *= (t - x[j]) / (x[i] - x[j]);
            }
        }
        result += term;
    }
    
    return result;
}

int main() {
    vector<double> x = {1, 2, 3, 4, 5};
    vector<double> y = {1, 4, 9, 16, 25};

    double t;
    while (cin >> t) {
        cout << f(x, y, t) << endl;
    }

    return 0;
}