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