exp4_4.m

使用matlab软件编制的计算程序

% exp4_4.m --- II型样条
% 这是我们编的II型样条的程序,仅供同学参考.
% II型样条是对两个端点附加二阶导数值的样条,也称曲率调整样条

function end_point_curvature_adjusted_spline

x=[0    1   2    3];
y=[0  0.5   2  1.5];

S = spline_II(x,y,0,0);   % 两个端点二阶导数为零（自然样条）

xx = linspace(0,3,101); 
yy = ppval(S,xx);         % 对分段多项式求值
plot(x,y,'o',xx,yy,'-');
text(0.9,1.5,'natural spline  \rightarrow','Color','r')

%-----------------   II 型样条  -----------------------
function sp = spline_II(X,Y,Dx2_1,Dx2_n)
% sp = spline_II(X,Y,Dx2_1,Dx2_n) -- II型三次样条插值
%         两端点二阶导数已知，也称端点曲率调整样条
% 输入:   X      ---  插值点横坐标(行向量)
%         Y      ---  插值点纵坐标(行向量)
%         Dx2_1  ---  第一个端点的二阶导数值
%         Dx2_n  ---  最后一个端点的二阶导数值
%                     常用Dx2_1=Dx2_n=0即自然样条
% 输出:   sp     ---  样条(即分段多项式,matlab格式)
% 参见 P93 三弯矩法
N=length(X);
M=zeros(1,N); M(1)=Dx2_1; M(N)=Dx2_n;  % M 存二阶导数
H=diff(X);                             % 求差分: H(k) = X(k+1) - X(k)
                                       % H 的维数比 X 少 1
G=diff(Y)./H;
S=zeros(N-1,4);                        % S 存分段多项式的系数

% 计算三对角方程组,B是主对角线,A是下一对角线,C是上一对角线,D是右端项
B=zeros(N-2,1);A=zeros(N-3,1);C=zeros(N-3,1);D=zeros(N-2,1);
for i=1:N-2
    if i ~= 1
	    A(i-1)=H(i);
    end
    if i ~= N-2
        C(i)=H(i+1);
    end
	B(i)=2*(H(i)+H(i+1));
	D(i)=6*(G(i+1)-G(i));
end
D(1)=D(1)-H(1)*M(1);
D(N-2)=D(N-2)-H(N-1)*M(N);

% 解三对角方程组(这里用稀疏矩阵解参见exp3_5.m,也可用追赶法)
A1 =  sparse(1:N-2,1:N-2,B,N-2,N-2)...
    + sparse(2:N-2,1:N-3,A,N-2,N-2)...
    + sparse(1:N-3,2:N-2,C,N-2,N-2);
M(2:N-1) = A1\D;
% [注] 用这种方法求解P94-95的III型样条是方便的,因为不是三对角方程组了.

% 计算分段多项式的四个系数
%   为了与matlab保持一致,分段多项式要写成下面样子
%   Sk(x) = S(k,1)(x-X(k))^3 + S(k,2)(x-X(k))^2 + S(k,3)(x-X(k)) + S(k,4)
for k=0:N-2
    S(k+1,1)=(M(k+2)-M(k+1))/(6*H(k+1));
    S(k+1,2)=M(k+1)/2;
    S(k+1,3)=G(k+1)-H(k+1)*(2*M(k+1)+M(k+2))/6;
    S(k+1,4)=Y(k+1);
end

sp = mkpp(X,S);            % 把分段多项式转化为matlab的格式

%-----------------   function spline_II   END  -----------------------------


% ********** 你的试验 **********
% 【实验一】
% 用P102实验三的数据(见下)来试试,其它做相应改动,再调整一下端点的二阶导数值看看
% x = [0 70 130 210 337 578 776 1012 1142 1462 1841];
% y = [0 57  78 103 135 182 214  244  256  272  275];
%
% ★【实验二】(此题是附加题,有能力的同学可以试试)
% 编 I 型样条的程序,与上面程序相比只是求解的三对角方程组不一样,其余同.
% (你可用matlab自带的程序验证是否正确)