Лекция Интерполация и прекарване на крива през точки

Лекция 6. Интерполация и прекарване на крива през точки

Известно е, че през три точки може да се прокара крива (полином) от втора степен. Има много начини да се определи този квадратичен полином. По метода на Лагранж полинома се записва по следния начин:

(1)

Ще създадем функция interpf( ), която да изчислява точките на интерполиращата крива по метода на Лагранж (1), минаваща през три зададени точки. Реализацията на функцията може да се планира по следния начин:

• 
  Входни стойности: x=[x₁, x₂, x₃], y=[ y₁, y₂ , y₃], x_i
  
• 
  Изходни стойности: y_i
  
• 
  Задача: функцията да изчислява зависимоста y_i=f(x_i) използвайки полинома на Лагранж.
  

function yi=interpf(xi,x,y)

yi=y(1)*(xi-x(2))*(xi-x(3))/((x(1)-x(2))*(x(1)-x(3)))...

+y(2)*(xi-x(1))*(xi-x(3))/((x(2)-x(1))*(x(2)-x(3)))...

+y(3)*(xi-x(1))*(xi-x(2))/((x(3)-x(1))*(x(3)-x(2)));

return;

• 
  Потребителят задава точките (x₁, y₁), (x₂, y₂), (x₃, y₃), които трябва да се интерполират с квадратичен полином.
  
• 
  Определяме границите на интерполационната крива (от x_min до x_max).
  
• 
  Определяме стойностите на y_i за желаните интерполационни точки x_i като използваме функцията interpf.
  
• 
  Построяваме графика на получената интерполационна крива (x_i, y_i) като маркираме отделно трите задаващи точки.
  

clear all; help interp; %clear memory and print header

disp('Enter data points as x,y pairs(e.g.[1 2])');

for i=1:3

temp=input('Enter data point: ');

x(i)=temp(1);

y(i)=temp(2);

end

xr=input('Enter range of x values as [x_min x_max]:');

nplot=100; % number of points for interpolation curve

% Calculate curve points

for i=1:nplot

xi(i)=xr(1) + (xr(2)-xr(1))*(i-1)/(nplot-1);

yi(i)=interpf(xi(i), x, y);

end

plot(x,y,'*',xi,yi,'-');

xlabel('x');

ylabel('y');

Сега ще разгледаме по-общо същия проблем. Нека имаме набор от експериментални данни, които трябва да опишем с някакъв теоретично обоснован модел. Ако има несъответствие между експерименталните и теоретичните данни това се дължи или на неточно измерени данни т.е.експериментални грешки или на неподходящо избран т.е. неточен модел.

Ако моделът е линейна функция говорим за линейна регресия. Ако моделът е полином от втора степен казваме че данните се напасват с квадратична регресия и т.н. В MATLAB съществува библиотечна функция polyfit(x,y,d), която връща полином от степен d, който най-добре напасва данните дефинирани с векторите x и y. Както вече споменахме в MATLAB полиномите се представят от вектор с координати равни на коефициентите на полинома.

Ще разгледаме следния пример записан във файла fit1.m. Първо генерираме данни на полином от втора степен, към които добавяме случаен шум. След това използваме функцията polyfit за да се открие полинома, който най-добре моделира данните.

%this program interpolates random quadratic data

x=linspace(0,2,256);%generate x-samples

y=1 - x + 0.5*x.^2; %generate quadratic data points

plot(x,y,’.’);

hold on; %retain graph when next graph is added

yy=y+ 0.15*randn(size(x));%add random noise

plot(x,yy,'.k');

hold on;

yyy=polyval(p,x);%evaluate p for every x

plot(x,yyy,'-k');

hold off;

В следващия пример ще разгледаме Гаусов импулс с постоянна компонента. Моделът се описва от формулата:

(6)
Трябва да оптимизираме 4 параметъра, зададени с вектора р както следва: p=(y₀,s,a,x₀)

Първо ще дефинираме функцията peak, която изчислява несъответствието между набора от данни и модела при фиксирани параметри:

% this file is peak.m

function mf=peak(p,x,y);

yy=p(1)+p(2)*exp(-p(3)*(x-p(4)).^2);

mf=norm(yy-y);

Тук ако x е вектор x=[x₁,x₂,…x_N] Функцията norm(x) връща резултата:

%this file is fit2.m

x=linspace(0,5,1024);

p=[3;1;4;2.5];%this is the vector of the parameters

ys=p(1)+p(2)*exp(-p(3)*(x-p(4)).^2);%this is the signal

yn=ys+0.5*randn(size(x));%add noise to signal

pp=fminsearch('peak',[2.5;0.5;3;2],[],x,yn);%find vector of optimized parameters

yr=pp(1)+pp(2)*exp(-pp(3)*(x-pp(4)).^2);%model with optimized parameters

plot(x,ys,x,yn,'.',x,yr);

Тук е използвана библиотечната оптимизираща функция fminsearch('peak',[2.5;0.5;3;2],[],x,yn), която търси минимума на функцията peak(p,x,yn) при начални стойности на оптимизираните параметри 2.5, 0.5, 3, 2. Знакът [] указва да се вземат стойностите по подразбиране при определяне критериите за търсене и за спиране процеса на оптимизация. Върнатата стойност рр е вектора на оптимизираните параметри.