Matlab数值计算（多项式插值）

多项式插值问题

拉格朗日插值多项式

例1：在某个化学反应过程中，在有限个时刻t(min)，测得生成物浓度y( $10^{^{-3}}$ g/ $cm^{3}$ )d的数据如下：

$t_{i}$	1	2	3	4	6	8	10	12	14	16
$y_{i}$	4.00	6.41	8.01	8.79	9.53	9.86	10.33	10.42	10.53	10.61

求在时刻t=5分，t=16.4分时的浓度是多少。
lagrange.m

function y=lagrange(x0,y0,x)
n=length(x0);m=length(x);
for k=1:m
 z=x(k);
 s=0.0;
 for j=1:n
 p=1.0;
 for i=1:n
 if i~=j
 p=p*(z-x0(i))/(x0(j)-x0(i));
 end
 end
 s=p*y0(j)+s;
 end
 y(k)=s;
end

主程序：

x0=[1 2 3 4 6 8 10 12 14 16]; % 给定的x坐标数据点
y0=[4.00 6.41 8.01 8.79 9.53 9.86 10.33 10.42 10.53 10.61]; % 对应的y坐标数据点

x=[1:0.1:16]; % 生成从1到16的步长为0.1的x坐标数据

y=lagrange(x0,y0,x); % 使用拉格朗日插值计算在给定x坐标下的y值

x1=5; % 指定的插值点x1
x2=16.4; % 指定的插值点x2

y1=lagrange(x0,y0,x1); % 在插值点x1处的拉格朗日插值结果
y2=lagrange(x0,y0,x2); % 在插值点x2处的拉格朗日插值结果

plot(x0,y0,'.k','markersize',20) % 绘制原始数据点
hold on
plot(x,y,'-r','markersize',30) % 绘制拉格朗日插值曲线
hold on 
plot(x1,y1,'*b','markersize',8) % 标记插值点x1处的结果
hold on 
plot(x2,y2,'*b','markersize',8) % 标记插值点x2处的结果
legend('原数值点','Lagrange曲线') % 添加图例，说明原数值点和Lagrange曲线

运行效果如下：

例2：已知函数 y=f(x) 的如下离散数据：

x	0	1	2
y	2	3	12

（1）用线性插值求函数在x= 1.5 处的近似值

（2）用二次插值求函数在x= 1.5 处的近似值

（3）将以上信息可视化

% 给定的数据点
x = [0 1 2]; % 给定的x坐标数据点
y = [2 3 12]; % 对应的y坐标数据点

% (1) 线性插值
x_interp = 1.5; % 想要进行插值的x坐标
y_lin_interp = interp1(x, y, x_interp, 'linear'); % 使用interp1函数进行线性插值计算
disp(['(1). 线性插值结果：', num2str(y_lin_interp)]); % 显示线性插值结果

% (2) 二次插值
n = length(x); % 数据点个数
L = zeros(1, n); % 初始化存储拉格朗日基函数的数组
for i = 1:n
    L(i) = prod((x_interp - x([1:i-1,i+1:end]))./(x(i) - x([1:i-1,i+1:end]))); % 计算拉格朗日基函数
end

y_quad_interp = sum(y .* L); % 计算二次插值的结果
disp(['(2). 二次插值结果：', num2str(y_quad_interp)]); % 显示二次插值结果

% 可视化结果
figure;
plot(x, y, 'ro', 'MarkerSize', 8); % 绘制原始数据点
hold on;
plot(x_interp, y_lin_interp, 'bs', 'MarkerSize', 8); % 绘制线性插值结果点
plot(x_interp, y_quad_interp, 'gd', 'MarkerSize', 8); % 绘制二次插值结果点
legend('原始数据点', '线性插值结果', '二次插值结果');
xlabel('x');
ylabel('y');
title('线性插值和二次插值');
grid on;

运行效果：

牛顿插值多项式

用牛顿插值多项式解答上述例2

% 给定的数据点
x = [0 1 2]; % 给定的x坐标数据点
y = [2 3 12]; % 对应的y坐标数据点

% (1) 线性插值
x_interp = 1.5; % 想要进行插值的x坐标
y_lin_interp = interp1(x, y, x_interp, 'linear'); % 使用interp1函数进行线性插值计算
disp(['(1). 线性插值结果：', num2str(y_lin_interp)]); % 显示线性插值结果

% (2) 牛顿插值
n = length(x); % 数据点个数
% 计算差商
f = y;
for i = 2:n
    for j = n:-1:i
        f(j) = (f(j) - f(j-1)) / (x(j) - x(j-i+1));
    end
end

% 使用差商计算插值结果
y_newton_interp = f(n);
for i = n-1:-1:1
    y_newton_interp = f(i) + (x_interp - x(i)) * y_newton_interp;
end
disp(['(2). 牛顿插值结果：', num2str(y_newton_interp)]); % 显示牛顿插值结果

% 可视化结果
figure;
plot(x, y, 'ro', 'MarkerSize', 8); % 绘制原始数据点
hold on;
plot(x_interp, y_lin_interp, 'bs', 'MarkerSize', 8); % 绘制线性插值结果点
plot(x_interp, y_newton_interp, 'gd', 'MarkerSize', 8); % 绘制牛顿插值结果点
legend('原始数据点', '线性插值结果', '牛顿插值结果');
xlabel('x');
ylabel('y');
title('线性插值和牛顿插值');
grid on;

运行结果：

分段线性插值多项式

例3：已知数据

0	1	2	3	4
2	1	4	3	5

用分段线性插值多项式求x=1.5的近似值

% 给定的数据点
x = [0 1 2 3 4]; % 给定的x坐标数据点
y = [2 1 4 3 5]; % 对应的y坐标数据点

% 想要进行插值的x坐标
x_interp = 1.5;

% 寻找 x_interp 所在的区间
index = find(x <= x_interp, 1, 'last');

% 计算分段线性插值
if index < length(x)
    x_left = x(index);
    x_right = x(index + 1);
    y_left = y(index);
    y_right = y(index + 1);
    
    y_interp = y_left + (y_right - y_left) / (x_right - x_left) * (x_interp - x_left);
else
    y_interp = y(end);
end

disp(['分段线性插值结果：', num2str(y_interp)]); % 显示分段线性插值结果

% 可视化结果
figure;
plot(x, y, 'ro', 'MarkerSize', 8); % 绘制原始数据点
hold on;
plot(x_interp, y_interp, 'bs', 'MarkerSize', 8); % 绘制插值点
legend('原始数据点', '插值点');
xlabel('x');
ylabel('y');
title('分段线性插值');
grid on;

运行结果：

三次样条插值

用三次样条插值解决上述例3

% 给定的数据点
x = [0 1 2 3 4]; % 给定的x坐标数据点
y = [2 1 4 3 5]; % 对应的y坐标数据点

% 使用三次样条插值
pp = spline(x, [0 y 0]); % 在端点处增加0值，进行自然边界条件插值

% 想要进行插值的x坐标
x_interp = 1.5;

% 计算插值结果
y_interp = ppval(pp, x_interp);

disp(['三次样条插值结果：', num2str(y_interp)]); % 显示三次样条插值结果

% 可视化结果
xx = linspace(0, 4, 100); % 生成更密集的x坐标
yy = ppval(pp, xx);
figure;
plot(x, y, 'ro', 'MarkerSize', 8); % 绘制原始数据点
hold on;
plot(xx, yy, 'b-'); % 绘制插值曲线
legend('原始数据点', '三次样条插值曲线');
xlabel('x');
ylabel('y');
title('三次样条插值');
grid on;

运行结果：