数值分析实验(一)之插值方法

实验目的

通过设计、编制、调试2~3个多项式插值、拟合曲线的程序，加深对其数值计算方法及有关的基础理论知识的理解。
实现拉格朗日（Lagrange）插值多项式、牛顿（Newton）插值、用线性函数拟合给定数据的程序。

问题
已知插值节点序列，用拉格朗日（Lagrange）插值多项式、牛顿（Newton）插值、用线性函数计算的函数在点的近似值。
我们选取函数y=3x+5
样本点X=0,1,2,3,4,5,6,7,8,9；Y=5,8,11,14,17,20,23,26,29,32

拉格朗日（Lagrange）插值多项式

lagrange.m 文件
用于编写lagrange插值函数

%lagrange.m 
%拉格朗日插值及其误差估计
%输入的量:X是n+1个节点(x_i,y_i)(i = 1,2, ... , n+1)横坐标向量，Y是纵坐标向量， 
%x是以向量形式输入的m个插值点，M在[a,b]上满足｜f~(n+1)(x)｜≤M %注：f~(n+1)(x)表示f(x)的n+1阶导数 
%输出的量：y为m个插值构成的向量，R是误差限 
function [y, R] = lagrange(X, Y, x, M)
n = length(X);
m = length(x); 
for i = 1:m 
    z = x(i); 
    s = 0.0; 
    for k = 1:n 
        p = 1.0; 
        q1 = 1.0; 
        c1 = 1.0; 
        for j = 1:n 
            if j~=k 
                p = p * (z - X(j)) / (X(k) - X(j)); 
            end
            q1 = abs(q1 * (z - X(j)));
            c1 = c1 * j; 
        end
        s = p * Y(k) + s; 
    end
    y(i) = s; 
    R(i) = M * q1 / c1; 
end

sy1_1.m文件
用于运行求解得出结果

X = [0 1 2 3 4 5 6 7 8 9];
Y = [5 8 11 14 17 20 23 26 29 32]; 
x0 = 4.5; 
x = linspace(0,10,50);
M = 1; 
[y0, R] = lagrange(X, Y, x0, M); 
y = 3*x+5; 
errorbar(x0,y0,R,'.g')
hold on 
plot(X, Y, 'or', x0, y0, '.k', x, y, '-b'); 
legend('误差','样本点','拉格朗日插值估算','y=3x+5');

得出结果如图：

牛顿（Newton）插值

newton.m 文件
用于编写牛顿插值函数

function yi=newton(x,y,xi)
%Newton基本插值公式
%x为向量，全部的插值节点
%y为向量，差值节点处的函数值
%xi为标量，是自变量
%yi为xi出的函数估计值
n=length(x);
m=length(y);
%计算均差表Y
Y=zeros(n);
Y(:,1)=y';
for k=1:n-1
    for i=1:n-k
        Y(i,k+1)=(Y(i+1,k)-Y(i,k))/(x(i+k)-x(i));
    end
end
%计算牛顿插值公式
yi=0;
for i=1:n
    z=1;
    for k=1:i-1
        z=z*(xi-x(k));
    end
    yi=yi+Y(1,i)*z;
end

sy1_2.m文件
用于运行求解得出结果

X = [0 1 2 3 4 5 6 7 8 9];
Y = [5 8 11 14 17 20 23 26 29 32]; 
x0 = 4.5; 
x = linspace(0,10,50);
y0 = newton(X, Y, x0); 
y = 3*x+5; 
R = y0-3*x0-5;
errorbar(x0,y0,R,'.g')
hold on 
plot(X, Y, 'or', x0, y0, '.k', x, y, '-b'); 
legend('误差','样本点','牛顿插值估算','y=3x+5');

得出结果如图：

线性函数拟合

sy1_3.m文件
用于运行求解得出结果

X = [0 1 2 3 4 5 6 7 8 9];
Y = [5 8 11 14 17 20 23 26 29 32]; 
x0 = 4.5;
x = linspace(0,10,50);
y = 3*x+5;
y1 = polyfit(X,Y,1);
y0 = x0*y1(1)+y1(2);
R = y0-3*x0-5;
errorbar(x0,y0,R,'.g')
hold on 
plot(X, Y, 'or', x0, y0, '.k', x, y, '-b'); 
legend('误差','样本点','线性函数拟合','y=3x+5');

得出结果如图：