spline.m

数值分析中的计算样条插值的matlab源程序。

%三弯矩三次样条插值，第一类边界条件
clear;
x=0:pi/20:pi/2;                                  %x为向量,全部的插值节点;
y=sin(x);                                        %y为向量,插值节点处的函数值;
y_d1=cos(0);y_dn=cos(pi/2);                      %y_d1  y_dn，两端的一阶导数值,
n=length(x); ny=length(y);
%x_in=pi/3;                                         %x_in为标量，自变量x;
x_in=[pi/120 pi/110 pi/100 pi/50 pi/40 pi/6 pi/15 pi/4 pi/3];
n_in=length(x_in);
h=zeros(1,n-1);lambda=ones(1,n-1);mu=ones(1,n-1);f=zeros(1,n-1); 
M=zeros(n,1);d=zeros(n,1);
for j=1:n-1 
    f(j)=(y(j+1)-y(j))/(x(j+1)-x(j)) ;
end; 

for j=1:n-1
    h(j)=x(j+1)-x(j);
end
for j=2:n-1
    lambda(j)=h(j)/(h(j-1)+h(j)); mu(j)=h(j-1)/(h(j-1)+h(j));
    d(j)=6*((f(j)-f(j-1))/(h(j-1)+h(j)));
end
d(1)=6/h(1)*(f(1)-y_d1);
d(n)=6/h(n-1)*(y_dn-f(n-1));
A=diag(2*ones(1,n));
for i=1:n
    if i==1
        A(1,2)=1;
        A(2,1)=mu(1);
    elseif i==n
            A(n,n-1)=1;
        else
    A(i,i+1)=lambda(i-1);A(i+1,i)=mu(i);
    end
end
%M=A\d;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%                 追赶法求三弯矩方程的解           %%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
g=zeros(n,1);
w=zeros(n-1,1);
g(1)=d(1)/A(1,1);
w(1)=A(1,2)/A(1,1);
for i=2:n-1
    g(i)=(d(i)-A(i,i-1)*g(i-1))/(A(i,i)-A(i,i-1)*w(i-1));
    w(i)=A(i,i+1)/(A(i,i)-A(i,i-1)*w(i-1));
end
    g(n)=(d(n)-A(n,n-1)*g(n-1))/(A(n,n)-A(n,n-1)*w(n-1));
   M(n)=g(n);
for i=n-1:-1:1
    M(i)=g(i)-w(i)*M(i+1);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

D=zeros(n_in,4);
for i=1:n_in
     for k=1:n-1                                                                           %S为x_in处的函数估计值
        if x(k)<=x_in(i) && x_in(i)<=x(k+1)
         S(i)=M(k)/6/h(k)*((x(k+1)-x_in(i))^3)+M(k+1)/6/h(k)*((x_in(i)-x(k))^3)+1/h(k)*(y(k)-M(k)*h(k)^2/6)*(x(k+1)-x_in(i))+1/h(k)*(y(k+1)-M(k+1)*h(k)^2/6)*(x_in(i)-x(k));
         dd(i)=abs(sin(x_in(i))-S(i));
          D(i,1) = vpa(x_in(i),4);                           %第一列放全部的插值节点
          D(i,2)=vpa(sin(x_in(i)),4);                             %第二列放节点处的函数值
          D(i,3)=vpa(S(i),4);
          D(i,4)=vpa(dd(i),4);
         
        break;
         
       
        end
         
       % display(S);
        %display(dd);
        %return;
    end
end

%D=zeros(n_in,4);
%for k = 1 : n_in
 %  D(k,1) = x_in(k);                           
  % D(k,2)=sin(x_in(k));                             
  % D(k,3)=S(k);
  % D(k,4)=vpa(dd(k),7);
%end
D=vpa(D,4);
display('   所求点       精确值       估计值       误差');
display(D);
 plot(x,y,'-b',x_in,S,'-r'),
 xlabel('x'),
 ylabel('y'),
 title('三弯矩三次样条插值图')