分段线性插值

应用的函数为：
y=interp1(x0,y0,x)或y=interp1(x0,y0,x,’linear’)
其中的参数表示为：

• x0，y0表示的初始的插值节点向量
  -** x表示要得到的插值节点对应的横坐标向量**
  -** y表示的是返回值，返回的是要求得的插值节点的纵坐标**
• 注意函数当中的最后一个字符是1不是l

三次样条插值

三次样条插值是要保证插值函数要在插值节点上的导数相同
应用的函数为：
**y=interp1(x0,y0,x,’spline’)**或 y=spline(x0,y0,x)
其中的参数表示为：

• x0，y0表示的初始的插值节点向量
  -** x表示要得到的插值节点对应的横坐标向量**
  -** y表示的是返回值，返回的是要求得的插值节点的纵坐标**

例子

对 y=1/(1+x^2)在[-5, 5]上, 用n=11个等距分点作分段线性插值和三次样条插值, 用m=21个插值点作图,比较结果

n=11, m=21;
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);
y2=interp1(x0,y0,x,'spline');
[x' y' y1' y2']
plot(x,z,'r',x,y,'k:',x,y1,'b',x,y2,'g');
gtext('Piece.-linear.'),gtext('Spline'),gtext('y=1/(1+x^2)');

结果为

首先：
几种常用的插值方法主要有：Lagrange插值、Newton插值、Hermite插值；分段插值方法主要有：分段线性插值、分段三次Hermite插值、三次样条插值。
接下来：
已知的插值公式：

已知的分段插值公式：

将以上插值公式使用Matlab算法实现：
X为x的行向量，Y为y的行向量，Z为y的一阶导数向量或者边界一阶导数行向量。不同的方法用到的部分代码可能相同。
代码：

Lagrange插值公式求解及作图代码：
function [L]=Lagrange(X,Y)
[C,L,L1,l]=Lagran(X,Y);
X1=1:0.01:5;
Y1=polyval(C,X1);
plot(X,Y,'*',X1,Y1,'r');
grid on;
xlabel('x');
ylabel('y');

function [C,L,L1,l]=Lagran(X,Y)
m=length(X);
L=ones(m,m);
for k=1:m
    v=1;
    for i=1:m
        if k~=i
            v=conv(v,poly(X(i)))/(X(k)-X(i));
        end
    end
    L1(k,:)=v;
    l(k,:)=poly2sym(v);
end
C=Y*L1;
L=Y*l;

Newton插值公式求解及作图代码：
function [L]=Newton(X,Y)
[A,C,L]=Newploy(X,Y);
X1=1:0.01:5;
Y1=polyval(C,X1);
plot(X,Y,'*',X1,Y1,'r');
grid on
xlabel('x');
ylabel('y');

function [A,C,L]=Newploy(X,Y)
n=length(X);
A=zeros(n,n);
A(:,1)=Y';
s=0.0;
p=1.0;
q=1.0;
c1=1.0;
for j=2:n
    for i=j:n
        A(i,j)=(A(i,j-1)-A(i-1,j-1))/(X(i)-X(i-j+1));
    end
    b=poly(X(j-1));
    q1=conv(q,b);
    c1=c1*j;
    q=q1;
end
C=A(n,n);
b=poly(X(n));
q1=conv(q1,b);
for k=(n-1):-1:1
    C=conv(C,poly(X(k)));
    d=length(C);
    C(d)=C(d)+A(k,k);
end
L(k,:)=poly2sym(C);

Hermite插值公式求解及作图代码：
function [L]=Hermite(X,Y,Z)
[A,C,L]=Herm1(X,Y,Z);
X1=1:0.01:5;
Y1=polyval(C,X1);
plot(X,Y,'*',X1,Y1,'r');
grid on
xlabel('x');
ylabel('y');

function [A,C,L]=Herm1(X,Y,Z)
s=length(X);
n=2*s;
X2=zeros(1,n);
Y2=zeros(1,n);
Z2=zeros(1,n);
A=zeros(n,n);
for m=1:s
    X2(2*m-1)=X(m);
    X2(2*m)=X(m);
    Y2(2*m-1)=Y(m);
    Y2(2*m)=Y(m);
    Z2(2*m)=Z(m);
end
A(:,1)=Y2';
s=0.0;
p=1.0;
q=1.0;
c1=1.0;
for j=2:n
    for i=j:n
        if X2(i)==X2(i-j+1)
            A(i,j)=Z2(i);
        else
        A(i,j)=(A(i,j-1)-A(i-1,j-1))/(X2(i)-X2(i-j+1));
        end
    end
    b=poly(X2(j-1));
    q1=conv(q,b);
    c1=c1*j;
    q=q1;
end
C=A(n,n);
b=poly(X2(n));
q1=conv(q1,b);
for k=(n-1):-1:1
    C=conv(C,poly(X2(k)));
    d=length(C);
    C(d)=C(d)+A(k,k);
end
L(k,:)=poly2sym(C);

分段线性插值公式求解及作图代码：
function FD(X,Y)
n=length(X);
X1=X;
Y1=Y;
for i=1:n-1
    x=X1(i):0.01:X1(i+1);
    y=Y1(i)*(x-X1(i+1))/(X1(i)-X1(i+1))+Y1(i+1)*(x-X1(i))/(X1(i+1)-X1(i));
    plot(X,Y,'*',X1,Y1,'r');
    hold on
end
xlabel('x');
ylabel('y');
syms x;
for i=1:n-1
    y=Y1(i)*(x-X1(i+1))/(X1(i)-X1(i+1))+Y1(i+1)*(x-X1(i))/(X1(i+1)-X1(i))
end

分段三次Hermite插值公式求解及作图代码：
function TH(X,Y,Z)
n=length(X);
X3=zeros(2);
Y3=zeros(2);
Z3=zeros(2);
for i=1:n-1
    X3(1)=X(i);
    X3(2)=X(i+1);
    Y3(1)=Y(i);
    Y3(2)=Y(i+1);
    Z3(1)=Z(i);
    Z3(2)=Z(i+1);
    [A,C,L]=Herm1(X3,Y3,Z3)
    X1=X3(1):0.01:X3(2);
    Y1=polyval(C,X1);
    plot(X,Y,'*',X1,Y1,'r');
    grid on
    hold on
end
xlabel('x');
ylabel('y');

三次样条插值公式求解及作图代码：
function SY(X,Y,Z)
n=length(X);
Y1=zeros(n,1);
A=2*eye(n);
for i=2:n-1
    A(i,i-1)=(X(i)-X(i-1))/(X(i+1)-X(i-1));
    A(i,i+1)=1-(X(i)-X(i-1))/(X(i+1)-X(i-1));
end
A(1,2)=1;
A(n,n-1)=1;
Y1(1)=6/(X(2)-X(1))*((Y(2)-Y(1))/(X(2)-X(1))-Z(1));
Y1(n)=6/(X(n)-X(n-1))*(Z(2)-(Y(n)-Y(n-1))/(X(n)-X(n-1)));
for i=2:n-1
    Y1(i)=6/(X(i+1)-X(i-1))*((Y(i+1)-Y(i))/(X(i+1)-X(i))-(Y(i)-Y(i-1))/(X(i)-X(i-1)));
end
U=A;
L=eye(n);
for k=1:n-1
    L(k+1:n,k)=U(k+1:n,k)/U(k,k);
    U(k+1:n,k+1:n)=U(k+1:n,k+1:n)-L(k+1:n,k)*U(k,k+1:n);
    U(k+1:n,k)=zeros(n-k,1);
end
for j=1:n-1
    Y1(j)=Y1(j)/L(j,j);
    Y1(j+1:n)=Y1(j+1:n)-Y1(j)*L(j+1:n,j);
end
Y1(n)=Y1(n)/L(n,n);
for j=n:-1:2
    Y1(j)=Y1(j)/U(j,j);
    Y1(1:j-1)=Y1(1:j-1)-Y1(j)*U(1:j-1,j);
end
Y1(1)=Y1(1)/U(1,1);
for i=1:n-1
    x=X(i):0.01:X(i+1);
    y=Y1(i)/6*(X(i+1)-x).^3/(X(i+1)-X(i))+Y1(i+1)/6*(x-X(i)).^3/(X(i+1)-X(i))+(Y(i)-Y1(i)/6*(X(i+1)-X(i)).^2)*(X(i+1)-x)/(X(i+1)-X(i))+(Y(i+1)-Y1(i+1)/6*(X(i+1)-X(i)).^2)*(x-X(i))/(X(i+1)-X(i));
    plot(X,Y,'*',x,y,'r');
    hold on
end
xlabel('x');
ylabel('y');
syms x;
for i=1:n-1
    y=Y1(i)/6*(X(i+1)-x).^3/(X(i+1)-X(i))+Y1(i+1)/6*(x-X(i)).^3/(X(i+1)-X(i))+(Y(i)-Y1(i)/6*(X(i+1)-X(i)).^2)*(X(i+1)-x)/(X(i+1)-X(i))+(Y(i+1)-Y1(i+1)/6*(X(i+1)-X(i)).^2)*(x-X(i))/(X(i+1)-X(i))
end

对同一例题采用以上方法分别求解：

Lagrange插值：

Newton插值：

Hermite插值：

分段线性插值：

分段三次Hermite插值：

三次样条插值：

解释：
Lagrange插值和Hermite插值所得结果相同的原因是：两者均采用n次多项式插值，同属于代数插值的范畴。