三次样条插值作业题(教育教学)

例1 设)(x f 为定义在[0,3]上的函数，有下列函数值表：
x i 0 1 2 3 y i
0.5
2
1.5
且2.0)('0=x f ，1)('3-=x f ，试求区间[0,3]上满足上述条件的三次样条插值函数)(x s
本算法求解出的三次样条插值函数将写成三弯矩方程的形式：
)
()6()()
6()(6)(6)(211123131j j j
j j j j
j j j j j
j j j
j x x h h M y x x h h M y x x h M x x h M x s --
+
--
+
-+-=
+++++其中，方程中的系数
j
j h M 6，
j
j h M 61+，j
j j j h h M y )6(2-
，
j
j
j j h h M y )
6(211++-
将由Matlab
代码中的变量Coefs_1、Coefs_2、Coefs_3以及Coefs_4的值求出。

以下为Matlab 代码：
%=============================
% 本段代码解决作业题的例1
%============================= clear all clc
% 自变量x 与因变量y ，两个边界条件的取值 IndVar = [0, 1, 2, 3]; DepVar = [0, 0.5, 2, 1.5];
LeftBoun = 0.2; RightBoun = -1;
% 区间长度向量，其各元素为自变量各段的长度 h = zeros(1, length(IndVar) - 1); for i = 1 : length(IndVar) - 1
h(i) = IndVar(i + 1) - IndVar(i); end
% 为向量μ赋值
mu = zeros(1, length(h));
for i = 1 : length(mu) - 1
mu(i) = h(i) / (h(i) + h(i + 1));
end
mu(i + 1) = 1;
% 为向量λ赋值
lambda = zeros(1, length(h));
lambda(1) = 1;
for i = 2 : length(lambda)
lambda(i) = h(i) / (h(i - 1) + h(i));
end
% 为向量d赋值
d = zeros(1, length(h) + 1);
d(1) = 6 * ( (DepV ar(2) - DepVar(1) ) / ( IndVar(2) - IndVar(1) ) - LeftBoun) / h(1);
for i = 2 : length(h)
a = ( DepVar(i) - DepVar(i - 1) ) / ( IndVar(i) - IndVar(i - 1) );
b = ( DepVar(i + 1) - DepVar(i) ) / ( IndVar(i + 1) - IndVar(i) );
c = (b - a) / ( IndVar(i + 1) - IndVar(i - 1) );
d(i) = 6 * c;
end
d(i + 1) = 6 *( RightBoun - ( DepVar(i + 1) - DepVar(i) ) / ( IndVar(i + 1) - IndVar(i) ) ) / h(i);
% 为矩阵A赋值
% 将主对角线上的元素全部置为2
A = zeros( length(d), length(d) );
for i = 1 : length(d)
A(i, i) = 2;
end
% 将向量λ的各元素赋给主对角线右侧第一条对角线
for i = 1 : length(d) - 1
A(i, i + 1) = lambda(i);
end
% 将向量d的各元素赋给主对角线左侧第一条对角线
for i = 1 : length(d) - 1
A(i + 1, i) = mu(i);
end
% 求解向量M
M =A \ d';
% 求解每一段曲线的函数表达式
for i = 1 : length(h)
Coefs_1 = M(i) / (6 * h(i));
Part_1 = conv( Coefs_1, ...
conv( [-1, IndVar(i + 1)], ...
conv( [-1, IndVar(i + 1)], [-1, IndVar(i + 1)] ) ) );
S_1 = polyval (Part_1, [IndVar(i) : 0.01 : IndVar(i + 1)]);
Coefs_2 = M(i + 1)/(6 * h(i));
Part_2 = conv( Coefs_2, ...
conv( [1, -IndVar(i)], ...
conv( [1, -IndVar(i)], [1, -IndVar(i)] ) ) );
S_2 = polyval (Part_2, [IndVar(i) : 0.01 : IndVar(i + 1)]);
Coefs_3 = (DepVar(i) - M(i) * h(i)^2 / 6) / h(i);
Part_3 = conv(Coefs_3, [-1, IndVar(i + 1)]);
S_3 = polyval (Part_3, [IndVar(i) : 0.01 : IndVar(i + 1)]);
Coefs_4 = (DepVar(i + 1) - M(i + 1) * h(i)^2 / 6) / h(i);
Part_4 = conv(Coefs_4, [1, -IndVar(i)]);
S_4 = polyval (Part_4, [IndVar(i) : 0.01 : IndVar(i + 1)]);
S = S_1 + S_2 + S_3 + S_4;
plot ([IndVar(i) : 0.01 : IndVar(i + 1)], S, 'LineWidth', 1.25)
% 在样条插值曲线的相应位置标注该段曲线的函数表达式
text(i - 1, polyval(Part_1, 3), ...
['\itS', num2str(i), '(x)=', num2str(Coefs_1), '(', num2str( IndVar(i + 1) ), '-x)^{3}+', ...
num2str(Coefs_2), '(x-', num2str( IndVar(i) ), ')^{3}+', num2str(Coefs_3), ...
'(', num2str( IndVar(i + 1) ), '-x)+', num2str(Coefs_4), '(x-', num2str( IndVar(i) ), ')'], ...
'FontName', 'Times New Roman', 'FontSize', 14)
hold on
end
% 过x=1和x=2两个横轴点作垂线%
line([1, 1], [2.5, -0.5], 'LineStyle', '--');
line([2, 2], [2.5, -0.5], 'LineStyle', '--');
% 为x轴和y轴添加标注
xlabel( '\itx', 'FontName', 'Times New Roman', ...
'FontSize', 14, 'FontWeight', 'bold');
ylabel( '\its(x)', 'FontName', 'Times New Roman', ...
'Rotation', 0, 'FontSize', 14, 'FontWeight', 'bold');
最终，三次样条插值函数s(x)表达式为：
[][][]⎪⎩
⎪
⎨⎧∈-+-+-+--∈-+-+---∈+-++--=.3,2,)2(44.1)3(62.2)2(06.0)3(62.0,2,1,)1(62.2)2(08.0)1(62.0)2(42.0,1,0,08.0)1(06.042.0)1(06.0)(333333x x x x x x x x x x x x x x x x s
曲线的图像如图所示：。