spline.m

matlab写的B样条函数的计算

function  spline(X,Y,dY,x0,m)
N = size(X,2);
s0 = dY(1);sN = dY(2);
interval = 0.025;
%disp('x0为插值点:')；
%x0
h = zeros(1,N-1);
for i = 1:N - 1
    h(1,i) = X(i+1) - X(i);
end
d(1,1) = 6*((Y(1,2)-Y(1,1))/h(1,1)-s0)/h(1,1);
d(1,N) = 6*(sN-(Y(1,N)-Y(1,N-1))/h(1,N-1))/h(1,N-1);
for i = 2:N-1
    d(1,i) = 6*((Y(1,i+1)-Y(1,i))/h(1,i)-(Y(1,i)-Y(1,i-1))/h(1,i-1))/(h(1,i)+h(1,i-1));
end
mu = zeros(1,N-1);md = zeros(1,N-1);u = zeros(1,N-1);
md(1,1) = 1;mu(1,N-1) = 1;
for i = 1:N-2
    u(1,i) = h(1,i)/(h(1,i+1)+h(1,i));
    mu(1,i) = u(1,i);
    md(1,i+1) = 1-u(1,i);
end
%for i = 1:N-2
%    md(1,i+1) = 1-u(1,i);
%end
mm = zeros(1,N);
for i = 1:N
    mm(1,i) = 2;
end


if  N-1==length(mu)
  for   i=N-1:-1:1
        mu(i+1)=mu(i);
       end 
end              %将mu设置为N维向量
                        
md(1)=md(1)/mm(1);     d(1)=d(1)/mm(1);
for    i=2:N-1
      mm(i)=mm(i)-mu(i)*md(i-1);
      md(i)=md(i)/mm(i);
      d(i)=(d(i)-mu(i)*d(i-1))/mm(i);
end
d(N)=(d(N)-mu(N)*d(N-1))/(mm(N)-mu(N)*md(N-1));
for     i=N-1:-1:1
      d(i)=d(i)-md(i)*d(i+1);
end
M = zeros(1,N);
M = d;
%fprintf('M为三对角方程的解\n');
%M
%fprintf('\n');
syms t;
digits(m);
for i = 1:N - 1;
    pp(i) = M(i)*(X(i + 1) - t)^3/(6*h(i)) + M(i+1)*(t - X(i))^3/(6*h(i)) + (Y(i) - M(i)*h(i)^2/6)*(X(i+1)-t)/h(i)...
        + (Y(i + 1) - M(i+1)*h(i)^2/6)*(t - X(i))/h(i);
    pp(i) = simplify(pp(i));
    coeff = sym2poly(pp(i));
    if length(coeff)== 3
        tt = coeff(1:3);
        coeff(1:4) = 0;
        coeff(2:4) = tt;
    end
    if length(coeff)== 2
        tt = coeff(1:2);
        coeff(1:4) = 0;
        coeff(3:4) = tt;
    end
    if length(coeff)== 1
        tt = coeff;
        coeff(1:4) = 0;
        coeff(4) = tt;
    end
    if x0 > X(i)&x0<X(i+1)
        L = i;
       % y0 = uint8(coeff(1)*x0^3 + coeff(2)*x0^2 +coeff(3)*x0 + coeff(4));
         y0 = coeff(1)*x0^3 + coeff(2)*x0^2 +coeff(3)*x0 + coeff(4);
    end
%     val = X(i):interval:X(i+1);
%     for k = 1:length(val)
%         fval(k) = coeff(1)*val(k)^3 + coeff(2)*val(k)^2 +coeff(3)*val(k)+ coeff(4);
%     end
%     if mod(i,2)==1
%         plot(val,fval,'r+');
%     else
%         plot(val,fval,'b.');
%     end
%     hold on;
%     clear val fval;
%     ans = sym(coeff,'d');
%     ans = poly2sym(ans,'t');
%    fprintf('三次样条函数S(%d)=',i);
%    pretty(ans);
 end
%fprintf('在区间[%f,%f]\n',X(L),X(L+1));
%fprintf('函数在插值点x0 = %f的值为\n',x0);
%y0;