sy_scyt.m

实验一：三次样条插值（P56

%
m=21;n=11;k=[0:10];
x=[-5,-4.8,-4.5,-4.3,-4,-3.8,-3.5,-3.3,-3,-2.8,-2.5,-2.3,-2,-1.8,-1.5,-1.3,-1,-0.8,-0.5,-0.3,0];

% f(x)
for i=1:1:m
    f(i)=1/(1+(x(i))^2);
end

% Ln(x)
l=ones(n,m);
Ln=zeros(1,m);
for i=1:1:n
    xk(i)=k(i)-5;
    fx(i)=1/(1+(xk(i))^2);
end
for k=1:1:m
    for i=1:1:n
        for j=1:1:n
            if i~=j
                l(i,k)=l(i,k)*(x(k)-xk(j))/(xk(i)-xk(j));
            end        
        end
    end
end
for i=1:1:m
    for j=1:1:n
        Ln(i)=Ln(i)+l(j,i)*fx(j);
    end
end

% Sn(x)
h=1;
mn=zeros(n,n);
v=ones(1,n-1);u=ones(1,n-1);
M=zeros(n,1);d=zeros(n,1);
v(1)=1;u(n-1)=1;
for i=2:1:(n-1)
    v(i)=0.5;
    u(i-1)=0.5;
end
d(1)=6/h*( ( fx(2)-fx(1) ) / ( xk(2)-xk(1) ) - 10/(26)^4 );
d(n)=6/h*( -( fx(n)-fx(n-1) ) / ( xk(n)-xk(n-1) ) - 10/(26)^4 );
for i=2:1:(n-1)
    d(i)=6/2/h*( ( fx(i+1)-fx(i) ) / ( xk(i+1)-xk(i) ) -( fx(i)-fx(i-1) ) / ( xk(i)-xk(i-1) ) );
end
for i=1:1:n
    for j=1:1:n
    if i==j
        mn(i,j)=2;
    end
    if i==(j+1)
        mn(i,j)=u(i-1);   
    end
    if j==(i+1)
        mn(i,j)=v(i);   
    end
    end
end 
M=(mn)^(-1)*d;
for i=1:1:5
    for j=1:1:4
        Sn(i*4-4+j)=M(i)*(xk(i+1)-x(i*4-4+j))^3/6/h + M(i+1)*(x(i*4-4+j)-xk(i))^3/6/h + (fx(i)-M(i)*h^2/6)*(xk(i+1)-x(i*4-4+j))/h + (fx(i+1)-M(i+1)*h^2/6)*(x(i*4-4+j)-xk(i))/h;
    end
end
Sn(21)=M(5)*(xk(6)-x(21))^3/6/h + M(6)*(x(21)-xk(5))^3/6/h + (fx(5)-M(5)*h^2/6)*(xk(6)-x(21))/h + (fx(6)-M(6)*h^2/6)*(x(21)-xk(5))/h;

% figure
plot(x,Ln,'b-',x,Sn,'g-',x,f,'r-');