spline.m

数值分析---部分算法 希望对需要的人有一点点用！

% -------------------------------------------------------------------------
% ---------------------------三次样条插值算法-------------------------------
% 函数 SX=Spline(matrix)
% 功能： 通过插值节点得到(n-1)个段上的三次样条插值多项式
% 输入参数：matrix  存有插值节点数据的矩阵,3行n列，n即为插值节点的个数
%                  其第一行为插值点Xi，第二行为函数在x＝Xi时的值f(Xi),
%                  第三行函数一阶倒数值(用到两端点值)
% 输出参数：  SX    1行(n-1)列的符合矩阵，即(n-1)个段上的三次样条插值多项式
%                  关于x的表达式
% 采用高斯列主元素消元法求解线性方程组
% -------------------------------------------------------------------------                 
function SX=Spline(matrix)
syms x SX
[m,n]=size(matrix);
B=zeros(4*n-4,4*n-4);
C=zeros(4*n-4,1);
for i=1:n-1                %注：Xn(即第n个插值节点)满足第(n-1)段的多项式。不注意会弄出n个多项式！！！
    t=sym(matrix(1,i));
    B(i,4*i-3:4*i)=[t^3,t^2,t,1];
    C(i,1)=matrix(2,i);
end
    t=sym(matrix(1,n));
    B(n,4*n-7:4*n-4)=[t^3,t^2,t,1];
    C(n,1)=matrix(2,n);
    t=sym(matrix(1,1));
    B(n+1,1:4)=[3*t^2,2*t,1,0];
    t=sym(matrix(1,n));
    B(n+2,4*n-7:4*n-4)=[3*t^2,2*t,1,0];
    C(n+1,1)=matrix(3,1);
    C(n+2,1)=matrix(3,n);
for i=n+3:2*n
    j=i-n-1;
    k=i-n-3;
    t=sym(matrix(1,j));
    B(i,4*k+1:4*k+4)=[t^3,t^2,t,1];
    B(i,4*k+5:4*k+8)=[-t^3,-t^2,-t,-1];
end
for i=2*n+1:3*n-2
    j=i-2*n+1;
    k=i-2*n-1;
    t=sym(matrix(1,j));
    B(i,4*k+1:4*k+4)=[3*t^2,2*t,1,0];
    B(i,4*k+5:4*k+8)=[-3*t^2,-2*t,-1,0];
end
for i=3*n-1:4*n-4
    j=i-3*n+3;
    k=i-3*n+1;
    t=sym(matrix(1,j));
    B(i,4*k+1:4*k+4)=[6*t,2,0,0];
    B(i,4*k+5:4*k+8)=[-6*t,-2,0,0];
end
% A=(B\C).';   %注：矩阵法求解
% j=1;
% for i=1:4:4*n-4
%     digits(6),SX(j)=vpa(simple(A(i)*x^3+A(i+1)*x^2+A(i+2)*x+A(i+3)));
%     j=j+1;
% end
B(1:4*n-4,4*n-3)=C(1:4*n-4,1);   %注：高斯消元列主元素法求解
A=Solution(B);
j=1;
for i=1:4:4*n-4
    digits(6),SX(j)=vpa(simple(A(i)*x^3+A(i+1)*x^2+A(i+2)*x+A(i+3)));
    j=j+1;
end
SX=simple(SX);