插值相关总结

　　插值的通俗解释就是一种用一些已知的数据去预测想要的数据的方法。

多项式插值

　　多项式插值是最常见的一种函数插值(插值函数为多项式)。

$${p_n}(x) = {a_0} + {a_1}x + {a_2}{x^2} +  cdots  + {a_n}{x^n}$$

　　从几何上看可以理解为：已知平面上n+1个不同点，要寻找一条n次插值多项式函数$p(x)$通过曲线$f(x)$上已知的这n+1个点。使$p(x)$接近$f(x)$。

　　而将n个点代入多项式函数，则可用方程组表示，即

$$left{egin{array}{l}{a_{0}+a_{1} x_{0}+a_{2} x_{0}^{2}+cdots+a_{n} x_{0}^{n}=y_{0}} \ {a_{0}+a_{1} x_{1}+a_{2} x_{1}^{2}+cdots+a_{n} x_{1}^{n}=y_{1}} \ {ldots ldots ldots ldots ldots ldots ldots ldots ldots ldots ldots ldots} \ {a_{0}+a_{1} x_{n}+a_{2} x_{n}^{2}+cdots+a_{n} x_{n}^{n}=y_{n}}end{array} ight.$$

　　当系数矩阵满秩时，有唯一解。系数矩阵的行列式为以范德蒙德行列式，因此只要$n+1$个点互不相同(意味着范德蒙德行列式不为0)，就一定能够解出系数方程，得到多项式函数的系数${a_0},{a_1},...,{a_n}$。

$$det (A) = left| {egin{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}} ight|{ m{ }}! = { m{ 0}}$$

 

　　插值多项式一般有两种常见的表达形式，一个是拉格朗日插值多项式，另一个是牛顿插值多项式。

拉格朗日插值

　　一般的多项式插值即求线性方程组的解即可，但是这样求解的复杂度比较高。

　　基于一般多项式插值的求解思想，拉格朗日插值法诞生了。

　　1779年英国数学家爱德华·华林于最早发现拉格朗日插值法。

　　不久后1783年，由莱昂哈德·欧拉再次发现。

　　直到1795年，拉格朗日在其著作《师范学校数学基础教程》中发表了这个插值方法。

　　为了降低解决问题的复杂度，引入了“基函数”的概念。

　　首先构造一组基函数：

$${l_i}(x) = frac{{left( {x - {x_0}} ight) cdots left( {x - {x_{i - 1}}} ight)left( {x - {x_{i + 1}}} ight) cdots left( {x - {x_n}} ight)}}{{left( {{x_i} - {x_0}} ight) cdots left( {{x_i} - {x_{i - 1}}} ight)left( {{x_i} - {x_{i + 1}}} ight) cdots left( {{x_i} - {x_n}} ight)}} = prodlimits_{egin{array}{*{20}{c}}
{j = 0}\
{j e 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$个基函数。如果对推导过程感兴趣，参考链接里有，很简单的。

　　${l_i}(x)$是$n$次多项式，满足

$${l_i}left( {{x_j}} ight) = left{ {egin{array}{*{20}{l}}
0&{j e i}\
1&{j = i}
end{array}} ight.$$

　　令

$${L_n}(x) = sumlimits_{i = 0}^n {{y_i}} {l_i}(x) = sumlimits_{i = 0}^n {{y_i}} left( {prodlimits_{egin{array}{*{20}{c}}
{j = 0}\
{j = i}
end{array}}^n {frac{{x - {x_j}}}{{{x_i} - {x_j}}}} } ight)$$

　　这就是$n$次的拉格朗日插值多项式，且$n+1$的$n$次拉格朗日插值多项式是惟一的。

　　拉格朗日插值法可以找到一个多项式，其恰好在各个观测的点取到观测到的值，这样的多项式称为拉格朗日插值多项式。

　　数学上来说，拉格朗日插值法可以给出一个恰好穿过二维平面上若干个已知点的多项式函数。

　　&&：拉格朗日插值法在插值区间内插值的精度远远大于区间外的精度，所以，在区间外，拉格朗日插值一般是不准确的。

　　MATLAB：对一维数据进行拉格朗日插值

function [ yi ] = lagrange_interp (X,Y,xi)
n=length(X);       %得到已知数据长度
m=length(xi);      %得到待插值数据长度
yi=zeros(size(xi));
for j=1:m          %待插值数据有m个，计算每个插值结果
    for i=1:n       %已知的n个数据构造中间值
temp=1;   %temp用于存储中间值
        for k=1:n
            if(i~=k)  %和自身标号相同的不相乘
                temp=temp*(xi(j)-X(k))/(X(i)-X(k));
            end
        end
        yi(j)=Y(i)*temp+yi(j);
    end
end
end

牛顿插值

　　用拉格朗日插值法，如果增加新的节点，需要从头重新计算，而牛顿差值法机智的引入了差商的概念，使其在差值节点增加时便于计算。

1. 差商

　　差商：设有函数$f(x)$，${x_0},{x_1},{x_2}, cdots $为一系列互不相同的点，称$frac{fleft(x_{i} ight)-fleft(x_{j} ight)}{x_{i}-x_{j}}(i eq j)$为$f(x)$关于点$x_{i}, x_{j}$的一阶差商，记为$fleft[ {{x_i},{x_j}} ight]$。

　　即$$fleft[x_{i}, x_{j} ight]=frac{fleft(x_{i} ight)-fleft(x_{j} ight)}{x_{i}-x_{j}}$$

　　二阶差商即一阶差商的差商，即

$$fleft[ {{x_i},{x_j},{x_k}} ight] = frac{{fleft[ {{x_i},{x_j}} ight] - fleft[ {{x_j},{x_k}} ight]}}{{{x_i} - {x_k}}}$$

　　同理，$k$阶差商为：

$$fleft[x_{0}, x_{1}, cdots, x_{k} ight]=frac{fleft[x_{0}, x_{1}, cdots, x_{k-1} ight]-fleft[x_{1}, x_{2}, cdots, x_{k} ight]}{x_{0}-x_{k}}$$

　　而且，差商中各点的顺序是可以轮换的：

$$egin{array}{l}{fleft[x_{i}, x_{j} ight]=fleft[x_{j}, x_{i} ight]} \ {fleft[x_{i}, x_{j}, x_{k} ight]=fleft[x_{i}, x_{k}, x_{j} ight]=fleft[x_{j}, x_{i}, x_{k} ight]}end{array}$$

2. 牛顿插值公式

　　①把差商定义式展开，得到差商公式表：

$$egin{array}{l}{f(x)=fleft(x_{p} ight)+fleft[x, x_{0} ight]left(x-x_{p} ight)} \ {fleft[x, x_{0} ight]=fleft[x_{0}, x_{1} ight]+fleft[x, x_{0}, x_{1} ight]left(x-x_{2} ight)} \ {fleft[x, x_{0}, x_{1} ight]=fleft{x_{0}, x_{1}, x_{2}, x_{2} ight]+fleft[x, x_{0}, x_{1}, x_{2} ight]left(x-x_{2} ight)} \ {cdots} \ {fleft[x, x_{0}, cdots, x_{n-1} ight]=fleft{x_{0}, x_{1}, cdots, x_{n} ight]+fleft[x, x_{0}, cdots, x_{n} ight]left(x-x_{n} ight)}end{array}$$

　　②自底向上归纳得到，

$$egin{array}{l}{f(x)=fleft(x_{1} ight)+fleft[x_{n}, x_{1} ight]left(x-x_{1} ight)} \ {+fleft[x_{n}, x_{1}, x_{2} ight]left(x-x_{1} ight)left(x-x_{1} ight)+cdots+fleft[x_{n}, x_{n} cdots, x_{n} ight]left(x-x_{1} ight)left(x-x_{1} ight) cdotleft(x-x_{n-1} ight)} \ {+fleft[x_{n}, x_{n}, cdots, x_{1} ight]left(x-x_{1} ight)left(x-x_{1} ight) cdotleft(x-x_{n} ight)}end{array}$$

　　③由于最后一项含有未知数$x$，因此去掉作为余项，便可得到牛顿插值公式。

$$egin{aligned} N_{n}(x) &=fleft(x_{0} ight)+left(x-x_{0} ight) fleft[x_{0}, x_{1} ight]+cdots \ &+left(x-x_{0} ight)left(x-x_{1} ight) cdotsleft(x-x_{n-1} ight) fleft[x_{0}, x_{1}, cdots, x_{n} ight] end{aligned}$$

　　其插值余项记为

$$egin{aligned} R_{n}(x) &=left(x-x_{0} ight)left(x-x_{1} ight) cdotsleft(x-x_{n} ight) fleft[x, x_{0}, x_{1}, cdots, x_{n} ight] \ &=omega_{n+1}(x) fleft[x, x_{0}, x_{1}, cdots, x_{n} ight] end{aligned}$$

　　MATLAB：对一维数据进行牛顿插值

function yi=newton_interp(X,Y,xi)
syms t;             %定义自变量t，用于字符公式
if(length(X)==length(Y))
    n=length(X);
    c(1:n)=0.0;
else
    disp('X和Y的维数不相等！');
    return;
end
f=Y(1);             %f用来记录得到的牛顿插值公式的字符串表达式
l=1;
for i=1:n-1
    y1=zeros(1,n-i);
    for j=i+1:n
        y1(j)=(Y(j)-Y(i))/(X(j)-X(i));
    end
    c(i)=y1(i+1);   %c记录差分
    l=l*(t-X(i));    %l记录(x-x0)(x-x1)……的值 
    f=f+c(i)*l;     %累加得到差分公式
    Y=y1;
end
f=simplify(f);       %简化得到的牛顿插值公式
m=length(xi);       %开始输出
for i=1:m
    yi(i)=subs(f,'t',xi(i));   % 根据公式计算需要的值
end
yi=double(yi);     % 转换为数值型，为返回值
end

　　对比拉格朗日插值法：

　　与拉格朗日插值法相比，牛顿插值法增加节点时，不需要重新计算，在某些情景下比拉格朗日插值法更好用，而且在计算余项时，牛顿插值的余项由于不需要导数，故$f(x)$是由离散点或者导数不存在时仍然适用，这是拉格朗日余项计算所不能比拟的。

　　差商与导数的关系为：

$$fleft[x_{0}, x_{1}, cdots, x_{n} ight]=frac{f^{(n)}(xi)}{n !}$$

　　其中，$xi in(alpha, eta), alpha=min _{0 leq i leq n}left{x_{i} ight}, eta=max _{0 leq i leq n}left{x_{i} ight}$

泰勒插值

　　此外，还有Taylor插值，不过，一般基本用不到。

　　泰勒多项式：

$${p_x}(x) = fleft( {{x_0}} ight) + {f^prime }left( {{x_0}} ight)left( {x - {x_0}} ight) + frac{{{f^{prime prime }}left( {{x_0}} ight)}}{{2!}}{left( {x - {x_0}} ight)^2} +  cdots  + frac{{{f^{(n)}}left( {{x_0}} ight)}}{{n!}}{left( {x - {x_0}} ight)^n}$$ 

　　而泰勒插值的条件需要满足$0-n$阶导数可导：$p_n^{(k)}left( {{x_0}} ight) = {f^{(k)}}left( {{x_0}} ight),quad k = 0,1, cdots ,n$

　　但是，满足n阶可导这个条件过于苛刻，所以泰勒插值并不常用。

埃尔米特插值

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

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

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

　　即满足条件：$Hleft(x_{j} ight)=y_{j}, quad H^{prime}left(x_{j} ight)=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$次多项式，且满足条件：

$$egin{array}{l}{alpha_{j}left(x_{k} ight)=delta_{j k}=left{egin{array}{ll}{0,} & {j eq k} \ {1,} & {j=k}end{array}, quad alpha_{j}^{prime}left(x_{k} ight)=0 ight.} \ {eta_{j}left(x_{k} ight)=0, quad eta_{j}^{prime}left(x_{k} ight)=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} eta_{j}(x) ight]$$

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

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

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

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

　　得到

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

　　整理得

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

　　解出

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

　　由于

$$l_{j}(x)=frac{left(x-x_{0} ight) cdotsleft(x-x_{j-1} ight)left(x-x_{j+1} ight) cdotsleft(x-x_{n} ight)}{left(x_{j}-x_{0} ight) cdotsleft(x_{j}-x_{j-1} ight)left(x_{j}-x_{j+1} ight) cdotsleft(x_{j}-x_{n} ight)}$$

　　两边取对数再求导，得

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

　　于是，

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

　　同理，得到

$$eta_{j}(x)=left(x-x_{j} ight) l_{j}^{2}(x)$$

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

$$left{ {egin{array}{*{20}{l}}
{{H_{2n + 1}}(x) = sumlimits_{j = 0}^n {left[ {{y_j}{alpha _j}(x) + {m_j}{eta _j}(x)} ight]} }\
{{y_j} = fleft( {{x_j}} ight),{m_j} = {f^prime }left( {{x_j}} ight)}\
{{alpha _j}(x) = left( {1 - 2left( {x - {x_j}} ight)sumlimits_{k = 0}^n {frac{1}{{{x_j} - {x_k}}}} } ight)l_j^2(x),k e j}\
{{eta _j}(x) = left( {x - {x_j}} ight)l_j^2(x)}
end{array}} ight..$$

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

$$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$)

　　插值基函数满足条件：

$$egin{array}{l}{alpha_{k}left(x_{k} ight)=1, quad alpha_{k}left(x_{k+1} ight)=0, quad alpha_{k}^{prime}left(x_{k} ight)=alpha_{k}^{prime}left(x_{k+1} ight)=0} \ {alpha_{k+1}left(x_{k} ight)=0, quad alpha_{k+1}left(x_{k+1} ight)=1, quad alpha_{k+1}^{prime}left(x_{k} ight)=alpha_{k+1}^{prime}left(x_{k+1} ight)=0} \ {eta_{k}left(x_{k} ight)=eta_{k}left(x_{k+1} ight)=0, eta_{k}^{prime}left(x_{k} ight)=1, quad eta_{k}^{prime}left(x_{k+1} ight)=0} \ {eta_{k+1}left(x_{k} ight)=eta_{k+1}left(x_{k+1} ight)=0, eta_{k+1}left(x_{k} ight)=0, quad eta_{k+1}^{prime}left(x_{k+1} ight)=1}end{array}$$

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

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

　　余项为：

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

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

分段插值

　　分段插值：对每一个分段区间$({x_i},{x_{i + 1}})$分别进行插值，则最后所得插值函数为分段函数。

　　很常用！

　　对于分段插值，我们主要探讨分段低次插值，因为高次插值存在病态性质，会出现大幅度波动现象，不能很好地拟合。

　　这就是个高次插值的图像，可以看出到了$x pm 5$附近时，拟合曲线偏离实际曲线很远。

分段线性插值

　　分段线性插值即把每两个相邻的节点用直线连起来，这样形成的折线就是分段线性插值函数，记作${I_n}(x)$，它满足$I_{n}left(x_{i} ight)=y_{i}$，且${I_n}(x)$在每个小区间$left[x_{i}, x_{i+1} ight]$上是线性函数。$(i=0,1, cdots, n)$

$$I_{n}(x)=sum_{i=0}^{n} y_{i} l_{i}(x)$$

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

　　用${I_n}(x)$计算$x$点的插值时，只用到$x$左右的两个节点，计算量与节点数$n$无关。但$n$越大，分段越多，插值误差越小。
　　MATLAB：

$$mathrm{y}= ext { interp } 1left(mathrm{x} 0, mathrm{y0}, mathrm{x}, ext { 'method }^{prime} ight)$$

　　method：

nearest	最近项插值
linear	线性插值
spline	逐段3次样条插值
cubic	保凹凸性3次插值

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')

分段三次埃尔米特插值

　　分段差值简单易行，又克服了Runge现象，但它却导致一阶导数不连续，为了克服线性插值一阶导数不连续的缺点，可以使用分段三次埃尔米特插值。

　　分段三次埃尔米特插值函数即导数连续的分段插值函数。
　　根据两点三次埃尔米特插值多项式，得到${I_h}(x)$在$left[x_{k}, x_{k+1} ight]$的表达式为：

$${I_h}(x) = {left( {frac{{x - {x_{k + 1}}}}{{{x_k} - {x_{k + 1}}}}} ight)^2}left( {1 + 2frac{{x - {x_k}}}{{{x_{k + 1}} - {x_k}}}} ight){f_k} + {left( {frac{{x - {x_k}}}{{{x_{k + 1}} - {x_k}}}} ight)^2} cdot left( {1 + 2frac{{x - {x_{k + 1}}}}{{{x_k} - {x_{k + 1}}}}} ight){f_{k + 1}} + {left( {frac{{x - {x_{k + 1}}}}{{{x_k} - {x_{k + 1}}}}} ight)^2}left( {x - {x_k}} ight){f^prime }left( {{x_k}} ight) + {left( {frac{{x - {x_k}}}{{{x_{k + 1}} - {x_k}}}} ight)^2}left( {x - {x_{k + 1}}} ight)f_{k + 1}^prime $$
　　若在整个区间$[a, b]$上定义一组分段三次插值基函数$alpha_{j}(x)$及$eta_{j}(x)$,$(j=0,1, cdots, n)$，则$I_{h}(x)$可表示为：

　　$$I_{h}(x)=sum_{j=0}^{n}left[f_{j} alpha_{j}(x)+f_{j}^{prime} eta_{j}(x) ight]$$

　　其中$alpha_{j}(x)$，$eta_{j}(x)$分别表示为

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

$${eta _j}(x) = left{ {egin{array}{*{20}{c}}
{{{left( {frac{{x - {x_{j - 1}}}}{{{x_j} - {x_{j - 1}}}}} ight)}^2}left( {x - {x_j}} ight){x_{j - 1}} le x le {x_j}}\
egin{array}{l}
{left( {frac{{x - {x_{j + 1}}}}{{{x_j} - {x_{j + 1}}}}} ight)^2}left( {x - {x_j}} ight){x_j} le x le {x_{j + 1}}\
{ m{ }}0
end{array}
end{array}} ight.$$

　　$I_{h}(x)$可表示为：

$$egin{array}{c}{I_{h}(x)=f_{k} alpha_{k}(x)+f_{k+1} alpha_{k+1}(x)+f_{k}^{prime} eta_{k}(x)+f_{k+1}^{prime} eta_{k+1}(x)} \ {left(x_{k} leq x leq x_{k+1} ight)}end{array}$$
　　MATLAB：来自官方文档 https://ww2.mathworks.cn/help/matlab/ref/pchip.html

　　将 spline 和 pchip 为两个不同函数生成的插值结果进行比较。

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')

　　为了避免高次插值可能出现的大幅度波动现象，在实际应用中通常采用分段低次插值来提高近似程度，比如可用分段线性插值或分段三次埃尔米特插值来逼近已知函数，但它们的总体光滑性较差。为了克服这一缺点，一种全局化的分段插值方法，即三次样条插值，下面我们单独说一下。

三次样条插值

　　分段插值的优点是误差小，稳定性高，但是曲线的光滑性不好，而许多实际问题需要插值曲线有较高阶的整体光滑性，理论上是可以使用高次埃尔米特插值或分段高次埃尔米特插值，但是必须知道每一个点的导数是不现实的。所以，样条插值诞生了。理论推导参考：https://www.cnblogs.com/duye/p/8671820.html

　　MATLAB：

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]。

　　

clear
clc 
x0=[0,3,5,7,9,11,12,13,14,15]; 
y0=[0,1.2,1.7,2.0,2.1,2.0,1.8,1.2,1.0,1.6]; 
t=0:0.05:15;
y1=spline(x0,y0,t); 
dy1=(spline(x0,y0,0.0001)-spline(x0,y0,0))/0.0001%x=0处斜率 
min1=min(spline(x0,y0,13:0.001:15))%13到15最小值 
figure
plot(x0,y0,'ro',t,y1);%画出曲线 
title('三次样条插值'); 

二维插值

　　二维插值是指被插值函数$z = fleft( {x,y} ight)$为二元函数，画出的图像是三维的。

网络节点插值

　　z=interp2(x0,y0,z0,x,y,'method')

　　其中，

x0，y0是节点坐标(分别为行向量和列向量)；
z0是节点的值；
method和前面所述相同。
$x$，$y$为插值点坐标，$z$为函数值。

如果是三次样条插值

pp=csape({x0,y0},z0,conds)

z=fnval(pp,{x,y})

“conds”与一维相同，一般默认。

clear
clc
x=100:100:400;
y=100:100:400;
z=[636,697,624,478;
    712,630,478,420;
    674,598,412,400;
    626,552,334,310];
pp=csape({x,y},z');  %z矩阵的行列对应向量x,y
xi=100:10:500;
yi=100:10:400;
cz1=fnval(pp,{xi,yi});
cz2=interp2(x,y,z,xi,yi','spline');
[i,j]=find(cz1==max(max(cz1)))
subplot(1,2,1);
surf(xi,yi,cz1');
shading interp;   %插入颜色插值
axis equal;
title('cz1');
subplot(1,2,2);
surf(xi,yi,cz2);
shading interp;
axis equal;
title('cz2');

散乱节点插值

 　　zi=griddata(x,y,z,xi,yi)

　　x，y，z为数据点，xi，yi为插值点，返回zi的值。

clear
clc
%样本点信息
x=[129,140,103.5,88,185.5,95,105,157.5,107.5,77,81,162,162,117.5];
y=[7.5,141.5,23,147,22.5,137.5,85.5,-6.5,-81,3,56.5,-66.5,84,-33.5];
z=-[4,8,6,8,6,8,8,9,9,8,8,9,4,9];
xi=75:200;
yi=-50:150;
zi=griddata(x,y,z,xi,yi','cubic');
subplot(1,2,1);
plot(x,y,'*');
title('xy');
subplot(1,2,2);
mesh(xi,yi,zi);
shading interp;
axis equal;
title('xyz');

　　

 

 

 

 

部分内容的参考链接：

[1]:https://www.cnblogs.com/BlueMountain-HaggenDazs/p/4271977.html

[2]:https://www.cnblogs.com/duye/p/8671820.html

[3]:https://wenku.baidu.com/view/53498f69492fb4daa58da0116c175f0e7cd1190c.html