插值的故事

插值有个英文名叫interpolation，在这个点上，我觉得我们汉字表达的更实在！关于插值的起源，留给历史学家去考虑吧。我只想举一个简单的小例子。
在你的脑海中想象今天上午成都的问题是10度，下午是20度，那中午是多少度呢？你脑子里很快说出可能是15度。15度隐含了线性变化的假设。如果我们用向上凸的抛物线近似，那就过有可能就变16度。 既然是插值，我们就有各种各样的方法。

今天，我要讲的就是插值。 什么插值呢？我们假设一个一维函数$f$被$N + 1$点均分成$N$等分，函数在这$N + 1$的值都已知。我们想要通过计算程序知道函数在其它任何位置的值。对于有N+1个值的函数，我们希望用一个N阶多项式来近似。

幸运的是，对这个问题，数学家已经证明存在唯一个多项式，使得多项式在插值点的值和函数值相等。并且插值函数的误差可以有个估值范围。

我们有了唯一的多项式，现在就是构造多项式的问题。选单项式作为基函数简单，但是这个线性系统是满阵，解起来很不方便。有没有更好的办法呢？我们接下来就要看大名鼎鼎的牛顿插值法。

牛顿插值法中的$P(x)$多项式是一个$N$阶的多项式，其多项式被人为的构造成形式如下：
\begin{eqnarray}\label{1}
P(x) &=& \gamma_0 + \gamma_1 (x - x_0) + \gamma_2 (x - x_0)(x - x_1) +...+  \gamma_n (x - x_0)(x - x_1)...(x - x_{n-1}).
\end{eqnarray}

我们可以观察一下，当 $x = x_0$时，其它含有$ (x - x_0)$项的多项式为$0$，多项式$P(x)$的值变为 $P(x_0)=  \gamma_0$,那么我们就能够求得$\gamma_0$.
同样，当$x = x_1$时，基函数$x-x_0$后面的所有项都变为零，我们能够求得
\begin{eqnarray}\label{1}
P(x_1) &=& \gamma_0 + \gamma_1 (x_1 - x_0)
\end{eqnarray}
把$\gamma_0$移到左边,并用$P(x_0)$代替，加以整理可得
\begin{eqnarray}\label{2}
\gamma_1 &=& \frac{P(x_1) - P(x_0)} {x_1 - x_0}.
\end{eqnarray}
上面的表达非常有趣，这不就是普通的直线的斜率的表示方式嘛！

然后依次类推，就可以得到插值多项式所有基函数上的权值
\begin{eqnarray}\label{2}
\gamma_k &=& \frac{P(x_k) - P_{k-1}(x_k)} {\prod_{i=0}^{i=k-1} (x_k - x_i)}.
\end{eqnarray}

科学计算的其中一个关键点就出来了，效率。计算机是很傻的机器，要想让它好好的为你工作，需要在解析上做一些工作，这样能够提高效率。怎么提高呢？可能要对原线性系统做一定的变化？ 这种简单的操作就是把三角矩阵化为对角线为一的矩阵。

你需要记住一个，差商是一个数值，经常用$f[x_0,x_1]$来表示. 最后权值就是N+1个差商值。

求得了差商，那么多项式就确立了。从而就可以根据这个已知的多项式来求任何一个点的值。 写的比较简约的插值多项式是
\begin{eqnarray}\label{2}
P(x) &=& \sum_{j=0}^n f[x_0,x_1,...,x_j] \prod_{i=0}^{j-1} (x - x_i).
\end{eqnarray}
有趣的是，含有所有点信息的差商$f[x_0,x_1,...,x_j]$可以用递推关系导出。我们来看看如何得出差商的迭代公式
\begin{eqnarray}\label{2}
f[x_0,x_1,...,x_n]  &=& f[x_1,...,x_n, x_0] \\
                    &=& \frac{f[x_1,...,x_{n-1}, x_0] - f[x_1,...,x_n]} {x_0 - x_n} \\
                    &=& \frac{f[x_1,...,x_n] - f[x_1,...,x_{n-1}, x_0] } {x_n - x_0} \\
                    &=& \frac{f[x_1,...,x_n] - f[x_0, x_1,...,x_{n-1}]} {x_n - x_0} .
\end{eqnarray}

有趣的是，还有n+1点信息的$f[x_0,x_1,...,x_n]$的差商通过递归转化为更小问题的差商还有n个点信息的$f[x_1,...,x_n]$ 和n点的$f[x_0, x_1,...,x_{n-1}]$，然后继续调用递归公司求值。

递归是一种思想，它把一个问题的计算转换为结构相同的，但是规模更小的计算。实际上，这种结构加上重复性也让我们看到了事物的本质！我们来看看如何求这个系数。

我们先计算第一层，把所有N+1个的插值点的值赋给差商.
\begin{eqnarray}\label{2}
for~ (j = 0; j<=n; j++) \\
    \gamma(j) = f(j);
\end{eqnarray}
下面开始一层一层的计算，倒着从后面开始
\begin{eqnarray}\label{2}
for~ (j = 1; j<=n; j++)  \\
    for~ (i = n; i >=j; i--) \\
    \gamma(i) &=&  (\gamma(i) -  \gamma(i-1))/(x(i) - x(i-j)) \\
  \end{eqnarray}
关键是脑子要清楚，整个算法的移动，一步都少不得！如果我们展开这个算法，空间延伸结构竟然类似于二叉树形结构！！！
    
求得了系数，那多项式也可以用迭代法构造出啦.
    \begin{eqnarray}\label{2}
P_k(x) = P_{k-1}+ \gamma_k (x-x_k)
\end{eqnarray}

接下来就是程序实现了。