Интерполирование функции при помощи математической системы MatLab, страница 3

Таким образом, на входе программы мы имеем точки заданной кусочно-линейной функции, а на выходе графики функции и погрешностей.

Текст программы

Hx.m

%Равномерное разбиение отрезка

function x=Hx(a,b,n) %a и b – концы отрезка, n – порядок интерполирования

h=(b-a)/n;

for i=1:n+1

    x(i)=a+h*(i-1); %узлы интерполирования

end

Chx.m

%Разбиение отрезка по Чебышеву

function x=Chx(n) %n - порядок интерполирования

for i=1:n+1

    x(i)=cos((2*(n+2-i)-1)*pi/(2*n+2)); %узлы интерполирования

end

Pspline.m

%Интерполирование сплайна

function Spl=Pspline(Sp,t) % Sp – матрица точек заданной функции, t – вектор узлов

k=length(Sp); %Порядок сплайна

x=Sp(1,1:k); %Вектор аргументов

y=Sp(2,1:k); %Вектор значений функции

B(1:k,1)=1;

for i=1:k

    B(i,2)=x(i);

    for j=3:k

        B(i,j)=abs(x(i)-x(j-1));

    end

end

a=B\y'; %Коэффициенты сплайна

Spl=a(1)+a(2)*t;

for i=3:k

    Spl=Spl+a(i)*abs(t-x(i-1)); %Значения сплайна в узлах t

end

RRasn.m

%Вычисление разделенных разностей

function rr=RRasn(n,f,x)

%n – порядок интерполирования, f – вектор значений функции в узлах, x – вектор узлов

for i=1:n+1

    rr(i,1)=f(i);

end

for j=2:n+1

    for i=1:n+2-j

        rr(i,j)=(rr(i+1,j-1)-rr(i,j-1))/(x(i+j-1)-x(i)); %Значения разделенных разностей

    end

end

Newton.m

%Полином Ньютона

function Nn=Newton(n,f,t,x)

%n – степень полинома, f – вектор значений функции в узлах, t – вектор узлов, x – узел

rr=RRasn(n,f,t);

fx=rr(1,1:n+1);

Nn=0; dtx=1;

for i=1:n+1

    Nn=Nn+dtx*fx(i); %Значения полинома в узле x

    dtx=dtx*(x-t(i));

end

GrSp.m

%Построение графика сплайна

load Splin.dat; %Загрузка сплайна

figure;

x=[-1:0.02:1];

y=Pspline(Splin,x);

plot(x,y,'m');

GrPogr.m

%Построение погрешностей для разного количества точек разбиения