Metody numeryczne - laboratorium.

Matlab 6.1

wersja 2.3

Kwiecień 2002

Spis treści

Matlab - podstawy

1.1

Wiadomości ogólne.

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1.2

Praca z Matlab -em . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1.3

Wyrażenia w Matlab -ie. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1.3.1

Zmienne.

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1.3.2

Liczby. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1.3.3

Operatory.

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1.3.4

Funkcje. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1.4

Pliki skryptowe i funkcjne. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1.5

Rysunki w Matlab -ie - polecenia plot i subplot. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1.6

Pliki danych . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1.6.1

Polecenie save. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1.6.2

Polecenie load. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1.7

Literatura . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Nieuniknione błędy.

2.1

Miary błędów.

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

2.2

Reprezentacja binarna liczby rzeczywistej. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

2.3

Zamieniamy całkowitą potęgę do układu binarnego . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

2.4

Zamieniamy rzeczywistą mantysę do układu binarnego . . . . . . . . . . . . . . . . . . . . . . . . . . . .

2.5

Liczby użyteczne i nieużyteczne.

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Macierze i operacje na nich, techniki przetwarzania macierzy.

3.1

Mnożenie macierzy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3.2

Potęgowanie i dzielenie macierzy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3.3

Transpozycja i sprzężenie macierzy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3.4

Rząd macierzy

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3.5

Wyznacznik macierzy

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Macierze specjalne . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3.6.1

Macierz diagonalna (przekątniowa) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3.6.2

Macierz jednostkowa . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3.6.3

Macierz permutacji . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3.6.4

Macierz trójkątna

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3.6.5

Macierz odwrotna

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3.6.6

Macierz symetryczna . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3.6.7

Macierz ortogonalna . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Układy równań liniowych.

4.1

Rozwiązanie formalne dla kwadratowej macierzy układu. . . . . . . . . . . . . . . . . . . . . . . . . . . .

4.2

Eliminacja Gaussa. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4.2.1

Rozwiązanie układu diagonalnego . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4.2.2

Rozwiązanie układu trójkątnego

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4.2.3

Eliminacja Gaussa . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4.2.4

Eliminacja Gaussa z wyszukiwaniem elementu głównego . . . . . . . . . . . . . . . . . . . . . . .

4.2.5

Rozwiązanie układu z użyciem opratora dzielenia lewostronnego. . . . . . . . . . . . . . . . . . .

4.3

Co robi operator dzielenia lewostronnego? . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Aproksymacja.

5.1

Aproksymacja linią prostą . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

5.2

Aproksymacja . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

5.3

Co jeszcze robi operator dzielenia lewostronnego? . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Interpolacja.

6.1

Interpolacja wielomianami sklejanymi. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

6.2

Interpolacja kawałkami liniowa. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

100

6.3

Interpolacja kawałkami sześcienna - Hermite. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

102

6.4

Interpolacja kawałkami sześcienna - splajny. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

103

6.5

Funkcje interpolacyjne Matlab -a. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104

Wstęp

1. W obliczniach numerycznych (w przeciwieństwie do obliczeń symbolicznych)

dopuszczamy błąd w reprezentacji liczby.

2. W celu zaoszczędzenia pamięci komputera stosuje się różne typy zmiennych.

3. Dla potrzeb nauki praktycznego zapoznania się z metodami numerycznymi

będziemy stosować program Matlab (MATrix LABoratory), który wykorzy-
stuje tylko jeden typ zmiennych - macierz.












. . . a

...

. . . a












Matlab - podstawy

1.1.

Wiadomości ogólne.

Matlab umożliwia wykorzystanie metod rachunku macierzowego za pomocą in-
teraktywnego interfejsu i prostego języka poleceń. Wszystkie polecenia są formu-
łowane w trybie tekstowym, a podstawową strukturą danych jest macierz, której
elementami mogą być liczby rzeczywiste lub zespolone.

W skład interfejsu wchodzą między innymi:

•

Command Window - okno poleceń

•

Command History - okno poprzednich poleceń

•

Launch Pad - okno dostępu do pakietów Matlab -a

•

Current Directory - okno katalogu roboczego

•

Workspace - inoformacje o zmiennych w przestrzeni roboczej

1.2.

Praca z Matlab -em

Praca z Matlab -em przypomina pracę w typowym systemie operacyjnym (DOS,
UNIX) i polega na wydawaniu poleceń, które po zatwierdzeniu są wykonywane
przez interpreter.

•

polecenie może być jedno, bądź jest to ciąg poleceń oddzielonych przecinkiem
lub średnikiem.

•

pliki, które są używane w pakiecie Matlab mają rozszerzenie *.m

•

w Matlab -ie są rozróżnialne małe i duże litery!

Podstawowe polecenia Matlab -a:

Enter

wykonanaj polecenie

↑

ponowne wyświetlenie polecenia na ekranie

help

pomoc na dowolny temat

help nazwa tematu pomoc na określony temat
exit

zakończenie pracy z Matlab -em (lub quit)

demo

pokaz przykładowych zastosowań Matlab -a

Polecenia DOS-a:

dir

wyświetlenie zawartości aktualnego katalogu

type nazwa pliku wyświetlenie zawartości pliku
!notepad

wywołanie ”Notatnika” (lub innego programu)

W Matlab -ie podstawowym poleceniem jest instrukcja przypisania:

zmienna = wartość

>> x=5

x =

>> z=6;

jeśli na końcu polecenia zostanie posta-
wiony średnik, to wartość zmiennej nie
będzie wyświetlona,

zmienna = wyrażenie

>> 5+4

ans =

>> z=5*7
z =

jeśli nie zostanie zdefiniowana zmienna
(nie nadana nazwa zmiennej), to Ma-
tlab przypisuje jej nazwę zmiennej ro-
boczej ans,

1.3.

Wyrażenia w Matlab -ie.

Matlab udostępnia użytkownikowi polecenia - komendy, które w odróżnieniu od
innych języków programowania związane są z macierzami. Polecenia wykonują
operacje na:

•

zmiennych,

•

liczbach,

•

operatorach,

•

funkcjach.

1.3.1.

Zmienne.

Matlab nie wymaga deklaracji typu i wymiarów zmiennych. Automatycznie two-
rzy zmienną i przydziela jej domyślny obszar pamięci. Jeżeli zmienna już istnieje,

Matlab zmienia jej zawartość i w sposób dynamiczny przyporządkowuje jej wy-
magany obszar pamięci. Polecenie:

>> num_students = 25

tworzy

macierz

1 × 1

nazwie

num_students

i przechowuje w niej

wartość 25.

Nazwa zmiennej musi zaczynać się od litery alfabetu, a po niej może występować
dowolny znak alfanumeryczny lub znak podkreślenia. Maksymalna długość nazwy
zmiennej wynosi 31 znaków. Przy wyborze nazw zmiennych należy pamiętać o
tym, że Matlab rozróżnia małe i duże litery.

1.3.2.

Liczby.

Matlab używa konwencjonalną notację dziesiętną, z kropką dziesiętną. W przy-
padku notacji naukowej litera e służy do określenia wykładnika potęgi dziesięć. W
przypadku liczb zespolonych używa się zarówno literę i, jak j przy części urojonej
liczby zespolonej. Przykłady liczb w programie Matlab :

-99

0.0001

9.6397238

1.60210e-20

6.02252e23

-3.14159j

3e5i

Wszystkie liczby są zachowywane w formacie długim (z ang. format long). Na
komputerach wykorzystujących 32-bitowy procesor odpowiada to dokładności do
16 cyfr znaczących po kropce dziesiętnej i zakresowi od 10

−308

do 10

+308

Polecenie format kontroluje format liczb wyświetlanych przez program MATLAB.
Polecenie jest jedynie odpowiedzialne za sposób wyświetlania liczb, a nie dokład-
ność obliczeń i sposób ich zapamiętywania.
Ćwiczenie:

Format liczby

>> a = 0.12345678901234567890

polecenie nadania wartości zmiennej a

>> a

wyświetlenie wartości zmiennej a

>> format short e

>> a

ustawienie formatu liczb na 0.12345 10

123

>> format short

>> a

ustawienie formatu liczb na 0.12345

>> format long
>> a

format 0.123456789012345

>> format long e
>> a

format 0.123456789012345 10

123

1.3.3.

Operatory.

Operatory arytmetyczne (macierzowe i tablicowe)

Symbol operacji Znaczenie

Symbol operacji Priorytet

macierzowej

tablicowej

operacji

potęgowanie

najwyższy

’

sprzężenie

nie ma

najwyższy

.’

transpozycja

nie ma

najwyższy

mnożenie

niższy

dzielnie prawostronne

niższy

dzielenie lewostronne

niższy

dodawanie

najniższy

odejmowanie

najniższy

Operatory logczne

Symbol Znaczenie

mniej niż

mniej niż lub równe

większe niż

większe niż lub równe

równa się

nie równa się

iloczyn logiczny (i)

suma logiczna (lub)

zaprzeczenie (nie)

Uwaga! W języku Matlab tylko liczba zero ma wartość logiczną false - fałsz,
pozostałe liczby mają wartość logiczną true - prawda.

1.3.4.

Funkcje.

Matlab udostępnia użytkownikowi szereg standardowych elementarnych funk-
cji obejmujących m.in. abs, sqrt, exp i sin. Większość tych funkcji akceptuje
zespolone argumenty.

>> help elfun

polecenie wyświetlające listę elementar-
nych funkcji matematycznych

>> help specfun
>> help elmat

polecenie wyświetlające listę zaawanso-
wanych matematycznych i macierzowych
funkcji

Niektóre funkcje takie, jak sqrt i sin, są wbudowane w program, dzięki czemu
są bardzo efektywne. Inne funkcje takie, jak gamma i sinh, są zewnętrznymi pli-
kami *.m. Użytkownik może przejrzeć kody tych funkcji i ewentualnie poddać je
modyfikacjom.

Matlab udostępnia również szereg funkcji specjalnych generujących pomocne
stałe:

3.14159265358979

liczba urojona

√

−1

liczba urojona (to samo znaczenie co i)

eps

odległość od 1.0 do następnej liczby

realmin najmniejsza dodatnia liczba
realmax największa ujemna liczba
Inf

nieskończoność (∞)

NaN

nie jest to liczba (z ang. not a number)

>> eps
ans =

2.220446049250313e-016

>>eps = 1.e-6

nazwy powyższych funkcji nie są zastrze-
żone, dzięki czemu jest możliwa zmiana
ich wartości

>>clear eps

domyślna wartość funkcji może być od-
tworzona przy pomocy polecenia

Przykład:

>> a=10/0

Warning: Divide by zero.
a =

Inf

>>b=0/0
Warning: Divide by zero.
b =

NaN

Ostrzeżenie: Dzielenie przez zero.

>>c=Inf-Inf

c =

NaN

Przykład:

Wyrażenia z wykorzystaniem wbudowanych funkcji.

>>rho = (1+sqrt(5))/2
rho =
1.6180

>>a = abs(3+4i)
a =
5

>>duzo = exp(log(realmax))
duzo =
1.7977e+308

>>zbytduzo = pi*duzo
zbytduzo =
Inf

funkcje sqrt() -

√

·, abs() - |·|, . . .

Ćwiczenie:

Wypełnij macierze. Sprawdź rezutalt poleceniem whos.

>> a=1.234;
>> size(a)

polecenie size zwraca wymiar macierzy,
zmienna a to macierz o wymiare 1x1

>> clear a

polecenie usunięcia zmiennej a z pamięci

>> a=[1 2 3 4 5; 6 7 8 9 0]

wypełnianie macierzy (; jest separato-
rem wierszy)

>> a=[1.1:2.2:7]

wypełnianie macierzy (min:skok::max),
domyślny skok to 1

>> who

polecenie who pokazuje wszystkie zmien-
ne w przestrzeni roboczej

>> whos

polecenie

whos

pokazuje

wszystkie

zmienne i ich rozmiary

Ćwiczenie:

Wykonaj polecenia - spróbuj przewidzieć wynik. Sprawdź rezutalt

poleceniem whos.

>> a=[1,0,0,0,0,1]

zapamiętaj rolę przecinka - oddziela on
koleje pola macierzy

>> b=[2;4;6;10]

zapamiętaj rolę średnika - to separator
wierszy macierzy

>> c=[5 3 5; 6 2 -3]

spacja pełni rolę taką jak przecinek

>> d=[3 4

5 7
9 10]

zapamiętaj rolę zmiany linii (Enter) w
trakcie wypełniania macierzy, jest ona ta-
ka sama jak rola średnika

>> e=[3 5 10 0; 0 0 ...
0 3; 3 9 9 8]

... to znak kontynuacji wiersza

Ćwiczenie:

Wykonaj polecenia - spróbuj przewidzieć wynik. Sprawdź rezutalt

poleceniem whos.

>> t=[4 24 9]

>> q=[t 0 t]

zauważ, że elementami macierzy mogą
być macierze

>> x=[3 6]

>> y=[d;x]
>> z=[x;d]

macierz d była wprowadzone wcześniej i
zapamiętana w przestrzeni roboczej

>> r=[c;x,5]

>> c

nie pomyl średnika z przecinkiem i ko-
niecznie sprawdź poleceniem whos co już
masz w pamięci, macierz c zdefinowano
wcześniej

>> c(2,1)
>> v=[c(2,1);b]

c(2,1) wskazuje konkretny element ma-
cierzy

>> a(2,1)=-3
>> c(1,5)=1

czy wynik tych poleceń jest zaskakujący
?

1.4.

Pliki skryptowe i funkcjne.

Pliki *.M są plikam ASCII, które mogą zawierać:

•

sekwencje poleceń,

•

wywołania innych plików M lub samych siebie,

•

skrypty, czyli ciągi poleceń - pliki-M skryptowe,

•

funkcje tworzone przez użytkownika - pliki-M funkcyjne.

Pliki-M funkcyjne zawierają definicje nowych funkcji, które działają na zmiennych
lokalnych i globalnych. Komunikują się z przestrzenią roboczą poprzez parametry
i zmienne globalne (polecenie global).
Pliki-M funkcyjne rozpoczynają się od słowa kluczowego function o następującej
składni:

function [lista param. wyjść] = nazwa funkcji (lista param. wejść)

Ćwiczenie:

1. Utwórz za pomocą dowolnego edytora ASCII plik o nazwie pit.m

2. Wprowadź do niego następujące polecenia:

function z = pit(x,y);
% moja funkcja
%
temp = x.^2 + y.^2;
z = sqrt(temp);

3. Zapisz plik w katalogu roboczym.

4. Wykonaj polecenia z wykorzystaniem nowej funkcji pit.

5. Spróbuj wykonać inną własną funkcję.

>> pit(3,4)
>> a=4, b=6, c=pit(a,b)

Wywołaj funkcję na różne znane ci sposo-
by. Wykonaj polecenie help pit. Zmo-
dyfikuj komentarz (tekst poprzedzony
znakiem %)

1.5.

Rysunki w Matlab -ie - polecenia plot i subplot.

Polecenie plot służy do sporządzania rysunków 2-D i może być wywołane z na-
stępującymi parametrami:

plot(y)
plot(x,y)
plot(x,y,s)
plot(x,y,s,x1,y1,s1,...)

y - wektor wrtości na osi pionowej, na osi

poziomej będą kolejne liczby natu-
ralne

x - wektor wartości na osi poziomej o

tym samym wymiarze co wektor y

s - ciąg co najwyżej trzech znaków okre-

ślających: kolor, rodzaj punktu, typ
lini.

niebieski

punkt

ciągła

zielony

koło

kropkowa

czerwony

krzyż

kreska-kropka

turkusowy

plus

przerywana

amarantowy

gwiazda

żółty

kwadrat

czarny

diament

v^<>

trójkąt

Przykład:

>> x=-pi:0.1:pi;
>> plot(x,sin(x),’+-’,x,cos(x),’*’)

pierwsze polecenie wypełnia macierz x w
zakresie od −π do π co 0.1, drugie rysuje
dwie funkcje
polecenie sin(x) generuje macierz war-
tości dla wszystkich warości macierzy x

−4

−3

−2

−1

−0.8

−0.6

−0.4

−0.2

0.2

0.4

0.6

0.8

Polecenie subplot dzieli okno graficzne na mniejsze prostokątne pola.

subplot(m,n,p)
subplot(mnp)

m,n - ilość pól pionowo i poziomo

p - numer aktywnego pola

Przykład:

>> x = -5:0.01:5;
>> subplot(2,2,2), plot(x,sin(x))
>> subplot(2,1,2), plot(x,cos(x))

−5

−4

−3

−2

−1

−0.5

0.5

−5

−1

−0.5

0.5

Polecenia dodatkowe do opisu rysunku:

polecenie opis

przykład

title

tytuł

title(’Moj rysunek’)

xlabel

opis osi poziomej

xlabel(’x w cm’)

ylabel

opis osi pionowej

ylabel(’f(x)’)

legend

utworzenie legendy

legend(’moja linia’)

text

wstawienie tekstu w punkcie (x, y)

text(1,1,’moj tekst’)

gtext

wstawienie tekstu przy pomocy myszy gtext(’tekst’)

grid

wyświetlanie siatki

axis

skalowanie i wygląd osi

axis off, axis on

Inne polecenia związane z rysunkiem:

hold

nowy rysunek nie usunie istniejącego

wydruk lub zapis rysunku do pliku

clf

wyczyszczenie okna graficznego

Przykład:

Rysunek 3-D.

>> [x,y] = meshgrid(-2:.2:2, -2:.2:2);
>> z = x .* exp(-x.^2 - y.^2);
>> mesh(z)

−0.5

0.5

1.6.

Pliki danych

1.6.1.

Polecenie save.

Polecenie save pozwala na zapamiętanie zawartości pamięci operacyjnej progra-
mu w pliku danych z rozszerzeniem *.MAT, które następnie można wczytać do
programu dzięki poleceniu load.

save

save ’nazwa pliku’
save ’nazwa pliku’ ’lista zmiennych’

’nazwa pliku’ - nazwa pliku do które-

go zapisżą się dane (domyślnie ’ma-
tlab.mat’),

’lista zmiennych’ - lista zmiennych do

zapisania (domyślnie wszystkie).

Dodatkowe opcje:

-ascii

w 8-cyfrowym formacie tekstowym

-double w 16-cyfrowym formacie tekstowym
-tabs

elementy macierzy rozdzielone zostają znakami tabulacji

-append dopisuje dane do istniejącego pliku *.MAT

Uwaga! W jednym pliku tekstowym (ASCII) powinna być zapamiętana jedna
zmienna.

Przykład:

>> clear

>> A=[1,2,3;4,5,6;7,8,9];
>> b=12;
>> save
Saving to: matlab.mat

zapis zmiennych A i b do pliku ’ma-
tlab.mat’

>> clear

>> who
>> load
Loading from: matlab.mat
>> who
Your variables are:
A

czyścimy wszystkie zmienne poleceniem
clear, upewniamy się, że to prawda po-
leceniem who i czytamy zapisane wcze-
śniej zmienne z pliku ’matlab.mat’

>> save zmiennaa A

>> clear
>> load zmiennaa
>> who
Your variables are:
A

zapis i odczyt zmiennej A z pliku ’zmien-
naa.mat’

Przykład:

>> A

A =
1

>> save zmiennaa.dat A -ascii -double -tabs
>> type zmiennaa.dat
1.0000000000000000e+000
2.0000000000000000e+000
4.0000000000000000e+000
5.0000000000000000e+000

>> whos

Name

Size

Bytes

Class

2x2

double array

Grand total is 4 elements using 32 bytes

Przykład:

>> B=A’

B =
1

>> clear A
>> whos

Name

Size

Bytes

Class

3x3

double array

Grand total is 9 elements using 72 bytes

>> load zmiennaa
>> who
Your variables are:
A

Przykład:

>> save zmienne A B

>> clear
>> load zmienne B
>> who
Your variables are:
B

1.6.2.

Polecenie load.

Polecenie load wczytuje pliki binarne oraz pliki tekstowe zawierające dane nume-
ryczne. Plik tekstowy powinien zawierać prostokątną tablicę liczb rozdzielonych
spacjami. W pliku tym linia zawiera jeden wiersz macierzy. Każdy z wierszy po-
siada równą liczbę elementów.

Ćwiczenie:

Utwórz plik tekstowy magik.dat zawierający następujące cztery

linie:

16.0

3.0

2.0

13.0

5.0

10.0

11.0

8.0

9.0

6.0

7.0

12.0

4.0

15.0

14.0

1.0

>> load magik.dat
>> magik

wczytuje plik danych i tworzy zmienna
magik zawierającą przykładową macierz.

1.7.

Literatura

1. Matlab 6 - poradnik użytkownika, Bogumił i Zbigniew Mrozek, ISBN 83-

7101-449-X, cena 38zł, 230 stron.

2. Matlab - obliczenia numeryczne i ich zastosowania, Andrzej Zalewski i Rafał

Cegieła, ISBN 83-85060-85-5, cena 30zł, 400 stron.

Nieuniknione błędy.

Przykład:

>> format long e

Polecenie zmiany formatu wyświetlania
liczby (help format).

>> 2.6 + 0.2
ans =

2.800000000000000e+000

>> ans + 0.2
ans =

3.000000000000000e+000

>> ans + 0.2
ans =

3.200000000000001e+000

>> 2.6 + 0.6
ans =

3.200000000000000e+000

Dlaczego 2.6 + 0.2 + 0.2 + 0.2 6= 2.6 + 0.6?

2.1.

Miary błędów.

błąd bezwzględny = |wartość otrzymana − wartość prawdziwa|

błąd względny =

wartość otrzymana − wartość prawdziwa

wartość prawdziwa

2.2.

Reprezentacja binarna liczby rzeczywistej.

•

Każdą liczbę można zamienić do następującej postaci, nazywanej znormalizo-
waną notacją naukową:

123.456 → 0.123456

}

mantysa

×10

•

Skończona liczba bitów jest przeznaczana na zapamiętanie liczby, przy czym
część bitów b

przeznaczana jest na znormalizowana mantysę z przedziału

(0.1,1.0), a część na całkowitą potęgę b

. Zarówno b

jak i b

mogą przyjmować

wartości 1 lub 0.

64 bity

{

|{z}

znak

· · · · · · · · · b

}

52 bity

}

mantysa

|{z}

znak

· · · b

}

10 bitów

}

potęga

2.3.

Zamieniamy całkowitą potęgę do układu binarnego

= p

−

k=10

p
k

k−1

= p

k−1

− b

p
k

k−1

gdzie p

= potęga

Przykład:

Niech potęga p

= 10.

k 2

k−1

p
k

= p

k−1

− b

p
k

k−1

...

4 8

3 4

2 2

1 1

2.4.

Zamieniamy rzeczywistą mantysę do układu binarnego

= m

−

k=1

m
k









= m

k−1

− b

m
k









gdzie m

= mantysa

Przykład:

Niech mantysa m

= 0.8125.

1
2

= m

k−1

− b

1
2

1 0.5

0.3125

2 0.25

0.0625

3 0.125

0.0625

4 0.0625

0.0000

...

Ćwiczenie:

1. Rozpakuj plik mn.zip w katalogu: ...\toolbox\

2. Ustaw

bieżący

katalog

(okno

Current

Directory)

na:

...\toolbox\mn\bledy.

3. Makra Matlab -a mają rozszerzenie *.m. Sprawdź, czy widoczne jest makro

man2bin.m.

4. Kliknij dwa razy w nazwę pliku man2bin.m. Zapoznaj się z jego zawartością:

(a) pierwsza linia zawiera zawsze nagłówek funkcji o takiej samej nazwie jak

nazwa pliku,

(b) linie poprzedzone znakiem % to komentarze,

5. Opis funkcji wywołuje się przez polecenie help "nazwa_funkcji" w oknie

poleceń (Command Window).

2.5.

Liczby użyteczne i nieużyteczne.

>> relmax

>> relmin

Sprawdź wartości tych zmiennych. Da-
lej ich wartości będziemy oznaczać odpo-
wiednio: R

max

i R

min

Zakresy dyskretnej reprezentacji liczb rzeczywistych:

przydatność liczb

zakres wartości

komentarz

nieużyteczne

(−∞, −R

max

)

przepełnienie (z ang. overflow)

użyteczne

(−R

max

, −R

min

)

pełna precyzja mantysy

część. użyteczne

(−R

min

, −R

min

− ε)

gorsza precyzja mantysy

nieużyteczne

(−R

min

− ε, 0, +R

min

− ε) zero

część. użyteczne

(+R

min

− ε, +R

max

)

gorsza precyzja mantysy

użyteczne

(+R

min

, +R

max

)

pełna precyzja mantysy

nieużyteczne

(+R

max

, +∞)

przepełnienie

Macierze i operacje na nich, techniki
przetwarzania macierzy.

Macierz A to tablica elementów a

na której można wykonać operacje macierzowe.

Wskaźnik i = 1, . . . m oznacza numer wiersza, a k = 1, . . . n numer kolumny w
której znajduje się element. Liczby m × n nazywane są rozmiarem macierzy.

A =












. . . a

...

. . . a












Poznaliśmy na wcześniejszych zajęciach symbole operacji macierzowych: ^ potę-
gowanie, ’ sprzężenie, .’ transpozycja, * mnożenie, / dzielenie prawostronne, \
dzielenie lewostronne, + odejmowanie, - odejmowanie. Teraz zajmiemy się nimi
bliżej.

3.1.

Mnożenie macierzy

Mnożenie macierzy A o rozmiarze m × k przez macierz B o rozmiarze r × n jest
możliwe jeśli k = r, tzn kiedy ilość kolumn (k) w pierwszej macierzy jest równa
ilości wierszy (r) w drugiej macierzy. Definicja mnożenia macierzy (symbol * w

Matlab -ie) jest następująca:

C = AB ⇐⇒ c

i,j

k=1

i,k

k,j

Każdy element i-tego wiersza pierwszej macierzy jest mnożony przez j-tą kolumnę
drugiej macierzy. Suma wyników mnożeń zapisywana jest w macierzy wynikowej
na pozycji (i, j).
Własności mnożenia macierzy:

A(BC) = (AB)C,

A(B + C) = AB + AC,

ale na ogół AB 6= BA.

Wydanie polecenia C=A*B w Matlab -ie spowoduje wykonanie następującego
algorytmu:

Algorytm 3.1 Mnożenie macierzy.

Wynik: C[m x n]
Dane:

A[m x r], B[r x n]

C = zeros(m,n)
for j = 1:m

for i = 1:n

for k = 1:r

C(i,j) = A(i,k) * B(k,j) + C(i,j);

end

3.2.

Potęgowanie i dzielenie macierzy

Operacja potęgowania (symbol ^ w Matlab -ie) to krótszy zapis mnożenia macie-
rzy przez samą siebie, czyli A

= AA, a A

= (AA)A, gdzie A musi być macierzą

kwadratową.

Operacja odwrotna do mnożenia to dzielenie macierzy. Jeśli C = AB to podobnie
jak dla wielkości skalarnych z dzielenia prawostronnego (symbol /) otrzymamy
A = C/B oraz przy dzieleniu lewostronnym (symbol \) B = A \ C. Wynik
dzielenia macierzy może być nieokreślony.

Ćwiczenie:

Pomnóż i podziel macierze

Niech

A =











−4

2 −2

0 −1

−1











w =








−2








x =








−3

2
4








Oblicz y = Aw i z = Ax. Następnie utwórz macierz B = [w, x] i wykonaj
mnożenie C = AB. Porównaj kolumny macierzy C z macierzami y i z.

Oblicz D = A \ C i porównaj z macierzą B. Podobnie sprawdź, czy E = C/B
równa się A?

Oblicz F = A \ A, oraz G = A/A. Zauważ, że F i G różnią się rozmiarami. Czy
wiesz jakie wartości będą miały elementy na przekątnej macierzy H = B/B?

3.3.

Transpozycja i sprzężenie macierzy

Dla dowolnej macierzy A o elementach a

istnieje macierz transponowana A

której kolumny są wierszami macierzy A. W Matlab -ie symbolem operatora
transpozycji są dwa znaki występujące jeden za drugim (bez spacji) .’

A =












. . . a

...

. . . a












→












. . . a

...

. . . a












Operator sprzężenia ’ działa tak samo jak operator .’ na macierzach o elementach
rzeczywistych, dlatego używany jest zamiennie. Jeśli elementami macierzy są liczby
zespolone to sprzężenie macierzy powoduje jej transpozycję i zamianę każdego z
elementów na sprzężony.

Przykład:

>> A=[1+2i 2 3 4; 5 6 7 8].’

A =

1.0000 + 2.0000i

5.0000

2.0000

6.0000

3.0000

7.0000

4.0000

8.0000

element a

= 1 + 2i jest liczbą zespoloną

>> A=[1+2i 2 3 4; 5 6 7 8]’

A =

1.0000 - 2.0000i

5.0000

2.0000

6.0000

3.0000

7.0000

4.0000

8.0000

liczba sprzężona do 1 + 2i to 1 − 2i

3.4.

Rząd macierzy

Rząd macierzy (z ang. rank) A o rozmiarze m × n to ilość wierszy lub kolumn li-
niowo niezależnych (wykłady). W Matlab -ie istnieje polecenie rank, które liczy
rząd macierzy (dalej oznaczany przez r). Podobnie jak wszystkie inne oblicze-
nia numeryczne, obliczenie rzędu macierzy są obarczone błędami nieunikniony-
mi. Rząd macierzy o rozmiarze m × n spełnia nierówność r ¬ min{m, n}. Jeśli
r = min{m, n} mówimy, że macierz ma rząd najwyższy.

Przykład:

>> A=[1 0 1; 2 1 0; 3 1 1]

A =

>> rank(A)
ans =

ta macierz ma niepełny rząd ponieważ
2 < min{3, 3}

3.5.

Wyznacznik macierzy

Dla każdej macierzy kwadratowej A o rozmiarze n × n istnieje liczba nazywana
wyznacznikiem, którą oznaczamy symbolem det(A). Jeśli det(A) 6= 0, to macierz
nazywa się nieosobliwą. Jeśli, na przykład, rozmiar macierzy jest 2 × 2 to wy-
znacznik

det(A) =

= a

− a

W ogólności definicję wyznacznika możemy zapisać poprzez minory M

, które są

wyzacznikami macierzy o rozmiarze n − 1 × n − 1 utworzonymi przez zmazanie
i-tego wiersza i j-tej kolumny

det(A) =

j=1

(−1)

i+j

gdzie i oznacza jeden z wierszy.

Wyznacznik w Matlab-ie obliczymy wydając polecenie det(A).

3.6.

Macierze specjalne

Ponieważ macierze mogą być dowolnymi tablicami dwuwymiarowymi o elemen-
tach liczbowych to istnieją macierze specjalne mające określone własności. Znając
własności tych macierzy można uprościć obliczenia. Matlab dostarcza zestaw
komend do ich tworzenia.

3.6.1.

Macierz diagonalna (przekątniowa)

Macierz diagonalna ma elementy jedynie na głównej przekątnej

C = diag(c

, c

, . . . , c

) =












. . .

...

. . . c












Jeżeli pomnożymy dowolną macierz przez macierz diagonalną to otrzymamy efekt
skalowania kolumn, a mnożąc w odwrotnej kolejności skalowania wierszy.

Przykład:

>> x = [1 -5];
>> A = diag(x)
A =

-5

Polecenie Matlab -a diag służy do two-
rzenia macierzy diagonalnej jeśli argu-
mentem wejściowym jest macierz, której
jeden z wymiarów jest równy 1 (wektor).

>> B = [1 -5; 6 12];
>> y = diag(B)
y =

Polecenie diag zwraca elemnty z przekąt-
nej, o ile argumentem jest macierz o wy-
miarach różnych od 1.

>> C = ones(2,2};

>> C*A
ans =

-5

Polecenie ones wypełnia macierz jedyn-
kami. Dalej pokazany jest efekt sklowa-
nia kolumn poprzez mnożenie z macierzą
diagonalna A. Mnożąc AC, otrzymamy
efekt sklowania wierszy.

Ćwiczenie:

Macierze wielo-diagonalne.

Poleceni diag(x,k) ma dwa parametry. Pierwszy z nich x to macierz elemntów
leżących na przekątnej, a drugi to numer przekątnej macierzy. Brak k oznacza do-
myślnie 0. Numer 0 przypisany jest do głównej przekątnej, liczby ujemne wskazuja
na przekątne poniżej, a dodatnie powyżej głównej przekątnej. Wykonaj polecenie:

A = diag(2*ones(4,1)) - diag(ones(3,1),1) - diag(ones(3,1),-1)

B = tridiag(-1,2,-1)

Tę samą macierz trójdiagonalną moż-
na otrzymać poleceniem tridiags, któ-
re znajdziesz w podkatalogu Matlab -a
...\toolbox\mn\macierze\

3.6.2.

Macierz jednostkowa

Macierz jednostkowa to macierz diagonalna o elementach na głownej przekątnej
równych 1

I = diag(1, 1, . . . , 1) =












1 0 . . . 0
0 1 . . . 0

... ... ... ...

0 0 . . . 1












Jeżeli pomnożymy dowolną macierz przez macierz jednostkową to otrzymamy tę
samą macierz.

A I = I A = A

>> I = eye(n)

Polecenie eye(n) tworzy macierz jed-
nostkową o rozmiarze n × n.

>> I = eye(m,n)

Polecenie eye(m,n) tworzy macierz po-
dobną do jednostkej o rozmiarze m × n.

3.6.3.

Macierz permutacji

Jeżeli w macierzy jednostkowej przestawimy wiersze lub kolumny to otrzymamy
macierz permutacji P , na przykład:

P =








0 1 0
0 0 1
1 0 0








Mnożąc dowolną macierz przez macierz pertmutacji otrzymamy taką samą macierz
z tym, że wiersze tej macierzy (lub kolumny) ulegną przestawieniu w taki sam
sposób jak przestawione są wiersze (lub kolumny) w macierzy permutacji.

Przykład:

>> A = [1 2 3 4; 5 6 7 8; 9 10 11 12];
>> P = [1 0 0; 0 0 1; 0 1 0];

Macierz dowolna A i macierz permutacji
P

>> B = P*A

B =

>> C = A’*P
C =

Realizacja efeku przestawienia wierszy i
kolumn przy uzyciu macierzy permutacji.

3.6.4.

Macierz trójkątna

Przez macierz trójkątną rozumiemy macierz postaci:

L =












. . .

...

. . . l












lub

R =












. . . r

...

. . . r












Macierz L nazywamy macierzą trójkątną lewą (lub dolną), a macierz R - prawą
(lub górną). Suma, iloczyn i odwrotność macierzy trójkątnej tego samego rodzaju
są dalej macierzami trójkątnymi. Macierze trójkątne mają wiele innych ciekawych
własności, które są wykorzystywane przy rozwiązywaniu układów równań.

3.6.5.

Macierz odwrotna

Macierz odwrotna do A onaczana jest przez A

−1

i jest to macierz spełniająca

następujący warunek:

−1

A = A A

−1

= I.

Macierz odwrotna jest związana z rozwiązywaniem układów rownań liniowych o
których będzie mowa na następnych zajęciach.

Polecenie Matlab -a inv(A) liczy macierz odwrotną A

−1

3.6.6.

Macierz symetryczna

Macierz symetryczna to taka macierz dla której spełniony jest następujący waru-
nek:

A = A

Jej elemnty są symetryczne względem głównej przekątnej a

= a

Przykład:

>> B=[1 2; 3 4; 5 6];

Macierz B nie jest symetryczna.

>> B*B’

ans =

>> B’*B
ans =

Jednym ze sposobów generacji macierzy
symetrycznej z dowolnej macierzy jest
pomnożenie jej przez macierz transpono-
waną.

3.6.7.

Macierz ortogonalna

Macierz ortogonalna Q to taka macierz dla której spełniony jest następujący wa-
runek:

Q = I.

Układy równań liniowych.

Układ n równań liniowych o n niewiadomych x

, x

, . . . , x

k=1

i,k

= b

(i = 1, 2, . . . , n)

lub

AX = B

ma rozwiązanie wtedy i tylko wtedy, gdy macierz układu A i macierz rozszerzona
R = [A, B] mają ten sam rząd w przeciwnym przypadku są sprzeczne

A =












. . . a

...

. . . a












R =












. . . a

...

. . . a












Jeśli rząd macierzy układu rank(A) = n, to dla dowolnej macierzy B istnieje
rozwiązanie układu X i jest ono jednoznaczne.

Przykład:
Poszukajmy współczyników równania kwadratowego, które aproksymuje krzywą
pompy.

0.5

1.5

100

105

110

115

120

q(m

/s) * 10

−3

h(m)

krzywa pompy

q(m

/s) h(m)

1 × 10

−4

116.2

8 × 10

−4

106.4

14 × 10

−4

82.4

Model (h - wysokość pompowania, q - wydatek pompy):

h = c

+ c

q + c

gdzie c

, c

są szukanymi niewiadomymi.

Podstawmy trzy wybrane punkty (q,h) (patrz rysunek) do modelu

116.2 =

1 × 10

−8

1 × 10

−4

+ c

106.4 =

64 × 10

−8

8 × 10

−4

+ c

82.4 = 196 × 10

−8

+ 14 × 10

−4

+ c

Po przepisaniu do formy macierzowej otrzymamy:

Ax = b

gdzie

A =








1 × 10

−8

1 × 10

−4

64 × 10

−8

8 × 10

−4

196 × 10

−8

14 × 10

−4








x =















b =








116.2
106.4

82.4








Ćwiczenie:
Zdefiniuj macierz A, b i wykonaj A

−1

b poleceniem: >> x = inv(A)*b. Otrzy-

masz c

= −2 × 10

, c

= 4 × 10

= 116. Wykonaj wykres funkcji h(q) w

zakresie q = (0, 1.5 × 10

−3

4.1.

Rozwiązanie formalne dla kwadratowej macierzy ukła-
du.

Formalne rozwiązanie układu równań A X = B to

X = A

−1

gdzie rozmiar A jest n × n.

Nie rozwiązuj A X = B szukając A

−1

i mnożąc przez B

Jeśli widzisz: X = A

−1

B, to szukaj: rozwiązania A X = B przez eliminacje

Gaussa lub inny podobny algorytm.

4.2.

Eliminacja Gaussa.

Celem eliminacji Gaussa jest zamiana dowolnego układu równań (n × n) na rów-
noważny układ, którego macierz jest trójkątna. Przyjrzyjmy się przykładom roz-
wiązań dla układów o macierzy diagonalnej i trójkątnej. Do każdego przypadku
prezentowany jest przykładowy algorytm, który należy szczgółowo przeanalizować.

4.2.1.

Rozwiązanie układu diagonalnego

Przykład:

Układ zdefinowany przez

A =








1 0 0
0 3 0
0 0 5








b =








−1

−15








jest równoważny z

−1

= −15

→

= −1

6
3

= 2

−15

= −3.

Algorytm 4.1 Układ diagonalny.

Wynik: x[n x 1]
Dane:

A[n x n], b[n x 1]

for i = 1:n

x(i) = b(i) / A(i,i)

end

W Matlab -ie:

>> A = ...

>> b = ...
>> x = b./diag(A)

A to macierz diagonalna. To jest jedy-
ny przypadek, gdzie dzielenie elementu
przez element (./) jest wykorzystane do
rozwiązywania układu równań.

4.2.2.

Rozwiązanie układu trójkątnego

Przykład:

Układ zdefinowany przez

A =








−2 1

0 3 −2
0 0








b =








−1








jest równoważny z

−2x

+ x

+ −2x

= −1

→

8
4

= 2

1
3

(−1 + 2x

) =

3
3

= 1

−2

(9 − x

− 2x

) =

−2

= −2

Algorytm 4.2 Układ trójkątny górny. Proces rozwiązywania polega na za-
stosowaniu metody podstawiania w odwrotnej kolejności poczynając od x

Wynik: x[n x 1]
Dane:

U[n x n], b[n x 1]

x(n) = b(n) / U(n,n)
for i = n-1:1

s = b(i)
for j = i+1:n

s = s - U(i,j) * x(j)

end
x(i) = s / U(i,i)

end

Algorytm 4.3 Układ trójkątny dolny. Proces rozwiązywania polega na za-
stosowaniu metody podstawiania w kolejności poczynając od x

do x

Wynik: x[n x 1]
Dane:

L[n x n], b[n x 1]

x(1) = b(1) / L(1,1)
for i = 2:n

s = b(i)
for j = 1:i-1

s = s - L(i,j) * x(j)

end
x(i) = s / L(i,i)

end

4.2.3.

Eliminacja Gaussa

Przykład:

Rozwiążemy układ o macierzy:

A =








−3

2 −1

6 −6

3 −4








b =








−1
−7
−6








W tym celu utworzymy macierz rozszerzoną

A = [Ab] =








−3

2 −1 −1

6 −6

7 −7

3 −4

4 −6








a następnie będziemy przekształcać wiersze macierzy ˜

A tak by otrzymać macierz

trójkątną.

Dodajmy 2 razy wiersz 1 do wierza 2, a następnie dodajmy 1 raz wiersz 1 do 3,
by poniżej głownej przekątnej w kolumnie 1 otrzymac zera.

(1)








−3

2 −1 −1

0 −2

5 −9

0 −2

3 −7








Teraz zajmiemy się kolumną drugą. Należy odjąć od wiersza 3 wiersz 2.

(2)








−3

2 −1 −1

0 −2

5 −9

0 −2








→

−2

= −1

−2

(−9 − 5x

) = 2

−3

(−1 − 2x

+ x

) = 2

Algorytm 4.4 Eliminacja Gaussa.

Wynik: U[n x n]
Dane:

A[n x n]

U=[A b]
for i = 1:n-1

for k = i+1:n

for j = i:n+1

U(k,j) = U(k,j) - (U(k,i) / U(i,i)) * U(i,j)

end

4.2.4.

Eliminacja Gaussa z wyszukiwaniem elementu głównego

Przykład:

Rozwiążemy inny układ o macierzy:

A =











4 −2 −2

4 −3

−3 −3

8 −2

−1

6 −3











b =











−4

5
7
7











Podobnie jak poprzednio utworzymy macierz rozszerzoną











4 −2 −2 −4

4 −3

−3 −3

8 −2

−1

6 −3











i tak jak poprzednio spróbujemy utworzyć macierz trójkątną.

Odejmiemy 1/2 raza wiersz 1 od wiersza 2, dodamy 3/2 raza 1 do 3 i 1/2 raza
wiersz 1 do 4.











2 4 −2 −2 −4
0 0

5 −2

0 3

5 −5

0 3

5 −4











okazuje się że element (22) na głównej przekatnej (element główny) jest równy
zero co uniemożliwia dalsze operacje. Dlatego musimy wykonać dodatkowe prze-
stawienie wiersza 2 i 4











2 4 −2 −2 −4
0 3

5 −5

0 3

5 −4

0 0

5 −2











Kontynuujemy eliminację: odejmujemy 1 raz wiersz 2 od 3. Ponownie na przekątnej
pojawia się zero zatem przestawiamy wiersz 3 z 4.











2 4 −2 −2 −4
0 3

5 −4

0 0

0 −1 −4

0 0

5 −2











→











2 4 −2 −2 −4
0 3

5 −4

0 0

5 −2

0 0

0 −1 −4











Dalej już postępujemy tak jak poprzednio.

Algorytm 4.5 Eliminacja Gaussa z wyszukiwaniem elementu głównego.

Wynik: U[n x n]
Dane:

A[n x n]

U=[A b]
for i = 1:n-1

for k=i:n

poszukaj numer wiersza j gdzie max(|U(j,i)|)>=max(|U(k,i)|)
przestaw wiersz U(j,:) z wierszem U(i,:)

end
for i = i+1:n

for j = i:n+1

U(k,j) = U(k,j) - (U(k,i) / U(i,i)) * U(i,j)

end

4.2.5.

Rozwiązanie układu z użyciem opratora dzielenia lewo-
stronnego.

Dla równania skalarnego wiemy, że na przykład

5x = 20 → x = (5)

−1

20 = 4

Analogicznie możemy zapisać dla układu równań

Ax = b → x = A

−1

gdzie A

−1

b jest formalnym rozwiązaniem równań.

W notacji Matlab -a polecenie rozwiązania układu ma postać następującą

x = A\b

4.3.

Co robi operator dzielenia lewostronnego?

Dla zadanej macierzy A o rozmiarze n × n i wektora b operator \ wykonuje szereg
testów macierzy A. Matlab wybiera metodę rozwiązania układu tak by błędy
numeryczne i ilość operacji były najmniejsze.

1. Sprawdza, czy macierz można przedstawić w postaci trójkątnej? Jeśli tak, to

przyjmuje metodę trójkątów.

2. Sprawdza, czy A jest symetryczna i dodatnia? Jeśli tak, to przyjmuje metodę

Choleskiego.

3. Jeśli metoda Choleskiego zawiedzie lub jeśli macierz jest niesymetryczna, to

Matlab stosuje metodę rozkładu macierzy A na dolną i górną LU .

Limity maszynowe:

•

zapotrzebowanie na RAM ∼ n

•

czas obliczeń ∼ n

Aproksymacja.

Zajmiemy się szukaniem krzywej y = F (x), dla zadanej liczby punktów (x

, y

)

dla i = 1, . . . , m

F (x) = c

(x) + c

(x) + · · · + c

(x),

gdzie c

to nieznane współczynniki, a funkcje f

(x) nazywane są funkcjami bazo-

wymi, gdzie j = 1, . . . , n.
Aproksymację stosujemy wówczas, gdy ilość zadanych punktów m jest większa
od ilości nieznanych współczynników n krzywej F (x). Zwykle nie można przepro-
wadzić krzywej przez wszystkie zadane punkty. Szuka się wówczas ”najlepszej”
krzywej, w sensie minimum kwadratu błędów.

5.1.

Aproksymacja linią prostą

Rozważmy najprostszy przypadek, kiedy f

(x) = x, f

(x) = 1, a pozostałe funkcje

(x) = 0. Wówczas dla kolejnych punktów (x

, y

) mamy

= y

...

= y

lub












... ...































...












Zatem w postaci macierzowej zapiszemy Ac = y. Ponieważ równanie to dla m > n
nie ma dokładnego rozwiązania możemy zapisać

r = y − Ac,

gdzie r jest wektorem pionowych odległości pomiędzy szukaną krzywą i zadanymi
punktami. Będziemy szukali takiego rowiązania dla którego

i=1

lub macierzowo

ma wartość minimalną.

Ponieważ r = y − Ac, to iloczyn r

r możemy zapisać w następującej formie:

r = (y − Ac)

(y − Ac)

= y

y − y

Ac − c

y + c

= y

y − 2y

Ac + c

Ac.

Iloczyn ten osiąga minimum jeśli

r) = 0

→

−A

y + A

Ac = 0,

co ostatecznie prowadzi do równaia

A)c = A

które nazywane jest równaniem aproksymacji.

Przykład:

Aproksymacja liniowa.

Szukamy linii prostej aproksymującej punkty (1 1), (2 2), (4 2), (5 3).

>> x=[1 2 4 5]; y=[1 2 2 3];
>> x=x’; y=y’;
>> A=[x ones(size(x))];
>> c=(A’*A)\(A’*y)
c =

0.4000
0.8000

Definiujemy wektory x i y ze współrzęd-
nymi punktów. Ponieważ dalej potrzeb-
ne są wektory kolumnowe, wykonujemy
transpozycję. Następnie wypełniamy ma-
cierz A, która w pierwszej kolumnie ma
wektor x, a w drugiej wektor jedynkowy
o rozmiarze macierzy x. Macierz współ-
czynników c liczymy wykorzystując rów-
nanie aproksymacji i operator dzielenia
lewostronnego \.

>> ya=c(1)*x+c(2)
>> plot(x,y,’o’,x,ya,’-’)
>> grid on
>> xlabel(’x’)
>> ylabel(’y=F(x)’)

Znając współczynniki c

, c

możemy ob-

liczyc nowe aproksymowane wartości ya
i przedstawić wyniki graficznie.

5.2.

Aproksymacja

Można wykazać, że równanie aproksymacji

A)c = A

jest prawdziwe dla dowolnej funkcji aproksymacji w postaci:

F (x) = c

(x) + c

(x) + · · · + c

(x),

gdzie c

to nieznane współczynniki, f

(x) to funkcje bazowe, a macierze A, c i y

mają następującą postać:

A =












) f

) . . . f

)

) f

) . . . f

)

...

) f

) . . . f

)












c =












...























...












Przykład:

Aproksymacja

Szukamy funkcji w postaci

y =

+ c

która aproksymuje następujące punkty:

x 0.955 1.380 1.854 2.093 2.674 3.006 3.255 3.940 4.060

y 5.722 4.812 4.727 4.850 5.011 5.253 5.617 6.282 6.255

Znajdziesz je w pliku ...\mn\dane\aprox.dat.

>> load aprox.dat
>> x=aprox(:,1)
>> y=aprox(:,2)

Ładowanie danych z pliku i kopiowanie
danych do odpowiednich wektorów ko-
lumnowych x, y.

>> A=[1./x x];
>> c=(A’*A)\(A’*y)
c =

4.2596
1.3008

>> xa=linspace(min(x),max(x),100);
>> xa=xa’;
>> Aa=[1./xa xa];
>> ya=Aa*c;
>> plot(x,y,’o’,xa,ya,’-’)
>> xlabel(’x’)
>> ylabel(’y=F(x)’)
>> legend(’dane’,’aproksymacja’)

Przygotowanie macierzy A i rozwiąza-
nia równania aproksymacji. Polecenie
linspace generyje wektor wierszowy o
elementach równoodleglych w zadanym
zakresie. Zmienne z literą ”a” w nazwie
odnoszą się do wartości aproksymowa-
nych.

5.3.

Co jeszcze robi operator dzielenia lewostronnego?

Rozwiązanie równania macierzowego Ax = y, polega na wyborze właściwej meto-
dy i jej zostosowaniu. Wiemy już, że dla zadanej macierzy A o rozmiarze n × n
i wektora y, operator \ wykonuje szereg testów macierzy A. Matlab wybiera
metodę rozwiązania układu tak by błędy numeryczne i ilość operacji były naj-
mniejsze.

Matlab dzięki wbudowanej wiedzy z algebry potrafi znacznie więcej.

x = A\y = (A’*A)\(A’*y)

W przypadku kiedy macierz układu A ma rozmiar m × n, gdzie m > n au-
tomatycznie zakłada, że użytkownik chce otrzymać aproksymację rozwiązania z
najmniejszym błędem w sensie minimum kwadratu błędów. Zatem w Matlab -
ie równanie aproksymacji jest wywoływane bez udziału użytkownika, wystaryczy
użyć magicznego dzielenia lewostronnego.

Interpolacja.

dane
aproksymacja
interpolacja

Interpolacja polega (podobnie jak aproksyma-
cja) na poszukiwaniu funkcji pomiędzy znany-
mi punktami. W przeciwieństwie do aproksy-
macji funkcja ta

•

przechodzi przez te punkty.

Jeżeli poszukujemy funkcji poza zakresem za-
danych punktów to mówimy o ekstrapolacji.

Wybrane metody interpolacji.

0.5

1.5

−20

100

120

Najblizszy sasiad

0.5

1.5

−20

100

120

Kawalkami liniowo

0.5

1.5

−20

100

120

Kawalkami kwadratowo

0.5

1.5

−20

100

120

Splajny kwadratowo

Ćwiczenie:

Interpolacja wielomianem.

n−1

(x) = c

n−1

+ c

n−2

+ · · · + c

n−1

x + c

>> rok=[1986 1988 1990 1992 1994 1996];
>> cena=[3.13 3.32 3.38 3.41 3.37 3.44];
>> plot(rok,cena,’o’)
>> A=vander(rok);% funkcja MATLABA [x^n x^n-1 .. 1], n wym. x
>> cena=cena’;

% macierz prawych stron układu

>> c=A\cena

% otrzymasz ostrzeżenie - złe skalowanie

>> xnowe=linspace(min(rok),max(rok),100);
>> ynowe=polyval(c,xnowe);
>> plot(rok,cena,’o’,xnowe,ynowe,’-’)

>> rokumowny=rok-mean(rok) % poprawienie skalowania macierzy A

Wykonaj ponownie wszystkie obliczenia dla zmiennej niezależnej (przesuniętej)
rokumowny. Na kartce ”Zadania do zaliczenia (5)” w polu ”Zadanie 2” wpisz
obliczone c

. W katalogu ...\toolbox\mn\interpolacja znajdziesz pomocną

funkcję Idemo_02.

Ćwiczenie:

Interpolacja wielomianem Lagrange-a.

n−1

(x) = y

(x) + y

(x) + · · · + y

(x),

gdzie

(x) =

k=1,k6=j

x − x

− x

W katalogu ...\toolbox\mn\interpolacja znajdziesz funkcję ILagrange.
Korzystając z tej funkcji, oblicz cenę benzyny w 1991 roku dla danych z poprzed-
niego przykładu.

Ćwiczenie:

Interpolacja wielomianem Newton-a.

n−1

(x) = c

(x−x

)+c

(x−x

)(c−x

)+· · ·+c

(x−x

)(x−x

) · · · (x−x

)

W katalogu ...\toolbox\mn\interpolacja znajdziesz funkcję INewton. Ko-
rzystając z tej funkcji, oblicz cenę benzyny w 1993 roku.

Funkcje Idemo_03 i Idemo_04 demonstrują wykorzystanie funkcji: ILagrange,
INewton.
Uwaga! Nie trzeba rozwiązywać układu równań. Jest to zaleta tych metod.

Porównanie interpolacji:

liczba interpolowanych punktów

Flops

wielomian
Lagrange
Newton

Ćwiczenie:

Zbadaj wpływ stopnia wielomianu użytego do interpolacji na jej

jakość. Skorzystaj z funkcji Idemo_05. Wniosk wpisz na kartce ”Zadania do zali-
czenia (5)” w polu ”Zadanie 3”

−5

6.1.

Interpolacja wielomianami sklejanymi.

Interpolacja wielomianami sklejanymi to praktyczne rozwiązanie konieczności
zwiększania stopnia wielomianu. Zamiast stosować jeden wielomian dla wszystkich
danych punktów, metoda sklejania polega na użyciu wielu wielomianów niskiego
stopnia P

(x) dla przedziału danych x

¬ x ¬ x

i+1

. Punkty łączenia wielomia-

nów będziemy nazywali węzłami. W węzłach możemy żądać spełnienia określonych
warunków np. na ciągłość pochodnych.

Interpolacja wielomianami sklejanymi.

i+1

funkcja interpolacyjna
wezly
dane punkty

i+1

i+2

6.2.

Interpolacja kawałkami liniowa.

−1

−0.8

−0.6

−0.4

−0.2

0.2

0.4

0.6

0.8

plot(x,sin(x))

x=linespace(0,2*pi,100)
x=linespace(0,2*pi,6)

Ćwiczenie:

Znajdź dla danych VT.dat temperaturę złącza dla napięcia 0.8mV,

oraz napięcie na złączu dla temperatury 37 stopni. Wykorzystaj funkcję ILiniowa.

>> load(’VT.dat’)
>> V=VT(:,1);
>> T=VT(:,2);

Ładowanie danych z pliku ”VT.dat” do
zmiennej VT, a następnie przepisanie da-
nych z odpowiedniej kolumny macierzy
VT do V (napięcie) i T (temperatura).
Wynik wpisz w polu ”Zadania 4”.

>> Tszukane=ILiniowa(V,T,0.8)

Wywołanie funkcji ILiniowa dla przy-
padku kiedy T = f (V ).

6.3.

Interpolacja kawałkami sześcienna - Hermite.

(x) = a

+ b

(x − x

) + c

(x − x

)

+ d

(x − x

)

gdzie a

, b

, c

, d

są szukanymi współczynnikami przy następujących warunkach w

węzłach:

•

ciągłość P

i−1

) = P

•

znane wartości pierwszej pochodnej i jej ciągłość P

i−1

) = P

Ćwiczenie:

Skorzystaj z funkcji IHermite i Idemo_06 by określić jakość inter-

polacji od ilości węzłów dla funkcji y = xe

−x

. W polu ”Zadania 5” wpisz wielkość

błędu dla 6 i 10 węzłów.

6.4.

Interpolacja kawałkami sześcienna - splajny.

Metod splajnów jest zbliżona do interpolacji Hermite’a. Nie wymaga ona znajo-
mości pochodnych we wszystkich punktach węzłowych. Natomiast wymaga dodat-
kowo:

•

ciągłości drugiej pochodnej P

i−1

) = P

•

znanej wartości pierwszej pochodnej (nachylenia) na końcach P

) i P

– Ustalone nachylenie: P

) = stała1 i P

) = stała2,

– Nachylenie naturalne: P

) = P