% 编写函数时候的测试

x=[27.7;28;29;30];
y=[4.1;4.3;4.1;3];
b=[3 -4];

x_num=size(x,1);
s=zeros(x_num,4);
diffy=diff(y);%y的差 x_num-1列向量
FirstOrder=diffy./diff(x);%1阶差商
for i=1:length(x)-2
    delt2x(i)=x(i+2)-x(i);%这样的方法会自动生成行向量
end
delt2y=diff(FirstOrder);
SecondOrder=delt2y./(delt2x');%2阶差商

% 获取M
    %计算x的差h x_num-1维列向量
    h=diff(x);
    %计算d x_num维列向量
    d=zeros(x_num,1);
    d(1)=6/h(1)*(FirstOrder(1)-b(1));

    d(2:x_num-1)=6*SecondOrder(1:x_num-2);

    d(end)=6/h(3)*(b(2)-FirstOrder(end));
    %计算λ x_num-1维列向量,计算μ x_num-1维列向量
    lambda=zeros(x_num-1,1);
    lambda(1)=1;
    miu=zeros(x_num-1,1);
    miu(end)=1;
    for i=1:size(miu,1)-1
        miu(i)=h(i)/(h(i)+h(i+1));
    end
    for i=2:size(lambda,1)  
        lambda(i)=h(i)/(h(i-1)+h(i));
    end
    
    %用diag拼起来
    A=diag(miu,-1)+diag(2*ones(x_num,1))+diag(lambda,1);
    %求解方程Ax=d 得到M
    M=A^(-1)*d;
%计算多项式系数
s=zeros(x_num-1,4);
syms t;
for i=1:x_num-1
    coef1=M(i)/6/h(i);
    coef2=M(i+1)/6/h(i);
    coef3=(y(i)-M(i)*h(i)^2/6)/h(i);
    coef4=(y(i+1)-M(i+1)*h(i)^2/6)/h(i);
    p=coef1*(-t+x(i+1))^3+coef3*(x(i+1)-t)+coef2*(t-x(i))^3+coef4*(t-x(i));
    p=expand(p)
    p=collect(p);
    p=vpa(p,6)
    s(i,:)=vpa((sym2poly(p)),100);
end

%
% 函数单元测试

clc,clear;
x=[27.7;28;29;30];
y=[4.1;4.3;4.1;3];
b=[3 -4];
s = cubic_spline_interplotion(x, y, b)
s2=cubic_spline_interpolation(x',y',b)

%
% 画图

clc;clear;
n=10;%numbers of interval
%n=20;
x=[-1:2/n:1];%interpolation data
y=[];%interpolation data
for i=1:1:length(x)
    y=[y 1/(1+25*x(i)^2)];
end
N=my_newton(x,y);
t=linspace(-1,1,1001);%plotting data
Y_Newton=polyval(N,t);
plot(t,Y_Newton,'r-');%newton
hold on;
plot(x,y,'k*');
b=[0.0740,-0.0740];
spline_matrix=cubic_spline_interpolation(x,y,b);
for i=1:1:n%loop to plot the spline func
    tmp_t=linspace(x(i),x(i+1),1001);
    plot(tmp_t,polyval(spline_matrix(i,:),tmp_t),'b-.')
    polyval(spline_matrix(i,:),[y(i) y(i+1)]')
    if i==1
        legend('Newton fit','points','spline fit')
    end
end

%
%
%
%