% 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自帶的程序驗證是否正確)