Podstawy Informatyki 3 – ćwiczenia10
Strona 1 z 8
Podstawy Informatyki
Programowanie w Matlab
Ć
wiczenia 10
Wybrane możliwości numeryczne pakietu Matlab – c.d.
Interpolacja
Zagadnienie interpolacji można sformułować następująco: na przedziale <a, b> danych jest
n+1
różnych punktów x0, x1, ..., xn, które nazywamy węzłami interpolacji, oraz wartości
pewnej funkcji y=f(x) w tych punktach f(x0)=y0, f(x1)=y1, ..., f(xn)=yn.
Zadaniem interpolacji jest wyznaczenie przybliżonych wartości funkcji w punktach nie
będących węzłami oraz oszacowanie błędu tych przybliżonych wartości.
W tym celu należy znaleźć funkcję F(x), zwaną funkcją interpolującą, która w węzłach
interpolacji przyjmuje takie same wartości co funkcja f(x).
Interpolacja jest w pewnym sensie zadaniem odwrotnym do tablicowania funkcji. Przy
tablicowaniu, mając analityczną postać funkcji budujemy tablicę wartości, przy interpolacji
natomiast, na podstawie tablicy wartości funkcji określamy jej postać analityczną.
Standardowe procedury Matlaba realizują interpolację za pomocą wielomianu pierwszego
stopnia (czyli funkcji kawałkami liniowej), wielomianu trzeciego stopnia oraz funkcji
sklejanych stopnia trzeciego.
Funkcje interpolujące
yi=interp1(x,y,xi,metoda)
Zwraca wektor yi, będący wartościami
funkcji jednej zmiennej y=f(x) w
punktach określonych wektorem xi;
węzły interpolacji określają wektory x i
y
; metoda to łańcuch znaków określający
metodę interpolującą
zi=interp2(x,y,z,xi,yi,metoda)
Zwraca macierz zi, zawierającą wartości
funkcji dwóch zmiennych z=f(x,y) w
punktach określonych wektorami xi i yi;
węzły interpolacji określają macierze
x
,y,z
vi=interp3(x,y,z,v,xi,yi,zi,metoda)
Interpolacja funkcją 3 zmiennych,
analogicznie jak interp2
vi=interpn(x1,x2,x3,...,v,y1,y2,y3,
...)
Interpolacja funkcją n zmiennych,
analogicznie jak interp2
Podstawy Informatyki 3 – ćwiczenia10
Strona 2 z 8
Metody interpolacji
‘linear’ – interpolacja liniowa,
‘spline’ – interpolacja funkcjami sklejanymi stopnia trzeciego,
‘cubic’ lub ‘pchip’ – interpolacja wielomianami stopnia trzeciego,
‘nearest’ – interpolacja wartością najbliższego sąsiada
Wszystkie metody interpolacji wymagają, aby ciąg x był monotoniczny.
Przykład 10.1
Napisz skrypt, w którym przetestujesz różne metody interpolacji funkcji jednej zmiennej.
Węzły interpolacji stanowią pary (x
i
, y
i
)
gdzie x
i
= 0, 1, 2, ..., 20, natomiast
y
i
= 2sin(x
i
)+sin(2x
i
)+sin(3x
i
+
π
)
%Przyklad10_1
%INTERPOLACJA skrypt testuje różne metody interpolacji 1-D
%
%----------interpolowana funkcja-----------------
f= inline(
'2.*sin(x)+sin(2.*x)+sin(3.*x +pi)'
);
%---------- węzły interpolacji-------------------
x=[0:20];
y=f(x);
%---------- metody interpolacji -----------------
x_int=[0:0.1:20];
y_int1=interp1(x,y,x_int,
'spline'
);
y_int2=interp1(x,y,x_int,
'pchip'
);
y_int3=interp1(x,y,x_int,
'nearest'
);
y_int4=interp1(x,y,x_int,
'linear'
);
%---------------- wykresy------------------------
subplot(4,1,1), plot(x,y,
'*'
,x_int,y_int1,
'-r'
)
hold
on
, fplot(f,[0,20])
title(
'Metoda funkcji sklejanych'
)
subplot(4,1,2),plot(x,y,
'*'
,x_int,y_int2,
'-r'
)
hold
on
, fplot(f,[0,20])
title(
'Metoda wielomianow trzeciego stopnia'
)
subplot(4,1,3),plot(x,y,
'*'
,x_int,y_int3,
'-r'
)
hold
on
, fplot(f,[0,20])
title(
'Metoda najblizszego sasiada'
)
subplot(4,1,4),plot(x,y,
'*'
,x_int,y_int4,
'-r'
)
hold
on
, fplot(f,[0,20])
title(
'Metoda funkcji liniowych'
)
Przykład 10.2
Przyjmując dane z tabeli znajdź wartości funkcji interpolującej w punktach -1; -0.8; -0.6 …. 3
x
-1
0
1
2
3
y
3
5
4
2
6
Podstawy Informatyki 3 – ćwiczenia10
Strona 3 z 8
%PRZYKLAD10_2
clear
all
% wyczyszczenie pamięci
%określenie węzłów interpolacji
x=[-1 0 1 2 3];
y=[3 5 4 2 6];
%określenie punktow xi w ktorych
% znajdujemy interpolowane wartości funkcji yi=f(xi);
xi=-1:0.2:3;
% interpolacja liniowa
yi_lin=interp1(x,y,xi,
'linear'
);
%interpolacja funkcjami sklejanymi
yi_spline=interp1(x,y,xi,
'spline'
);
% wykresy
plot(x,y,
'*'
, xi, yi_lin,
':'
,xi,yi_spline,
'--'
)
legend(
'wezly interpolacji'
,
'interpolacja liniowa'
,
'funkcja sklejana'
)
Interpolacja Lagrange’a
Interpolacja za pomocą wielomianu polega na znalezieniu wielomianu
)
( x
L
n
stopnia co
najwyżej n, by wartości tego wielomianu i funkcji interpolowanej w węzłach były sobie
równe czyli:
n
n
n
x
a
x
a
x
a
a
x
L
+
+
+
+
=
...
)
(
2
2
1
0
Zadanie interpolacji wielomianowej posiada jednoznaczne rozwiązanie, czyli istnieje tylko
jeden wielomian spełniający warunek. Wartości współczynników wielomianu wyliczane są ze
wzoru Lagrange’a:
Szukany wielomian ma postać:
)
)...(
)(
)...(
)(
(
)
)...(
)(
)...(
)(
(
)
(
1
1
1
0
1
1
1
0
0
n
j
j
j
j
j
j
j
n
j
j
n
j
j
n
x
x
x
x
x
x
x
x
x
x
x
x
x
x
x
x
x
x
x
x
y
x
L
−
−
−
−
−
−
−
−
−
−
=
+
−
+
−
=
∑
przyjmując, że:
)
)...(
)(
(
)
(
1
0
n
n
x
x
x
x
x
x
x
−
−
−
=
ω
wzór Lagrange’a ma wtedy postać:
∑
=
−
=
n
j
j
n
j
n
j
n
x
x
x
x
y
x
L
0
'
)
(
)
(
)
(
)
(
ω
ω
gdzie: y
j
= y(x
j
)
a ω
n
’
(x
j
)
jest wartością pochodnej wielomianu ω(x) w punkcie x
j
. Pozwala to
na napisanie procedury numerycznej do obliczania wartości wielomianu stopnia n w
dowolnym punkcie x leżącym w przedziale <a, b>.
Przykład 10.3
Napisz funkcję lagrange.m obliczającą wartość wielomianu interpolacyjnego Lagrange’a
stopnia n w dowolnym punkcie x (x może być wektorem). Argumentami funkcji powinny być
również n+1 wymiarowe wektory punktów węzłowych xj i yj.
Przetestuj napisaną funkcję interpolując punkty węzłowe leżące na krzywej będącej
wykresem funkcji y = 1/(1+25x
2
).
Przeprowadź obliczenia dla 3, 5 i 9 punktów leżących w
równych odległościach w przedziale <-1; 1>.
Podstawy Informatyki 3 – ćwiczenia10
Strona 4 z 8
%PRZYKLAD10_3 skrypt sprawdza interpolację wielomianami Lagrange'a
funk=inline(
'1./(1+25*x.^2)'
);
x=linspace(-1,1,201); y=funk(x);
ia=linspace(1,201,3);
ib=linspace(1,201,5);
ic=linspace(1,201,9);
xia=x(ia); yia=y(ia);
x1a=x; x1a(ia)=[];
ya=lagrange(xia,yia,x1a,2);
xib=x(ib); yib=y(ib);
xlb=x; xlb(ib)=[];
yb=lagrange(xib,yib,xlb,4);
xic=x(ic); yic=y(ic);
xlc=x; xlc(ic)=[];
yc=lagrange(xic,yic,xlc,8);
subplot(2,2,1)
fplot(funk,[-1,1])
ylim([-1.5 1.5]), title(
'funkcja interpolowana: 1/(1+25x^2)'
)
subplot(2,2,2)
fplot(funk,[-1,1]), hold
on
, plot(xia, yia,
'o'
), plot(x1a,ya,
'--'
)
ylim([-1.5 1.5]), title(
'funkcja interpolujaca: wielomian rzedu 2'
)
subplot(2,2,3)
fplot(funk,[-1,1]), hold
on
, plot(xib, yib,
'o'
), plot(xlb,yb,
'--'
)
ylim([-1.5 1.5]), title(
'funkcja interpolujaca: wielomian rzedu 4'
)
subplot(2,2,4)
fplot(funk,[-1,1]), hold
on
, plot(xic, yic,
'o'
), plot(xlc,yc,
'--'
)
ylim([-1.5 1.5]), title(
'funkcja interpolujaca: wielomian rzedu 8'
)
gdzie funkcja lagrange.m ma postać:
function
[y]=lagrange(xj,yj,x,n)
%LAGRANGE wzór interpolacyjny Lagrange'a
%LAGRANGE(xj,yj,x,n) oblicza wartości wielomianu interpolacyjnego stopnia n
%w dowolnym punkcie x dla węzłów xj i odpowiadających im wartości funkcji
%interpolowanej yj.
%np.: y=lagrange([1 2 3], [1 2 3],1.7,2) daje y=1.7
omega_w=1;s=0;
for
j=1:n+1
omega_w=omega_w.*(x-xj(j));
omega_p=1;
for
i=1:n+1
if
i~=j
omega_p=omega_p.*(xj(j)-xj(i));
end
end
s=s+yj(j)./(omega_p.*(x-xj(j)));
end
y=omega_w.*s;
Podstawy Informatyki 3 – ćwiczenia10
Strona 5 z 8
Aproksymacja
Aproksymacja jest innym zagadnieniem niż interpolacja, chociaż sformułowanie jest
podobne. Załóżmy, że w przedziale <a, b> danych jest n+1 punktów x
0
= a, x
1
, x
2
, ..., x
n
= b
oraz wartości pewnej funkcji y=f(x) w tych punktach. Poszukiwana jest funkcja F(x) dobrze
przybliżająca f(x) w przedziale <a, b>, z tym, że, w odróżnieniu od interpolacji, nie żąda się,
by w podanych punktach obydwie funkcje miały takie same wartości. Pojawia się zatem błąd
aproksymacji zdefiniowany jako różnica pomiędzy wartościami funkcji aproksymującej F i
funkcji aproksymowanej f w punktach x
i
, i=0, 1,..., n. Zadaniem metody numerycznej jest w
tym przypadku minimalizacja błędu aproksymacji. Aproksymacja jest nazywana również
dopasowywaniem krzywej do danego zbioru danych. Jeżeli punkty (x
i
, y
i
) pochodzą z badań
eksperymentalnych nad procesami fizycznymi, to zwykle ich wartości są obarczone błędem
pomiaru. Interpolacja, czyli wymóg znalezienia funkcji przechodzącej dokładnie przez te
punkty, w tym przypadku nie ma sensu.
Jedną z najpopularniejszych metod aproksymacji jest aproksymacja średniokwadratowa,
zwana także metodą najmniejszych kwadratów. Polega ona na minimalizacji błędu
aproksymacji zdefiniowanego jako suma kwadratów odchyłek we wszystkich punktach
obliczeniowych:
[
]
2
0
)
(
)
(
)
(
i
i
n
i
i
x
f
x
F
x
w
−
∑
=
Do poprawy przybliżenia w wybranych węzłach może być wykorzystana funkcja wagowa
w(x
i
)
. Zwykle jednak jest ona tożsamościowo równa jedności, zatem wszystkie odchyłki są
jednakowo istotne.
Podstawową funkcją jaką oferuje Matlab do rozwiązania zadanie aproksymacji jest funkcja
polyfit
.
a = polyfit(x, y, n)
Znajduje współczynniki wielomianu n-tego
stopnia, który w sensie najmniejszych
kwadratów najlepiej pasuje do serii danych
pomiarowych (x, y)
Uwaga: do obliczania wartości wielomianu:
W(x,a)=a
1
x
n
+a
2
x
n-1
+...+a
n
x+a
n+1
,
o współczynnikach a w punkcie x można wykorzystać funkcję polyval(a, x), gdzie
a(1)
jest współczynnikiem stojącym przy x
n
, a(2)
jest współczynnikiem stojącym przy x
n-1
,
...., a(n+1) jest wyrazem wolnym.
Przykład 10.4
Przypuśćmy, że mamy następujące dane x
i
oraz y
i
i chcielibyśmy znaleźć najlepsze liniowe
dopasowanie do tych danych.
x
i
2
4
7
12
20
y
i
3
5
15
22
42
%PRZYKLAD10_4 skrypt dopasowuje krzywą do danych pomiarowych
%
x=[2 4 7 12 20];
y=[3 5 15 22 42];
figure, plot(x,y,
'o'
);
a=polyfit(x,y,1);
Podstawy Informatyki 3 – ćwiczenia10
Strona 6 z 8
xt=0:25;
yt=polyval(a,xt);
hold
on
, plot(xt,yt,
'r-'
)
legend(
'dane pomiarowe'
,
'dopasowanie'
);
Poza funkcją polyfit Matlab daje użytkownikowi narzędzie Basic Fitting,
uruchamiane z menu Tools w oknie graficznym, pozwalające na: automatyczne dopasowanie
krzywych sklejanych (spline) oraz wielomianowych (do dziesiątego stopnia), wyświetlenie
równania dopasowanej krzywej, obliczenie odchyleń w punktach danych oraz wartości błędu
aproksymacji (pierwiastek z sumy kwadratów odchyleń). Narzędzie to pozwala w łatwy
sposób porównać różne dopasowania i wybrać najlepsze z nich.
Przykład 10.5
Załóżmy, że obserwujemy pewien proces w czasie t∈<0, 10> , który powinien mieć charakter
sinusoildalny. Pomiary zawierają jednak pewien błąd losowy. Za pomocą narzędzia Basic
Fitting
dopasować najlepszą krzywą do tych danych.
%PRZYKLAD10_5 wykorzystanie narzędzia Basic Fitting
%
t=linspace(0,10,10);
%momenty pomiaru 0,1,2,...,10
y=sin(0.5*t)+0.1*randn(size(t));
%wygenerowane dane pomiar. z błędem losowym
plot(x,y,
'o'
);
%wykres surowych danych pomiarowych
Po narysowaniu surowych danych przejdź do okna graficznego i z menu Tools wybierz
polecenie Basic Fitting. Zaznacz kolejno opcje quadratic, cubic oraz 4th degree polynomial, a
następnie zaznacz opcje: Show equations, Plot residuals oraz Show norm of residuals.
Wynik tych poleceń pokazuje poniższy rysunek.
Jak widać, spośród tych trzech krzywych najlepsze dopasowanie, mierzone sumą kwadratów
odchyłek, wykazuje wielomian czwartego stopnia o podanym równaniu.
Podstawy Informatyki 3 – ćwiczenia10
Strona 7 z 8
Korzystając z funkcji polyfit, możemy dopasowywać do danych również inne krzywe, np.
wykładnicze:
y = ae
bx
czy potęgowe
y = qx
p
Wymaga to jednak wstępnego przekształcenia równań oraz danych. Najczęściej należy
przekształcić dane do skali logarytmicznej.
Zauważmy, że logarytmując stronami powyższe równania otrzymamy:
- dla pierwszego równania: ln(y) = ln(a)+bx czyli dostaniemy przekształconą krzywą:
y’ = a
0
+a
1
x
, gdzie y’ = ln(y), a
0
= ln(a)
oraz a
1
= b
.
- dla drugiego równania: ln(y) = ln(q)+pln(x) czyli dostaniemy przekształconą krzywą:
y’ = a
0
+a
1
x’
, gdzie y’ = ln(y), a
0
= ln(q)
, a
1
= p
oraz x’=ln(x)
Możemy zatem w obu przypadkach zastosować, do tak zmienionych danych, aproksymację
liniową za pomocą funkcji polyfit.
Przykład 10.6
W tabeli przedstawiono dane dotyczące wahań ciśnienia w czasie, odczytane z pompy
podciśnieniowej. Dopasuj krzywą P(t) = P
0
e
-t/T
i znajdź stałe P
0
oraz T.
t
i
0
0.5
1.0
5.0
10.0
20.0
P
i
760
625
528
85
14
0.16
Obliczając logarytm obu stron równania otrzymujemy:
ln(P) = ln(P
0
) – t/T
czyli P’ = a
1
t + a
0
, gdzie P’=ln(P), a
1
= -1/T
, oraz a
0
= ln(P
0
).
Po wyznaczeniu zatem współczynników a
0
i a
1
za pomocą funkcji polyfit, będziemy
mogli z łatwością obliczyć stałe P
0
oraz T.
%PRZYKLAD10_6 dopasowanie krzywej wykładniczej
t=[0 0.5 1 5 10 20];
P=[760 625 528 85 14 0.16];
% pierwotne dane
Pp=log(P);
% dane przekształcone
a=polyfit(t,Pp,1);
% dopasowanie liniowe do przekształconych danych
P0= exp(a(2));
%obliczenie stałej P0
T = -1/a(1);
%obliczenie stałej T
disp([
'stała P0 = '
,num2str(P0),
', stała T = '
,num2str(T)])
%narysowanie wykresu w skali liniowej
twyk=linspace(0,20,20);
% dane do wykresu
Pwyk=P0*exp(-twyk/T);
figure, subplot(2,1,1)
plot(t,P,
'o'
,twyk, Pwyk,
'r-'
), grid
on
xlabel(
'czas t (s)'
), ylabel(
'Cisnienie'
)
title(
'Wykres cisnienia w pompie'
)
%narysowanie wykresu w skali półlogarytmiczenj
Pwyk_log= exp(polyval(a,twyk));
%dane do wykresu półlogarytmicznego
subplot(2,1,2)
semilogy(t,P,
'o'
,twyk, Pwyk_log,
'r-'
), grid
on
xlabel(
'czas t (s)'
), ylabel(
'Cisnienie'
)
title(
'Wykres semilogarytmiczny cisnienia w pompie'
)
Podstawy Informatyki 3 – ćwiczenia10
Strona 8 z 8
0
2
4
6
8
10
12
14
16
18
20
0
200
400
600
800
czas t (s)
C
is
n
ie
n
ie
Wykres cisnienia w pompie
0
2
4
6
8
10
12
14
16
18
20
10
-1
10
0
10
1
10
2
10
3
czas t (s)
C
is
n
ie
n
ie
Wykres semilogarytmiczny cisnienia w pompie
Zadania do ćwiczeń 10
10.1 Poniższa tabela przedstawia dane pomiaru wilgotności z roku 1962 dla stacji
monitoringowej Kasprowy Wierch. Znajdź ciągłą funkcję zmian wilgotności wykorzystując
znane metody interpolacji.
Miesiąc
I
II III IV V VI VII VIII IX X XI XII
wilgotność
75 90 89 82 92 90 88
79
78 68 91 77
10.2 Wykorzystując narzędzie Basic Fitting znajdź najlepiej dopasowaną krzywą do
danych pomiarowych dotyczących stężenia ozonu w atmosferze pochodzących z pewnej
stacji meteorologicznej. Badania prowadzone były przez miesiąc co godzinę, a następnie
obliczono średnie stężenie ozonu dla jednego dnia.
godz.
0
1
2
3
4
5
6
7
8
9
10
11
wartość
67
61
50
35
33
38
70
90
105 121 130 130
godz.
12
13
14
15
16
17
18
19
20
21
22
23
wartość
129 117 101
91
100 105
95
83
66
74
78
60
10.3 Wykorzystując funkcję polyfit znajdź wielomian drugiego stopnia, który najlepiej
aproksymuje funkcję y=sin(x) w przedziale <0, π>. Narysuj wykresy funkcji y=sin(x) oraz
aproksymującego ją wielomianu.
10.4 Znajdź parametry modelu Malthusa wzrostu populacji N(t)=N
0
e
rt
wiedząc, że w
chwilach t = 0, 2, 4, 6, 8, 10 wielkość populacji wynosiła odpowiednio: 100, 150, 220,
330, 500, 740