例1 设)(x f 为定义在[0,3]上的函数，有下列函数值表：
且2.0)('0=x f ，1)('3-=x f ，试求区间[0,3]上满足上述条件的三次样条插值函数)(x s
本算法求解出的三次样条插值函数将写成三弯矩方程的形式：
)
()6()()
6()(6)(6)(211123
13
1j 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));
% 为向量d赋值
d = zeros(1, length(h) + 1);
d(1) = 6 * ( (DepVar(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
曲线的图像如图所示：
例2 已知函数值表：
试求在区间[1,5]上满足上述函数表所给出的插值条件的三次自然样条插值函数)(x s
本算法求解出的三次样条插值函数将写成三弯矩方程的形式：
)
()6()()
6()(6)(6)(211123
13
1j 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 代码：
%============================= % 本段代码解决作业题的例2
%============================= clear all clc
% 自变量x 与因变量y 的取值 IndVar = [1, 2, 4, 5]; DepVar = [1, 3, 4, 2];
% 区间长度向量，其各元素为自变量各段的长度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) = 0;
% 为向量λ赋值
lambda = zeros(1, length(h));
lambda(1) = 0;
for i = 2 : length(lambda)
lambda(i) = h(i) / (h(i - 1) + h(i));
end
% 为向量d赋值
d = zeros(1, length(h) + 1);
d(1) = 0;
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) = 0;
% 为矩阵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, polyval(Part_1, 5), ...
['\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=2和x=4两个横轴点作垂线 % line([2, 2], [4.5, 0.5], 'LineStyle', '--'); line([4, 4], [4.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)表达式为：
[][][]??
?
??∈-+-+--∈-+-+----∈-+-+--=.5,4,)4(2)5(375.4)5(375.0,4,2,)2(75.2)4(75.1)2(1875.0)4(0625.0,2,1,)1(125.3)2()1(125.0)(3333x x x x x x x x x x x x x x s
曲线的图像如图所示：
[][][][]???
????∈-+-+--∈-+-+----∈-+-+----∈-+-+--=.53.0,45.0,)45.0(1.9)53.0(3987.8)53.0(1442.2,45.0,39.0,)39.0(1903.11)45.0(417.10)39.0(859.2)45.0(399.2,39.0,30.0,)3.0(9518.6)39.0(1137.6)3.0(5993.1)39.0(4806.3,
30.0,25.0,)25.0(9697.10)3.0(10)25.0(2652.6)(33
33
33x x x x x x x x x x x x x x x x x x x s
试求在区间[0.25,0.53]上满足上述函数表所给出的插值条件的三次自然样条插值函数s(x)
求解出的三次样条插值函数将写成三弯矩方程的形式：
)
()6()()
6()(6)(6)(211123
13
1j 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 --
+
--
+
-+
-=
+++++本题采用和例2基本相同的Matlab 代码，只改变初始条件。 最终，三次自然样条插值函数s(x)表达式为：
曲线的图像如图所示：
求x y =的二次插值式)(2x P ，使：
);144(12)144(,11)121(,10)100(222===P P P
并计算115的近似值并估计误差。
本题采用拉格朗日二次插值法进行计算： 以下为Matlab 代码：
%============================= % 本段代码解决课本第2章习题与思考题第6题 %============================= clear all clc
% 自变量x 与因变量y 的取值 IndVar = [100, 121, 144]; DepVar = [10, 11, 12];
% 构造拉格朗日插值函数 Coefs_1 = DepVar(1) /...
( ( IndVar(1) - IndVar(2) ) * ( IndVar(1) - IndVar(3) ) ); Part_1 = conv( Coefs_1,...
conv( [1, -IndVar(2)], [1, -IndVar(3)] ) ); f_1 = polyval (Part_1, [IndVar(1) : 0.01 : IndVar(3)]);
Coefs_2 = DepVar(2) /...
( ( IndVar(2) - IndVar(1) ) * ( IndVar(2) - IndVar(3) ) );
Part_2 = conv( Coefs_2,...
conv( [1, -IndVar(1)], [1, -IndVar(3)] ) );
f_2 = polyval (Part_2, [IndVar(1) : 0.01 : IndVar(3)]);
Coefs_3 = DepVar(3) /...
( ( IndVar(3) - IndVar(1) ) * ( IndVar(3) - IndVar(2) ) );
Part_3 = conv( Coefs_3,...
conv( [1, -IndVar(1)], [1, -IndVar(2)] ) );
f_3 = polyval (Part_3, [IndVar(1) : 0.01 : IndVar(3)]);
f = f_1 + f_2 + f_3;
plot ([IndVar(1) : 0.01 : IndVar(3)], f, 'LineWidth', 1.25);
% 在样条插值曲线的相应位置标注该段曲线的函数表达式
text(110, polyval(Part_1, 110) + polyval(Part_2, 110) + polyval(Part_3, 110), ... ['\itP_2(x)=', num2str(Coefs_1), '(x-', num2str(IndVar(2)), ...
')(x-', num2str(IndVar(3)), ')+', num2str(Coefs_2), '(x-', num2str(IndVar(1)), ... ')(x-', num2str(IndVar(3)), ')+', num2str(Coefs_3), '(x-', num2str(IndVar(1)), ...
')(x-', num2str(IndVar(2)), ')'], ...
'FontName', 'Times New Roman', 'FontSize', 14)
% 为x 轴和y 轴添加标注
xlabel( '\itx', 'FontName', 'Times New Roman', ... 'FontSize', 14, 'FontWeight', 'bold');
ylabel( '\itP_2(x)', 'FontName', 'Times New Roman', ... 'Rotation', 0, 'FontSize', 14, 'FontWeight', 'bold');
% 求115平方根的估计值
fun115 = polyval(Part_1, 115) + polyval(Part_2, 115) + polyval(Part_3, 115)
% 求115平方根的准确值 sqrt(115)
% 求误差，单位为%
est = 100 * (fun115 - sqrt(115)) / fun115
最终，拉格朗日插值函数)(2x P 的表达式为：
)
121)(100(011858.0)144)(100(022774.0)144)(121(010823.0)(2--+-----=x x x x x x x P 从计算结果可以看出，115的准确值为10.7238，而通过拉格朗日插值法求出的近似值为10.7228，误差为0.0098%
)(2x P 的图像为：