treespline.m

三次样条插值

clc
clear
%function TreeSpline(x0,y0)
x0=[0.9,1.3,1.9,2.1,2.6,3.0,3.9,4.4,4.7,5.0,6.0,7.0,8.0,9.2,10.5,11.3,11.6,12.0,12.6,13.0,13.3];
y0=[1.3,1.5,1.85,2.1,2.6,2.7,2.4,2.15,2.05,2.1,2.25,2.3,2.25,1.95,1.4,0.9,0.7,0.6,0.5,0.4,0.25];	
   N=length(x0);
   h(1)=x0(2)-x0(1);
   u(1)=1;
   g(1)=3*(y0(2)-y0(1))/h(1);
  % A(1)=-u(1)/2;
   %B(1)=g(1)/2;
   x=[];
   s=[];
   for i=1:N-1
       h(i)=x0(i+1)-x0(i);
   end
   for i=2:N-1
       u(i)=h(i-1)/(h(i-1)+h(i));
       l(i)=h(i)/(h(i)+h(i-1));
       g(i)=3*(l(i)*(y0(i)-y0(i-1))/h(i-1)+u(i)*(y0(i+1)-y0(i))/h(i));
   end
   g0=3*(y0(2)-y0(1))/h(1);
   gN=3*(y0(N)-y0(N-1))/h(N-1);
   Y=zeros(N,N);
   Y=Y+2*eye(N);
   Y(1,1)=2;
   Y(1,2)=1;
   Y(N,N-1)=1;
   for i=2:N-1
       Y(i,i-1)=l(i);
       Y(i,i+1)=u(i);
   end
   g(1)=g0;
   gs=[g,gN];
   m=inv(Y)*gs'
       for i=1:N-1
           xx=linspace(x0(i),x0(i+1),20);
           M=length(xx);
           for k=1:M
               ss(k)=(h(i)+2*(xx(k)-x0(i)))*((xx(k)-x0(i+1))^2)*y0(i)/(h(i)^3)+(h(i)-2*(xx(k)-x0(i+1)))*((xx(k)-x0(i))^2)*y0(i+1)/(h(i)^3)+(xx(k)-x0(i))*((xx(k)-x0(i+1))^2)*m(i)/(h(i)^2)+(xx(k)-x0(i+1))*((xx(k)-x0(i))^2)*m(i+1)/(h(i)^2);
            end
            s=[s,ss];
            x=[x,xx];
         end
     plot(x0,y0,'*',x0,y0,x,s,'--k')