Matlab 数值计算之函数插值

微积分中对于连续或可导函数的研究虽有各路方法各显神通，但对于一个实际问题，若仅给出离散的数据点，通常首先要做的工作是根据这些离散点拟合出便于研究的连续函数。这类问题通常称作「数据拟合」，而其中的特殊情形便是函数插值。

多项式拟合插值中最「直接」的算法，就是对于 $n+1$ 个给定的数据点对 $(x_i,y_i),i=0,1,\dots ,n$，假设 $n$ 次多项式 $p_n\left(x\right)=\sum _{i=0}^n{a}_i{x}^i$ 过所有的数据点，即 $p_n\left(x_i\right)=y_i$，即：

$$\underset{V}{\underset{}{\left(\begin{array}{cccccc}1& x_0^1& x_0^2& \cdots & x_0^{n-1}& x_0^n\\ 1& x_1^1& x_1^2& \cdots & x_1^{n-1}& x_1^n\\ 1& x_2^1& x_2^2& \cdots & x_2^{n-1}& x_2^n\\ ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\\ 1& x_{n-1}^1& x_{n-1}^2& \cdots & x_{n-1}^{n-1}& x_{n-1}^n\\ 1& x_n^1& x_n^2& \cdots & x_n^{n-1}& x_n^n\end{array}\right)}}\underset{a}{\underset{}{\left(\begin{array}{c}{a}_0{x}_0^{n-1}\\ a_1{x}_1^{n-1}\\ a_2{x}_2^{n-1}\\ ⋮\ddots ⋮\\ a_{n-1}{x}_{n-1}^{n-1}\\ a_n{x}_n^{n-1}\end{array}\right)}}=\underset{y}{\underset{}{\left(\begin{array}{c}{y}_0{x}_0^{n-1}\\ y_1{x}_1^{n-1}\\ y_2{x}_2^{n-1}\\ ⋮\ddots ⋮\\ y_{n-1}{x}_{n-1}^{n-1}\\ y_n{x}_n^{n-1}\end{array}\right)}}$$

显然 $V$ 为 Vandermonde 矩阵。关于 Vandermonde 行列式的计算，这里列出两种计算方法。

$n$ 阶 Vandermonde 行列式 $$V_n=\left|\begin{array}{cccccc}1& x_1& x_1^2& \cdots & x_1^{n-2}& x_1^{n-1}\\ 1& x_2& x_2^2& \cdots & x_2^{n-2}& x_2^{n-1}\\ ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\\ 1& x_n& x_n^2& \cdots & x_n^{n-2}& x_n^{n-1}\end{array}\right|$$ 的值为 $$V_n=\prod _{1⩽i<j⩽n}\left(x_j-x_i\right)$$

以下证明过程摘自 Vandermonde Determinant。

方法一

从 $V_n$ 第二行开始，以后每行减去第一行，由初等变换得到的行列式值不变： $$\begin{array}{cc}{V}_n& =\left|\begin{array}{cccccc}1& x_1& x_1^2& \cdots & x_1^{n-2}& x_1^{n-1}\\ 1& x_2& x_2^2& \cdots & x_2^{n-2}& x_2^{n-1}\\ 1& x_3& x_3^2& \cdots & x_3^{n-2}& x_3^{n-1}\\ ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\\ 1& x_{n-1}& x_{n-1}^2& \cdots & x_{n-1}^{n-2}& x_{n-1}^{n-1}\\ 1& x_n& x_n^2& \cdots & x_n^{n-2}& x_n^{n-1}\end{array}\right|\\ & =\left|\begin{array}{cccccc}1& x_1& x_1^2& \cdots & x_1^{n-2}& x_1^{n-1}\\ 0& x_2-x_1& x_2^2-x_1^2& \cdots & x_2^{n-2}-x_1^{n-2}& x_2^{n-1}-x_1^{n-1}\\ 0& x_3-x_1& x_3^2-x_1^2& \cdots & x_3^{n-2}-x_1^{n-2}& x_3^{n-1}-x_1^{n-1}\\ ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\\ 0& x_{n-1}-x_1& x_{n-1}^2-x_1^2& \cdots & x_{n-1}^{n-2}-x_1^{n-2}& x_{n-1}^{n-1}-x_1^{n-1}\\ 0& x_n-x_1& x_n^2-x_1^2& \cdots & x_n^{n-2}-x_1^{n-2}& x_n^{n-1}-x_1^{n-1}\end{array}\right|\\ & =\left|\begin{array}{cccccc}1& 0& 0& \cdots & 0& 0\\ 0& x_2-x_1& \left(x_2-x_1\right)x_2& \cdots & \left(x_2-x_1\right)x_2^{n-3}& \left(x_2-x_1\right)x_2^{n-2}\\ 0& x_3-x_1& \left(x_3-x_1\right)x_3& \cdots & \left(x_3-x_1\right)x_3^{n-3}& \left(x_3-x_1\right)x_3^{n-2}\\ ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\\ 0& x_{n-1}-x_1& \left(x_{n-1}-x_1\right)x_{n-1}& \cdots & \left(x_{n-1}-x_1\right)x_{n-1}^{n-3}& \left(x_{n-1}-x_1\right)x_{n-1}^{n-2}\\ 0& x_n-x_1& \left(x_n-x_1\right)x_n& \cdots & \left(x_n-x_1\right)x_n^{n-3}& \left(x_n-x_1\right)x_n^{n-2}\end{array}\right|\\ & =\prod _{k=2}^n\left(x_k-x_1\right)\left|\begin{array}{cccccc}1& 0& 0& \cdots & 0& 0\\ 0& 1& x_2& \cdots & x_2^{n-3}& x_2^{n-2}\\ 0& 1& x_3& \cdots & x_3^{n-3}& x_3^{n-2}\\ ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\\ 0& 1& x_{n-1}& \cdots & x_{n-1}^{n-3}& x_{n-1}^{n-2}\\ 0& 1& x_n& \cdots & x_n^{n-3}& x_n^{n-2}\end{array}\right|\\ & =\prod _{k=2}^n\left(x_k-x_1\right)V_{n-1}\end{array}$$ 以此类推， $$V_n=\prod _{k=2}^n\left(x_k-x_1\right)\prod _{k=3}^n\left(x_k-x_2\right)V_{n-2}=\cdots =\prod _{1⩽i<j⩽n}\left(x_j-x_i\right)=\prod _{1⩽i<j⩽n}\left(x_j-x_i\right)$$

---

方法二

将行列式中 $V_n$ 中的所有 $x_n$ 替换为未知量 $x$，则 $V_n$ 可视为关于 $x$ 的 $n-1$ 次多项式，记为 $P\left(x\right)$，并按最后一行展开： $$\begin{array}{cc}P\left(x\right)& =\left|\begin{array}{cccccc}1& x_1& x_1^2& \cdots & x_1^{n-2}& x_1^{n-1}\\ 1& x_2& x_2^2& \cdots & x_2^{n-2}& x_2^{n-1}\\ 1& x_3& x_3^2& \cdots & x_3^{n-2}& x_3^{n-1}\\ ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\\ 1& x_{n-1}& x_{n-1}^2& \cdots & x_{n-1}^{n-2}& x_{n-1}^{n-1}\\ 1& x& x^2& \cdots & x^{n-2}& x^{n-1}\end{array}\right|\\ & =\left|\begin{array}{ccccc}{x}_1& x_1^2& \cdots & x_1^{n-2}& x_1^{n-1}\\ x_2& x_2^2& \cdots & x_2^{n-2}& x_2^{n-1}\\ x_3& x_3^2& \cdots & x_3^{n-2}& x_3^{n-1}\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ x_{n-1}& x_{n-1}^2& \cdots & x_{n-1}^{n-2}& x_{n-1}^{n-1}\end{array}\right|+\left|\begin{array}{ccccc}1& x_1^2& \cdots & x_1^{n-2}& x_1^{n-1}\\ 1& x_2^2& \cdots & x_2^{n-2}& x_2^{n-1}\\ 1& x_3^2& \cdots & x_3^{n-2}& x_3^{n-1}\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 1& x_{n-1}^2& \cdots & x_{n-1}^{n-2}& x_{n-1}^{n-1}\end{array}\right|x\\ & =+\cdots +\left|\begin{array}{ccccc}1& x_1& x_1^2& \cdots & x_1^{n-2}\\ 1& x_2& x_2^2& \cdots & x_2^{n-2}\\ 1& x_3& x_3^2& \cdots & x_3^{n-2}\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 1& x_{n-1}& x_{n-1}^2& \cdots & x_{n-1}^{n-2}\end{array}\right|x^{n-1}\end{array}$$ 显然 $P\left(x_1\right)=P\left(x_2\right)=\cdots =P\left(x_{n-1}\right)=0$，由多项式分解定理可知 $$P\left(x\right)=V_{n-1}\prod _{i=1}^{n-1}(x-x_i)$$ 将 $x$ 替换成 $x_n$，由递推可得： $$V_n=\prod _{1⩽i<n}(x_n-x_i)V_{n-1}=\prod _{1⩽i<n}(x_n-x_i)\prod _{1⩽i<n-1}(x_{n-1}-x_i)V_{n-2}=\cdots =\prod _{1⩽i<j⩽n}(x_j-x_i)$$

---

但是解这个线性方程组是很不明智的，如果数据点较近，$V$ 会接近奇异矩阵。而拉格朗日插值法则不需要硬着头皮解方程。

拉格朗日插值

$n$ 次 Lagrange 插值多项式系数 $l_0\left(x\right),l_1\left(x\right),\dots ,l_n\left(x\right)$ 定义为 $l_i\left(x_j\right)=\{\begin{array}{cc}1,& i=j\\ 0,& i\ne j\end{array}$，则多项式可表示为

$$P\left(x\right)=y_0{l}_0\left(x\right)+y_1{l}_1\left(x\right)+\cdots +y_n{l}_n\left(x\right)=\sum _{i=0}^n{y}_i{l}_i\left(x\right)$$

注意这里的每个 $l_i\left(x\right)$ 均为 $n$ 次多项式，且由定义可知 $l_i\left(x\right)$ 有 $n$ 个根

$$x_0,x_1,x_2,\dots ,x_{i-1},x_{i+1},\dots ,x_n$$

所以 $l_i\left(x\right)$ 可表示成 $l_i\left(x\right)=C_i\prod _{j\ne i}(x-x_j)$，进一步由 $l_i\left(x\right)$ 定义，

$$1=l_i\left(x_i\right)=C_i\prod _{j\ne i}(x_i-x_j)\Rightarrow C_i=\frac{1}{\prod _{j\ne i}(x_i-x_j)}\Rightarrow l_i\left(x\right)=\frac{\prod _{j\ne i}(x-x_j)}{\prod _{j\ne i}(x_i-x_j)}$$

所以 Lagrange 插值多项式为

$$p\left(x\right)=\sum _{i=0}^n{y}_i\frac{\prod _{j\ne i}(x-x_j)}{\prod _{j\ne i}(x_i-x_j)}=\sum _{i=0}^n{y}_i\frac{(x-x_0)(x-x_1)\cdots (x-x_{i-1})(x-x_{i+1})\cdots (x-x_n)}{(x_i-x_0)(x_i-x_1)\cdots (x_i-x_{i-1})(x_i-x_{i+1})\cdots (x_i-x_n)}$$

从形式上看，Lagrange 插值法是容易实现的，但如果加入新的数据点，则需要重新计算多项式系数，即算法随着数据点的增多不能复用以前的计算结果。

function [f, f0] = interp_lagrange(x, y, x0, prec)
    if nargin <= 3, prec = 4;

    if length(x) == length(y), n = length(x) - 1;
    else error('Incompatible dimension of x and y !!');
    end

    syms t;
    f = 0;

    for k = 1:n
        xs = [x(1:k-1) x(k+1:n)];
        f = f + y(k) * prod((t-xs)./(x(k)-xs));
    end

    f = vpa(f, prec);
    if nargin > 2, f0 = vpa(subs(f, t, x0), prec); end;
end

艾特肯插值

首先，由点 $(x_0,y_0)$ 到 $(x_j,y_j),j=1,2,\dots ,n$ 可以确定 $n$ 组一次多项式，即穿过两点的直线；第二步，根据 $$，进一步生成 $n-1$ 组二次多项式。次过程共重复 $n$ 次，最终得到一个 $n$ 次多项式。

具体而言，由点 $(x_0,y_0)$ 及 $(x_1,y_1)$ 得到的一次多项式为

$$p_{01}\left(x\right)=\frac{x-x_1}{x_0-x_1}{y}_0-\frac{x-x_0}{x_1-x_0}{y}_1=\frac{1}{x_1-x_0}\left|\begin{array}{cc}{y}_0& x_0-x\\ y_1& x_1-x\end{array}\right|$$

类比之，由 $x_0$ 及 $x_j$ 确定的一次多项式表示为

$$p_{0j}\left(x\right)=\frac{1}{x_j-x_0}\left|\begin{array}{cc}{y}_0& x_0-x\\ y_j& x_j-x\end{array}\right|,j=1,2,\dots ,n$$

进一步，由点 $(x_0,y_0),(x_1,y_1),(x_j,y_j)$ 确定的二次多项式为

$$p_{01j}\left(x\right)=\frac{1}{x_j-x_1}\left|\begin{array}{cc}{p}_{01}\left(x\right)& x_1-x\\ p_{0j}\left(x\right)& x_j-x\end{array}\right|,j=2,3,\dots ,n$$

一般而言，对于 $k+2$ 个数据点 $(x_0,y_0),(x_1,y_1),\dots ,(x_k,y_k)$ 及 $(x_j,y_j)$，$k+1$ 次插值多项式为

$$p_{012\cdots kj}\left(x\right)=\frac{1}{x_j-x_k}\left|\begin{array}{cc}{p}_{012\cdots k}\left(x\right)& x_k-x\\ p_{012\cdots (k-1)}\left(x\right)& x_j-x\end{array}\right|j=k+1,\dots ,n$$

整个计算过程如下表

$$\begin{array}{ccccccc}{x}_j& y_j& p_{0j}& p_{01j}& p_{012j}& p_{0123j}& x_j-x\\ x_0& y_0& & & & & x_0-x\\ x_1& y_1& p_{01}& & & & x_1-x\\ x_2& y_2& p_{02}& p_{012}& & & x_2-x\\ x_3& y_3& p_{03}& p_{013}& p_{0123}& & x_3-x\\ x_4& y_4& p_{04}& p_{014}& p_{0124}& p_{01234}& x_4-x\\ \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots \\ x_n& y_n& p_{0n}& p_{01n}& p_{012n}& p_{0123n}& x_n-x\end{array}$$

function [f, f0] = interp_aitken(x, y, x0, prec)
    if nargin <= 3, prec = 4;
    if length(x) == length(y), n = length(x) - 1;
    else error('Incompatible dimension of x and y !!');
    end

    syms t;
    f = 0;
    lastp = y;
    xs = x - t;

    for k = 1:n
        for j = k+1:n+1
            curp(j) = (lastp(k)*xs(j) - lastp(j)*xs(k)) / (xs(j)-xs(k));
        end
        lastp = curp;
    end

    f = vpa(simplify(lastp(end)), prec);
    if nargin > 2, f0 = vpa(subs(f, t, x0), prec); end;
end