用matlab进行函数插值的几种方法

方法一：Lagrange插值

function [y]=lagrange(x0,y0,x)
    n=length(x0);y=zeros(1,length(x));
    for k=1:n
        t=ones(1,length(x));
        for i=1:n
            if i~=k
                t=t.*(x-x0(i))/(x0(k)-x0(i));
            end
        end
        y=y+t*y0(k);
    end
end

 

 

方法二：Newton插值

function [y]=newton(x0,y0,x)
    n=length(x0);f=zeros(n,n);
    for k=1:n      
        f(k,1)=y0(k);
    end
    for i=2:n
        for k=i:n
            f(k,i)=(f(k,i-1)-f(k-1,i-1))/(x0(k)-x0(k+1-i));  
        end
    end
    y=f(1,1);          
    for k=2:n
        t=ones(1,length(x));
        for j=1:k-1
            t=t.*(x-x0(j)); 
        end
        y=f(k,k)*t+y;
    end
end

 

 

方法二*：牛顿前插公式

function [y]=newton_forward(x0,h,xn,Y,x)
    n=(xn-x0)/h+1;f=zeros(n,n);
    for k=1:n      
        f(k,1)=Y(k);
    end
    for i=2:n
        for k=i:n
            f(k,i)=f(k,i-1)-f(k-1,i-1);  
        end
    end
    t=(x-x0)/h; y=f(1,1);    
    for k=2:n
        p=ones(1,length(x));
        for j=1:k-1
            p=p.*(t-j+1)/j;
        end
        y=f(k,k)*p+y;
    end
end
    

 

方法三：分段线性插值

function [y]=linear(x0,y0,x)
    n=length(x0);m=length(x);y=zeros(1,m);
    for t=1:m
        for i=1:n-1
            if (x(t)<=x0(i+1))&&(x(t)>=x0(i))
              y(t)=y0(i)*(x(t)-x0(i+1))/(x0(i)-x0(i+1))+y0(i+1)*(x(t)-x0(i))/(x0(i+1)-x0(i));
            end
        end
    end
end

 

 

方法四：Hermite插值

function [y]=hermite(x0,y0,m,x)
    n=length(x0);h=length(x);y=zeros(1,h);
    for t=1:h
        for i=1:n-1
            if (x(t)<=x0(i+1))&&(x(t)>=x0(i))
              y(t)=y0(i)*(1+2*((x(t)-x0(i))/(x0(i+1)-x0(i))))*((x(t)-x0(i+1))/(x0(i)-x0(i+1)))^2 ...
                  +y0(i+1)*(1+2*((x(t)-x0(i+1))/(x0(i)-x0(i+1))))*((x(t)-x0(i))/(x0(i+1)-x0(i)))^2 ...
                  +m(i)*(x(t)-x0(i))*((x(t)-x0(i+1))/(x0(i)-x0(i+1)))^2 ...
                  +m(i+1)*(x(t)-x0(i+1))*((x(t)-x0(i))/(x0(i+1)-x0(i)))^2;
            end
        end
    end
end

 

 

方法五：三次样条插值（自然边界条件）

function [y]=threesimple(X,Y,M0,x)
     n=length(X); h=zeros(1,n-1);
     for i=1:n-1
         h(i)=X(i+1)-X(i);
     end
     u=zeros(1,n-1);v=zeros(1,n-1);d=zeros(1,n-1);
     for i=2:n-1
         u(i)=h(i-1)/(h(i-1)+h(i)); 
         v(i)=h(i)/(h(i-1)+h(i));
         d(i)=6*(((Y(i+1)-Y(i))/h(i)-(Y(i)-Y(i-1))/h(i-1)))/(h(i-1)+h(i));
     end
     A=zeros(n-2,n-2);b=zeros(n-2,1);M=zeros(1,n);
     for i=1:n-2
         A(i,i)=2; b(i)=d(i+1);
     end
     b(1)=b(1)-u(2)*M0; b(n-2)=b(n-2)-v(n-1)*M0;
     for i=2:n-2
         A(i,i-1)=u(i+1); A(i-1,i)=v(i);
     end
     M(1)=M0;M(n)=M0;M(2:n-1)=A\b;    
     m=length(x); y=zeros(1,m);   
     for t=1:m
        for i=1:n-1
           if (x(t)<=X(i+1))&&(x(t)>=X(i))
              p1=M(i)*(X(i+1)-x(t))^3/(6*h(i));
              p2=M(i+1)*(x(t)-X(i))^3/(6*h(i));
              p3=(Y(i)-M(i)/6*(h(i))^2)*(X(i+1)-x(t))/h(i);
              p4=(Y(i+1)-M(i+1)/6*(h(i))^2)*(x(t)-X(i))/h(i);
              y(t)=p1+p2+p3+p4; 
           end
        end
     end
end

 

使用范例：

clear
clc
%Q1
x0=[25 40 50 60];y0=[95 75 63 54];x=[70];
disp('answer for Q1');
y=lagrange(x0,y0,x);
y=roundn(y,-1);
disp('lagrange interpolation:y(70)=');disp(y);
%Q2
x0=[1994 1995 1996 1997 1998 1999 2000 2001 2002 2003];
y0=[67.052 68.008 69.803 72.024 73.400 72.063 74.669 74.487 74.065 76.777];
disp('answer for Q2');
X=1990:0.01:2005;
Y=newton(x0,y0,X);plot(x0,y0,'o',X,Y);saveas(gcf,'newton','png');
Y=newton_forward(1994,1,2003,y0,X);plot(x0,y0,'o',X,Y);saveas(gcf,'newton_forward','png');
X=1994:0.01:2003;
Y=linear(x0,y0,X);plot(x0,y0,'o',X,Y);saveas(gcf,'linear','png');
m=[1 2 -7 6 2 -2 2 3 4 -3];
Y=hermite(x0,y0,m,X);plot(x0,y0,'o',X,Y);saveas(gcf,'hermite','png');
x=2010;
y=newton(x0,y0,x);y=roundn(y,-3);
disp('lagrange interpolation:y(2010)=');disp(y);
%Q3
X=[0 1 2 3];Y=[3 5 4 1];M0=1;x=0:0.01:3;
[y_lagrange]=lagrange(X,Y,x);
[y_linear]=linear(X,Y,x);
[y_hermite]=hermite(X,Y,[2 1 1.5 0],x);
[y_threesimple]=threesimple(X,Y,M0,x);
plot(X,Y,'r.',x,y_lagrange,'k-',x,y_linear,'k--',x,y_hermite,'k:',x,y_threesimple,'k-.');
legend('original point','lagrange','linear','hermite','threesimple');
saveas(gcf,'plot of the function','png');