题目：牛顿插值&第一边界条件下三次样条插值

拉格朗日插值函数

function [ c,l ] = lagrancz( x,y)
%本函数的功能得到 n 次拉格朗日插值多项式的系数
%x 为 n+1 个节点的横坐标所组成的向量,y 为纵坐标所组成的向量
%c 为所得插值函数的系数组成的向量
%例如,在命令窗口输入 x=[1 2.25 -3];y=[0 2.25 0];lagran(x,y)则输出结果 ans 为
%1 2 -3,即所求的多项式为 x^2+2x-3
w=length(x);
l=zeros(w,w);
for k=1:w
v=1;
 for j=1:w
if k~=j
 v=conv(v,poly(x(j)))/(x(k)-x(j));
 end
 end
 l(k,:)=v;
end
c=y*l;
end

调用拉格朗日插值函数并绘图

x=[0 1 4 9 16 25 36 49 64];
y=[0 1 2 3 4 5 6 7 8];
[c,l]=lagrancz(x,y);
t=length(c);
s=0;
z=0:0.2:64;for k=1:t
 s1=c(k)*z.^(t-k);
 s=s+s1;
end
plot(x,y,’’ob ,z,s,’r’); 

第一边界条件下的三次样条插值求 M 向量函数

% Function：求解第一类边界条件下三弯矩方程组对应的矩阵形式的 M A d h
function [M,h] =get_M1(x,y,diff_y0,diff_yn)
[n,~] = size(x);
n = n-1;
% 自变量相邻两点差值向量
h = diff(x);
% 初始化三对角矩阵 A
A = 2 * eye(n+1,n+1);
% 求解三对角矩阵 A 的下对角参数 μ
u = zeros(n,1);
u(n) = 1;
for i = 1:(n-1)
u(i) = h(i) / (h(i) + h(i+1));
end
% 求解三对角矩阵 A 的上对角参数 λ
lambda = zeros(n,1);
lambda(1) = 1;
for i = 2:n
lambda(i) = h(i) / (h(i-1) + h(i));
end
% 求解矩阵方程右端向量 d
d = zeros(n+1,1);
% d(1) = 6 / h(1) * ((y(2) - y(1))/(x(2) - x(1)) - diff_y0);
d(1) = 6 / h(1) * ((y(2) - y(1))/h(1) - diff_y0);
% d(n+1) = 6 / h(n) * (diff_yn - ((y(end) - y(end-1))/(x(end) - x(end-1))));
d(n+1) = 6 / h(n) * (diff_yn - ((y(end) - y(end-1))/h(n)));
for i = 2:n
d(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(i1)));
end
% 求解三对角矩阵 A
A = A + diag(u,-1) + diag(lambda,1);% 矩阵方程左乘 A 的逆矩阵求解出 M
M = A^(-1)*d;
end

调用求 M 向量函数，绘制图像

x=[0;1;4;9;16;25;36;49;64];
y=[0;1;2;3;4;5;6;7;8];
diff_y0=1; %左端点导数值
diff_yn=1/15; %右端点导数值
 % M：三次样条插值函数中的弯矩
 % h：自变量相邻两点差值向量
 % 求解 M h
[M,h] =get_M1(x,y,diff_y0,diff_yn);
% 将 M、h 带入到三次样条插值函数表达式中，
%从而得出各小区间上的三次样条插值函数
t = 0:0.5:64;
[~,m] = size(t);
S=zeros(1,m);
F=zeros(1,m);
[n,~] = size(x);
n=n-1;
for i=1:n
 tt=(t >=x(i)) & (t <= x(i+1));
 S=( M(i)*((x(i+1) - t).^3./(6*h(i))) ...
 + M(i+1)*((t - x(i)).^3./(6*h(i))) ...
+ (y(i) - M(i)*(h(i)^2)/6).*((x(i+1) - t)/h(i)) ...
 + (y(i+1) - M(i+1)*(h(i)^2)/6).*((t - x(i))./h(i)) ).* tt;
 for j = 1:m %
if tt(j)
 F(j) = S(j);
 end
 end
end
figure;
plot(x,y,’b+’);
hold on
plot(t,F,’r’);
hold off
xlabel(‘x’),ylabel(‘y’);
title(‘Spline 插值’),legend(‘original points’,’spline points’);
grid on;

分析：显然，八次拉格朗日插值龙格现象严重，在[0,64]区间上三次样条插值更准确

改写程序，在[0,1]区间绘图结果

分析：观察图形斜率变化可知在[0,1]区间上拉格朗日插值更准确，可分别取抽样点均方
误差验证。