Изучение численного метода интерполяции кубическими сплайнами с различными видами граничных условий, страница 6

function M=Progon(A,n)

alfa(1)= - A(1,2)/A(1,1);

beta(1)=A(1,n-1)/A(1,1);

if n>5

    for i= 2 : n-3

        alfa(i)=-A(i,i+1)/(A(i,i-1)*alfa(i-1)+A(i,i));

        beta(i)=(A(i,n-1)-A(i,i-1)*beta(i-1))/(A(i,i-1)*alfa(i-1)+A(i,i));

    end;

end;

M(n-2)=(A(n-2,n-1)-A(n-2,n-3)*beta(n-3))/(A(n-2,n-3)*alfa(n-3)+A(n-2,n-2));

for i=n-3:1;

    M(i)=alfa(i)*M(i+1)+beta(i);

end;

Функции выполняют построение столбца , с учетом разных граничных условий, построение кубического полинома, и построение графика сплайна.

Смыкающийся сплайн:

function G=Splain(x,y,Deriv)

for i=1:length(x)-1

    h(i)=x(i+1)-x(i);

end;

for i=1:length(y)-1

    d(i)=(y(i+1)-y(i))/h(i);

end;

for i=1:length(d)-1

    u(i)=6*(d(i+1)-d(i));

end;

M1=Progon(koefSplain(h, u, d, Deriv,length(x)),length(x));

n=length(x);

M(1)=3*(d(1)-Deriv(1))/h(1)-M1(1)/2;

M(n)=3*(Deriv(2)-d(n-1))/h(n-1)-M1(n-2)/2;

for i =2:n-1

    M(i)=M1(i-1);

end;

for i=1:length(x)-1

    G(i,1)=(M(i+1)-M(i))/(6*h(i));

    G(i,2)=M(i)/2;

    G(i,3)=d(i)-h(i)*(2*M(i)+M(i+1))/6;

    G(i,4)=y(i);

end;

M1=Progon(koefSplain1(h, u, d,length(x)),length(x));

n=length(x);

M(1)=M1(1)-h(1)*(M1(2)-M1(1))/h(2);

M(n)=M1(n-2)+h(n-1)*(M1(n-2)-M1(n-3))/h(n-2);

for i =2:n-1

    M(i)=M1(i-1);

end;

for i=1:length(x)-1

    g(i,1)=(M(i+1)-M(i))/(6*h(i));

    g(i,2)=M(i)/2;

    g(i,3)=d(i)-h(i)*(2*M(i)+M(i+1))/6;

    g(i,4)=y(i);

end;

x1=x(1):0.01:x(2);

x2=x(2):0.01:x(3);

x3=x(3):0.01:x(4);

y1=polyval(g(1,:),x1-x(1));

y2=polyval(g(2,:),x2-x(2));

y3=polyval(g(3,:),x3-x(3));

plot(x1,y1,x2,y2,x3,y3,x,y,'ko');

grid on;

hold off;

title('График смыкающегося сплайна');

end

Естественный сплайн:

function g=Splain3(x,y)

for i=1:length(x)-1

    h(i)=x(i+1)-x(i);

end;

for i=1:length(y)-1

    d(i)=(y(i+1)-y(i))/h(i);

end;

for i=1:length(d)-1

    u(i)=6*(d(i+1)-d(i));

end;

M1=Progon(koefSplain3(h, u, d,length(x)),length(x));

n=length(x);

M(1)=0;

M(n)=0;

for i =2:n-1

    M(i)=M1(i-1);

end;

for i=1:length(x)-1

    g(i,1)=(M(i+1)-M(i))/(6*h(i));

    g(i,2)=M(i)/2;

    g(i,3)=d(i)-h(i)*(2*M(i)+M(i+1))/6;

    g(i,4)=y(i);

end;

x1=x(1):0.01:x(2);

x2=x(2):0.01:x(3);

x3=x(3):0.01:x(4);

y1=polyval(g(1,:),x1-x(1));

y2=polyval(g(2,:),x2-x(2));

y3=polyval(g(3,:),x3-x(3));

plot(x1,y1,x2,y2,x3,y3,x,y,'ko');

grid on;

hold off;

title('График естественного сплайна');

end

Екстраполяционный сплайн:

function g=Splain1(x,y)

for i=1:length(x)-1

    h(i)=x(i+1)-x(i);

end;

for i=1:length(y)-1

    d(i)=(y(i+1)-y(i))/h(i);

end;

for i=1:length(d)-1

    u(i)=6*(d(i+1)-d(i));

end;

M1=Progon(koefSplain1(h, u, d,length(x)),length(x));

n=length(x);

M(1)=M1(1)-h(1)*(M1(2)-M1(1))/h(2);

M(n)=M1(n-2)+h(n-1)*(M1(n-2)-M1(n-3))/h(n-2);

for i =2:n-1

    M(i)=M1(i-1);

end;

for i=1:length(x)-1

    g(i,1)=(M(i+1)-M(i))/(6*h(i));

    g(i,2)=M(i)/2;

    g(i,3)=d(i)-h(i)*(2*M(i)+M(i+1))/6;

    g(i,4)=y(i);

end;

x1=x(1):0.01:x(2);

x2=x(2):0.01:x(3);

x3=x(3):0.01:x(4);

y1=polyval(g(1,:),x1-x(1));

y2=polyval(g(2,:),x2-x(2));

y3=polyval(g(3,:),x3-x(3));

plot(x1,y1,x2,y2,x3,y3,x,y,'ko');

grid on;

hold off;

title('График экстраполяционного сплайна');

end

Сплайн заканчивающийся параболой:

function G=Splain2(x,y)

for i=1:length(x)-1

    h(i)=x(i+1)-x(i);

end;

for i=1:length(y)-1

    d(i)=(y(i+1)-y(i))/h(i);

end;

for i=1:length(d)-1

    u(i)=6*(d(i+1)-d(i));

end;

M1=Progon(koefSplain2(h, u, d, length(x)),length(x));

n=length(x);

M(1)=M1(1);

M(n)=M1(n-2);

for i =2:n-1

    M(i)=M1(i-1);

end;

for i=1:length(x)-1

    G(i,1)=(M(i+1)-M(i))/(6*h(i));

    G(i,2)=M(i)/2;

    G(i,3)=d(i)-h(i)*(2*M(i)+M(i+1))/6;

    G(i,4)=y(i);

end;

x1=0:0.01:1;

x2=1:0.01:2;

x3=2:0.01:3;

y1=polyval(G(1,:),x1-x(1));

y2=polyval(G(2,:),x2-x(2));

y3=polyval(G(3,:),x3-x(3));

plot(x1,y1,x2,y2,x3,y3,x,y,'ko');

grid on;

hold off;

title('График сплайна,заканчивающегося парабалой');

end