sanwan.m

数值分析常用几种方法

function sanwan(x,y,ydot)
%三次样条插值三弯矩法第二类边界条件；
%x为向量，全部的插值节点；
%y为向量，插值节点处的函数值；
%ydot为向量，端点处的导数值；
%xi为标量，自变量x;
n=length(x);
h=zeros(1,n);a=ones(1,n);r=ones(1,n);
M=zeros(n,1);b=zeros(n,1);
for k=2:n
h(k)=x(k)-x(k-1);
end
for k=2:n-1
    a(k)=h(k+1)/(h(k)+h(k+1));
    r(k)=1-a(k);
    b(k)=6/(h(k)+h(k+1))*((y(k+1)-y(k))/h(k+1)-(y(k)-y(k-1))/h(k));
end
b(1)=6/h(2)*(y(2)-y(1)/h(2)-ydot(1));
b(n)=6/h(n)*(ydot(2)-(y(n)-y(n-1))/h(n));
A=diag(2*ones(1,n));
for i=1:n-1
    A(i,i+1)=a(i);A(i+1,i)=r(i+1);
end
M=A\b;
for k=2:n
    t=x(k-1):0.1:x(k);
    s=zeros(1,length(t));
    for j=1:length(t)
        s(j)=M(k-1)/6/h(k)*(x(k)-t(j))^3+M(k)/6/h(k)*(t(j)-x(k-1))^3+1/h(k)*(y(k)-M(k)*h(k)^2/6)*(t(j)-x(k-1))+1/h(k)*(y(k-1)-M(k-1)*h(k)^2/6)*(x(k)-t(j));      
    end  
    plot(t,s,'k-');
    hold on
end
 x1=[-5:0.1:5];    
 y1=1./(1+x1.^2); 
 figure(1)
 plot(x,y,'*',x1,y1,'r-');
 legend('插值函数','原函数')  
 title('三弯矩法三次样条与原函数曲线')