数学建模算法（基于matlab和python）之 Lagrange插值、Newton插值(1/10)

实验目的及要求：

1、了解多项式插值公式的存在唯一性条件及其余项表达式的推导。

2、了解拉格朗日插值多项式的构造、计算及其基函数的特点，牛顿插值多项式的构造与应用，差商、差分的计算及基本性质。

实验内容：

1、编写Lagrange插值法、Newton插值法通用子程序，依据数据表

 构造Lagrange、Newton插值多项式，画出函数图像，并计算f(5.137)的近似值并估计误差。

1. 拉格朗日插值的MATLAB主程序

被调用的lagrange.m文件：

function Y=lagrange(x,y,X)

n=length(x); m=length(X);

for i=1:m

    z=X(i);s=0;

    for k=1:n

        p=1;

        for j=1:n

            if j~=k

                p=p*(z-x(j))/(x(k)-x(j));

            end

        end

        s=p*y(k)+s;

    end

    Y(i)=s;

End

调用lagrange.m文件的cha_zhi.m执行文件：

x=0:10;

y=[0 0.79 1.53 2.19 2.71 3.03 3.27 2.89 3.06 3.19 3.29];

X=5.137;

Y=lagrange(x,y,X)

plot(x,y,'b+',X,Y,'r*')

hold on

%semilogy(x,y)%折线图

values = spcrv([[x(1) x x(end)];[y(1) y y(end)]],3);

plot(values(1,:),values(2,:))

%光滑的曲线图图

%关于spcrv用法

%values = spcrv([[t(1) t t(end)];[y(1) y y(end)]],3);

%plot(values(1,:),values(2,:));%直接plot只能生成折线图

%标点

text(0.1,-0.1,'(0, 0)')

text(1.1,0.8,'(1, 0.79)')

text(1.9,1.4,'(2, 1.53)')

text(2.8,2.1,'(3, 2.19)')

text(3.8,2.6,'(4, 2.71)')

text(3.5,3.0,'(5, 3.03)')

text(3.8,3.25,'(5.137, 3.081)')

 text(6.1,3.35,'(6, 3.27)')

 text(6.6,2.7,'(7, 2.89)')

 text(7.7,2.9,'(8, 3.06)')

 text(8.7,3,'(9, 3.19)')

 text(9.75,3.15,'(10, 3.29)')

%关于text用法

%text( , ,'')前部分填写坐标，后部分''内为显示内容

Newton插值

算法分析

算法流程图

牛顿插值的MATLAB主程序

被调用的Newton.m文件：

function Y=Newton(x,y,X)

x=x(:);

y=y(:);

n=length(x);   

m=length(X); 

A=zeros(n,n); 

Y=0*X;      

A(:,1)=y;     

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

%A(1,1)=f(x1),A(2,2)=f(x1,x2), …,A(n,n)=f(x1,x2, …xn)

end

end

%用n-1次Newton插值多项式计算X中每个元素的近似值

for j=1:m

     for k=2:n

         t=1;

         for i=1:k-1

             t=(X(j)-x(i))*t;

          end

          Y(j)=Y(j)+A(k,k)*t;

      end

      Y(j)=Y(j)+A(1,1);

end

调用Newton.m文件的cha_zhi2.m执行文件：

x=0:10;

y=[0 0.79 1.53 2.19 2.71 3.03 3.27 2.89 3.06 3.19 3.29];

X=5.137;

Y=Newton(x,y,X)

plot(x,y,'b+',X,Y,'r*')

hold on

%semilogy(x,y)%折线图

values = spcrv([[x(1) x x(end)];[y(1) y y(end)]],3);

plot(values(1,:),values(2,:))

%光滑的曲线图图

%关于spcrv用法

%values = spcrv([[t(1) t t(end)];[y(1) y y(end)]],3);

%plot(values(1,:),values(2,:));%直接plot只能生成折线图

%标点

text(0.1,-0.1,'(0, 0)')

text(1.1,0.8,'(1, 0.79)')

text(1.9,1.4,'(2, 1.53)')

text(2.8,2.1,'(3, 2.19)')

text(3.8,2.6,'(4, 2.71)')

text(3.5,3.0,'(5, 3.03)')

text(3.8,3.25,'(5.137, 3.081)')

 text(6.1,3.35,'(6, 3.27)')

 text(6.6,2.7,'(7, 2.89)')

 text(7.7,2.9,'(8, 3.06)')

 text(8.7,3,'(9, 3.19)')

 text(9.75,3.15,'(10, 3.29)')

%关于text用法

%text( , ,'')前部分填写坐标，后部分''内为显示内容

1. 误差估计（事后估计法）

由以上实验已经求得近似值Y=3.080849054732755，设再取x4，x6为节点，用线性插值可以求得另外一个插值Y1，再由书第21页公式（15）可以得到误差值。

x=0:10;

y=[0 0.79 1.53 2.19 2.71 3.03 3.27 2.89 3.06 3.19 3.29];

X=5.137;

Y=lagrange(x,y,X)

%Y1=y4+((y6-y4)/(x6-x4))/(5.137-x4)

Y1=2.71+((3.27-2.71)/(6-4))*(5.137-4)

%Z=((X-x5)/(x6-x5))*(Y1-Y)

Z=((X-5)/(6-5))*(Y1-Y)

得到如下结果：

所以误差为：-0.007191000498387

可以得到新的近似值：Y2=Y+Z=3.073658054234368

运行结果：

1.拉格朗日插值法的运行结果:

 2.牛顿插值法的运行结果

结果分析与讨论：

1. 一般来说，如果我们有n个点（x1,y1）...(xn,yn),其中xi各不相同，则可以得到拉格朗日基函数。从而定义n-1次多项式：  
2. 差商具有对称性，任意调换节点的次序，差商的值不变。
3. 牛顿插值具有承袭性。当增加节点和插值节点的信息时，不需要像拉格朗日方法那样完全重新计算拉格朗日系数，只需要在原来插值函数的基础上增加高次项便可得到，计算量远小于拉格朗日插值方法。