Интерполирование, аппроксимация и сглаживание эксперимент данных
Найти приближенное значение функции y, заданной таблицей данных, при указанном значении аргумента, используя:
1. Интерполяционный полином Лагранжа (или любой подходящий для табличной функции).
2. Метод наименьших квадратов 2-го порядка.
3. Сплайн 3-го порядка.
Сделать выводы по работе (наиболее приемлемый способ аппроксимации заданной табличной функции, рациональный порядок функции отклика с учетом наличия/отсутствия эффекта Рунге на краях области определения аргумента).
Решение
Полученная табличная функция:
i
0 1 2 3 4
xi
0,131 0,136 0,141 0,146 0,151
yi
1,35207 1,36773 1,38357 1,39959 1,41579
1. Полином Лагранжа
Px=y0L0+y1L1+y2L2+y3L3+y4L4
Lj=i=04x-xixj-xi, i≠j
Многочлены Lj как и сам многочлен Лагранжа имеет 4-ю степень.
В Matlab создаём m-функцию, которая возвращает значения интерполяционного многочлена Лагранжа в промежуточных точках отрезка по заданным узлам интерполирования:
function yy=lagrange(x,y,xx)
% вычисление интерполяционного полинома в форме Лагранжа
% x - массив координат узлов
% y - массив значений интерполируемой функции
% xx - массив значений аргумента, для которых надо вычислить значения полинома
% yy - массив значений полинома в точках xx
% узнаем число узлов интерполяции (N=n+1)
N=length(x);
% создаем нулевой массив значений интерполяционного полинома
yy=zeros(size(xx));
% в цикле считаем сумму по узлам
for k=1:N
% вычисляем произведения, т.е. функции L_k
t=ones(size(xx));
for j=[1:k-1, k+1:N]
t=t.*(xx-x(j))/(x(k)-x(j));
end
% накапливаем сумму
yy = yy + y(k)*t;
end
2. Метод сплайнов.
Рассмотрим кубический сплайн.
На каждом i–м отрезке [xi–1, xi], i=1, 2,…, N, (N=3) будем искать полином третьей степени по формуле:
Six=ai+bix-xi+cix-xi22+dix-xi36 (2)
Требуется построить интерполянту — функцию f(x), совпадающую с функцией y(x) в узлах сетки:
fxi=yi
Накладываются также условия непрерывности первой и второй производной в узлах:
f'xi-0=f'xi+0, f''xi-0=f''xi+0, i=1,…,n-1
На границе при x0 и xn=x7.
f''x0=f''xn=0
В итоге необходимо решить СЛАУ вида относительно неизвестных в (2) ci:
hici+2hi+hi+1ci+1+hi+1ci+2=3yi+1-yihi+1-yi-yi-1hi, i=1,…,n-1
hi=xi-xi-1
С краевыми условиями
c1=0, cn+1=0
Данную систему можно решить методом прогонки. Соответственно для приведения к виду Aixi-1+Cixi+Bixi+1=Fi примем:
Ai=hi, i=2,…,n,
Bi=hi+1, i=1,…,n-1
Ci=2hi+hi+1, i=1,…,n-1
Fi=3yi+1-yihi+1-yi-yi-1hi, i=1,…,n-1
Таким образом, мы получим значения ci, i=1,…,n
. После чего находим последовательно:
ai=yi-1;
di=ci+1-ci3hi
bi=yi-1;
Реализуем данный алгоритм в MATLAB (отметим, что в MATLAB индексы массива начинаются от 1, а не от нуля, поэтому в данных формулах в тексте программы индексы при xi и yi будут на 1 больше).
% вычисление кубического сплайна
% x - массив координат узлов
% y - массив значений интерполируемой функции
% xx - массив значений аргумента, для которых надо вычислить значения
% yy - массив значений полинома в точках xx
function yy = Spline3(x,y,xx)
n = length(x)-1;
A(1)=0; B(n)=0;
for i=1:n
h(i)=x(i+1)-x(i);
end
for i=2:n
A(i)=h(i);
end
for i=1:(n-1)
B(i)=h(i+1);
end
for i=1:(n-1)
C(i)=2*(h(i)+h(i+1));
F(i)=3*((y(i+2)-y(i+1))/h(i+1)-(y(i+1)-y(i))/h(i));
end
c = zeros(n+1,1);
c(2:n) = Progonka(A,C,B,F);
for k=1:n
a(k)=y(k);
d(k)=(c(k+1)-c(k))/(3*h(k));
b(k)= (y(k+1)-y(k))/h(k)-h(k)*(c(k+1)+2*c(k))/3;
end
% создаем нулевой вектор для записи значений сплайна
yy = zeros(size(xx));
% вычисляем значения сплайна
% в цикле проходим по отрезкам между узлами интерполирования
for k = 1:n
% находим индексы элементов массива xx, значения которых
% лежат между x(k) и x(k+1)
ind = find((xx >= x(k)) & (xx<= x(k+1)));
% вычисляем значения соответствующего полинома
yy(ind) = a(k) + b(k)*(xx(ind)-x(k)) + ...
c(k)*(xx(ind)-x(k)).^2 + d(k)*(xx(ind)-x(k)).^3;
end
end
Метод прогонки задан отдельной функцией:
function [ X ] = Progonka(A,C,B,F)
n = length(C);
%Определяем прогоночные коэффициенты:
alpha(2)= -B(1)/C(1);
beta(2)= F(1)/C(1);
%Остальные коэффициенты определяются по рекуррентным формулам:
for i=2:(n-1)
alpha(i+1)=-B(i)/(A(i)*alpha(i)+C(i));
beta(i+1)=(F(i)-A(i)*beta(i))/(A(i)*alpha(i)+C(i));
end
%находим X(n):
X(n)=(F(n)-A(n)*beta(n))/(A(n)*alpha(n)+C(n));
%Остальные значения определяем по рекуррентной формуле
for j=(n-1):-1:1
X(j)=alpha(j+1)*X(j+1)+beta(j+1);
end
X=X';%ответ в виде вектор-столбца
end
3
Задать вопрос
на новый заказ в Автор24.
