插值相关总结

多项式插值

$$\det (A) = \left| {\begin{array}{*{20}{c}}
1&{{x_0}}&{x_0^2}& \cdots &{x_0^n}&{}\\
1&{{x_1}}&{x_1^2}& \cdots &{x_1^n}&{}\\
{}& \cdots & \cdots & \cdots & \cdots & \cdot \\
1&{{x_n}}&{x_n^2}& \cdots &{x_n^n}&{}
\end{array}} \right|{\rm{ }}! = {\rm{ 0}}$$

$${l_i}(x) = \frac{{\left( {x - {x_0}} \right) \cdots \left( {x - {x_{i - 1}}} \right)\left( {x - {x_{i + 1}}} \right) \cdots \left( {x - {x_n}} \right)}}{{\left( {{x_i} - {x_0}} \right) \cdots \left( {{x_i} - {x_{i - 1}}} \right)\left( {{x_i} - {x_{i + 1}}} \right) \cdots \left( {{x_i} - {x_n}} \right)}} = \prod\limits_{\begin{array}{*{20}{c}}
{j = 0}\\
{j \ne i}
\end{array}}^n {\frac{{x - {x_j}}}{{{x_i} - {x_j}}}} ,\quad (i = 0,1, \cdots ,n)$$

　　注：$\frac{{x - {x_j}}}{{{x_i} - {x_j}}}$就是求解$p(x)$的$n+1$个基函数。如果对推导过程感兴趣，参考链接里有，很简单的。

function y=lagrange(x0,y0,x)
%n个数据以数组x0,y0输入，m个插值点以数组x输入，输出数组y为m个插值。
n=length(x0);m=length(x);
for i=m;
    z=x(i);
    s=0.0;
    for k=1:n
        p=1.0;
        for j=1:n;
            if j~=k
            p=p*(z-x0(j))/(x0(k)-x0(j));
            end
        end
        s=p*y0(k)+s;
    end
    y(i)=s;
end

牛顿插值

埃尔米特插值

　　不少实际的插值问题不但要求在节点上的函数值相等，而且还要求对应的导数值也相等，甚至要求高阶导数也相等，满足这种要求的插值多项式就是埃尔米特插值多项式。

　　这意味着，拟合曲线与实际曲线不仅都过点$\left(x_{i}, y_{i}\right),\left(x_{i+1}, y_{i+1}\right)$，而且两点处还有相同的切线。

　　设在节点$a \leq x_{0}<x_{1}<\cdots<x_{n} \leq b$上，$y_{j}=f\left(x_{j}\right)$，$m_{j}=f^{\prime}\left(x_{j}\right) \quad(j=0,1, \cdots, n)$

　　即满足条件：$H\left(x_{j}\right)=y_{j}, \quad H^{\prime}\left(x_{j}\right)=m_{j}, \quad(j=0,1, \cdots, n)$

　　这里共有$2n+2$个插值条件，可唯一确定一个次数不超过$2n+1$的多项式$H_{2 n+1}(x)=H(x)$，其形式为：

$$H_{2 n+1}(x)=a_{0}+a_{1} x+\cdots+a_{2 n+1} x^{2 n+1}$$

　　采用拉格朗日插值多项式的基函数方法的思想：

　　先求出$2n+2$个插值基函数，每个基函数都是$2n+1$次多项式，且满足条件：

$$\begin{array}{l}{\alpha_{j}\left(x_{k}\right)=\delta_{j k}=\left\{\begin{array}{ll}{0,} & {j \neq k} \\ {1,} & {j=k}\end{array}, \quad \alpha_{j}^{\prime}\left(x_{k}\right)=0\right.} \\ {\beta_{j}\left(x_{k}\right)=0, \quad \beta_{j}^{\prime}\left(x_{k}\right)=\delta_{j k} \quad(j, k=0,1, \cdots, n)}\end{array}$$

　　用基函数表示插值多项式：

$$H_{2 n+1}(x)=\sum_{j=0}^{n}\left[y_{j} \alpha_{j}(x)+m_{j} \beta_{j}(x)\right]$$

　　由插值基函数所满足的条件，有$H_{2 n+1}\left(x_{k}\right)=y_{k}, \quad H_{2 n+1}^{\prime}\left(x_{k}\right)=m_{k}, \quad(k=0,1, \cdots, n)$

　　然后求出基函数$\alpha_{j}(x)$与$\beta_{j}(x)$

　　利用拉格朗日插值基函数$l_{j}(x)$，

　　令$\alpha_{j}(x)=(a x+b) l_{j}^{2}(x)$

　　得到

$$\begin{array}{l}{\alpha_{j}\left(x_{j}\right)=\left(a x_{j}+b\right) l_{j}^{2}\left(x_{j}\right)=1} \\ {\alpha_{j}^{\prime}\left(x_{j}\right)=l_{j}\left(x_{j}\right)\left[a l_{j}\left(x_{j}\right)+2\left(a x_{j}+b\right) l_{j}^{\prime}\left(x_{j}\right)\right]=0}\end{array}$$

　　整理得

　　$$\left\{\begin{array}{l}{a x_{j}+b=1} \\ {a+2 l_{j}^{\prime}\left(x_{j}\right)=0}\end{array}\right.$$

　　解出

$$a=-2 l_{j}^{\prime}\left(x_{j}\right), \quad b=1+2 x_{j} l_{j}^{\prime}\left(x_{j}\right)$$

　　由于

$$l_{j}(x)=\frac{\left(x-x_{0}\right) \cdots\left(x-x_{j-1}\right)\left(x-x_{j+1}\right) \cdots\left(x-x_{n}\right)}{\left(x_{j}-x_{0}\right) \cdots\left(x_{j}-x_{j-1}\right)\left(x_{j}-x_{j+1}\right) \cdots\left(x_{j}-x_{n}\right)}$$

　　两边取对数再求导，得

$$l_{j}^{\prime}\left(x_{j}\right)=\sum_{k=0 \atop k \neq j}^{n} \frac{1}{x_{j}-x_{k}}$$

　　于是，

$$\alpha_{j}(x)=\left(1-2\left(x-x_{j}\right) \sum_{k=0 \atop k \neq j}^{n} \frac{1}{x_{j}-x_{k}}\right) l_{j}^{2}(x)$$

　　同理，得到

$$\beta_{j}(x)=\left(x-x_{j}\right) l_{j}^{2}(x)$$

　　综上所述，埃尔米特插值多项式为：

$$\left\{ {\begin{array}{*{20}{l}}
{{H_{2n + 1}}(x) = \sum\limits_{j = 0}^n {\left[ {{y_j}{\alpha _j}(x) + {m_j}{\beta _j}(x)} \right]} }\\
{{y_j} = f\left( {{x_j}} \right),{m_j} = {f^\prime }\left( {{x_j}} \right)}\\
{{\alpha _j}(x) = \left( {1 - 2\left( {x - {x_j}} \right)\sum\limits_{k = 0}^n {\frac{1}{{{x_j} - {x_k}}}} } \right)l_j^2(x),k \ne j}\\
{{\beta _j}(x) = \left( {x - {x_j}} \right)l_j^2(x)}
\end{array}} \right..$$

　　仿照拉格朗日插值余项的证明方法，得到埃尔米特插值余项：

$$R(x)=f(x)-H_{2 n+1}(x)=\frac{f^{(2 n+2)}(\xi)}{(2 n+2) !} \omega_{n+1}^{2}(x)$$

　　$n$取1时，是一个非常重要的特例，即两点三次埃尔米特插值，比较典型和常用(4个基函数均为三次多项式，$n=1$)

　　插值基函数满足条件：

$$\begin{array}{l}{\alpha_{k}\left(x_{k}\right)=1, \quad \alpha_{k}\left(x_{k+1}\right)=0, \quad \alpha_{k}^{\prime}\left(x_{k}\right)=\alpha_{k}^{\prime}\left(x_{k+1}\right)=0} \\ {\alpha_{k+1}\left(x_{k}\right)=0, \quad \alpha_{k+1}\left(x_{k+1}\right)=1, \quad \alpha_{k+1}^{\prime}\left(x_{k}\right)=\alpha_{k+1}^{\prime}\left(x_{k+1}\right)=0} \\ {\beta_{k}\left(x_{k}\right)=\beta_{k}\left(x_{k+1}\right)=0, \beta_{k}^{\prime}\left(x_{k}\right)=1, \quad \beta_{k}^{\prime}\left(x_{k+1}\right)=0} \\ {\beta_{k+1}\left(x_{k}\right)=\beta_{k+1}\left(x_{k+1}\right)=0, \beta_{k+1}\left(x_{k}\right)=0, \quad \beta_{k+1}^{\prime}\left(x_{k+1}\right)=1}\end{array}$$

　　两点三次埃尔米特插值多项式为：

$$H_{3}(x)=y_{k} \alpha_{k}(x)+y_{k+1} \alpha_{k+1}(x)+m_{k} \beta_{k}(x)+m_{k+1} \beta_{k+1}(x)$$

　　余项为：

$$R_{3}(x)=\frac{1}{4 !} f^{(4)}(\xi)\left(x-x_{k}\right)^{2}\left(x-x_{k+1}\right)^{2}, \quad \xi \in\left(x_{k}, x_{k+1}\right)$$

function y-hermite(x0,y0,y1,x);
%x0,y0为样本点数据，y1为导数指，m个插值点以数组x输入，输出数组y为m个插值
n=length(x0);
m=length(x);
for k=1:m;
    yy=0.0;
    for i=1:n
        h=1.0;
        a=0.0;
        for j=i:n
            if j~=i
                h=h*((x(k)-x0(j))/(x0(i)-x0(j)))^2; 
                a=1/(x0(i)-x0(j))+a;
            end
        end
        yy=yy+h*((x0(i)-x0(k))*(2*a*y0(i)-y(i))+y0(i));
    end
    y(k)=yy;
end

分段插值

$${l_i}(x) = \left\{ {\begin{array}{*{20}{l}}
{\frac{{x - {x_{i - 1}}}}{{{x_i} - {x_{i - 1}}}},}&{x \in \left[ {{x_{i - 1}},{x_i}} \right]{\rm{ }}i = 0}\\
{\frac{{x - {x_{i + 1}}}}{{{x_i} - {x_{i + 1}}}},}&{x \in \left[ {{x_i},{x_{i + 1}}} \right]{\rm{ }}i = n}\\
{0,}&{else}
\end{array}{\rm{ }}} \right.$$

　　method：

n=11;  m=61;
x=-5:10/(m-1):5;
y=1./(1+x.^2);
z=0*x;
x0=-5:10/(n-1):5;
y0=1./(1+x0.^2);
y1=interp1(x0, y0, x);%interp1(x0,y0,x)为Matlab中现成的分段线性插值程序
plot(x, z, ‘r‘, x, y, ‘k:‘, x, y1, ‘r‘)
gtext(‘Piece. –linear.‘), gtext(‘y=1/(1+x^2)‘)
title(‘Piecewise Linear‘)

分段三次埃尔米特插值

$${\alpha _j}(x) = \left\{ {\begin{array}{*{20}{c}}
{{{\left( {\frac{{x - {x_{j - 1}}}}{{{x_j} - {x_{j - 1}}}}} \right)}^2}\left( {1 + 2\frac{{x - {x_j}}}{{{x_{j - 1}} - {x_j}}}} \right){x_{j - 1}} \le x \le {x_j}}\\
\begin{array}{l}
{\left( {\frac{{x - {x_{j + 1}}}}{{{x_j} - {x_{j + 1}}}}} \right)^2}\left( {1 + 2\frac{{x - {x_j}}}{{{x_{j + 1}} - {x_j}}}} \right){x_j} \le x \le {x_{j + 1}}\\
{\rm{ }}0
\end{array}
\end{array}} \right.$$

$${\beta _j}(x) = \left\{ {\begin{array}{*{20}{c}}
{{{\left( {\frac{{x - {x_{j - 1}}}}{{{x_j} - {x_{j - 1}}}}} \right)}^2}\left( {x - {x_j}} \right){x_{j - 1}} \le x \le {x_j}}\\
\begin{array}{l}
{\left( {\frac{{x - {x_{j + 1}}}}{{{x_j} - {x_{j + 1}}}}} \right)^2}\left( {x - {x_j}} \right){x_j} \le x \le {x_{j + 1}}\\
{\rm{ }}0
\end{array}
\end{array}} \right.$$

x = -3:3; 
y = [-1 -1 -1 0 1 1 1]; 
xq1 = -3:.01:3;
p = pchip(x,y,xq1); %分段三次埃尔米特插值
s = spline(x,y,xq1); %逐段3次样条插值
plot(x,y,‘o‘,xq1,p,‘-‘,xq1,s,‘-.‘) 
legend(‘Sample Points‘,‘pchip‘,‘spline‘,‘Location‘,‘SouthEast‘)

三次样条插值

1、y=interp1(${x_0}$,${y_0}$,$x$,‘spline‘);
% spline改成linear，则变成线性插值
2、y=spline(${x_0}$,${y_0}$,$x$);
% 这个是根据己知的x，y数据，用样条函数插值出${x_i}$处的值。
% 由己知的x，y数据，求出它的样条函数表达式，不过该表达式不是用矩阵直接表示，要求点x`的值，要用函数
3、pp=csape(x,y,conds); y=ppval(pp,$x$);
% 各种边界条件的三次样条插值
conds是指边界条件，边界类型可为:
①‘complete‘： 给定边界一阶导数，默认的边界条件
②‘not-a-knot‘：非扭结条件，不用给边界值.
③‘periodic‘： 周期性边界条件，不用给边界值.
④‘second‘： 给定边界二阶导数.
⑤‘variational‘： 自然样条(边界二阶导数为[0,0]。