三次样条插值作业题

例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));

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');