naturalsp1.m

数值计算实践是我们必备的一门学科！这是我在学习之余利用MATLAB编写的程序

function [m,H,lambda,mu,D,A,dY,sk]=Naturalsp(X,Y)
%X为节点(xi,yi)的横坐标;
%Y为节点(xi,yi)的纵坐标;
%m为X的维数;
%H,lambda,mu,D分别自然样条插值函数创建中的向量;
%A是系数矩阵;
%dY是线形方程组的解向量;
%sk是自然样条函数.
m=length(X);A=zeros(m,m);n=m-1;                  
H=zeros(1,n);lambda=zeros(1,n);
mu=zeros(1,n);lambda(1)=1;A(1,1)=2;A(1,2)=lambda(1);
D=zeros(1,n);H(1)=X(2)-X(1);mu(1)=1;D(1)=3*(Y(2)-Y(1));
for k=1:n
    hk=X(k+1)-X(k);H(k+1)=hk;
end
H=H(2:n+1);
for k=1;n-1
    lambdak=H(k)/(H(k)+H(k+1));lambda(k+1)=lambdak;
    muk=1-lambda(k+1);mu(k)=muk;
    dk=3*((mu(k).*(Y(k+1)-Y(k))./H(k))+(lambda(k+1).*(Y(k+2)-Y(k+1))./H(k+1)));
    D(k+1)=dk;
end
D(m)=3*(Y(m)-Y(m-1))/H(m-1);mu(n)=1;n;H;lambda;mu;D;
for i=1:m-1
    A(i,i)=2;A(m,m)=2;A(i,i+1)=lambda(i);
    A(i+1,i)=mu(i);
end
dY=A\D';     %求出线形方程组的解向量;
syms x
m=length(X);S=zeros(m-1,1);
for k=2:m
   
sk=Y(k-1)*((H(k-1)-2*X(k-1)+2*x)*(x-X(k))^2)/(H(k-1)^3)+Y(k)*((H(k-1)+2*X(k)
-2*x)*(x-X(k-1))^2)/(H(k-1)^3)+dY(k-1)*((x-X(k-1))*(x-X(k))^2)/(H(k-1)^2)
+dY(k)*((x-X(k))*(x-X(k-1))^2)/(H(k-1)^2)
%生成自然样条函数;
end