Zadania i scenariusze zaj¦¢ z laboratorium

komputerowego do wykªadu z Matematyki

Obliczeniowej

Leszek Marcinkowski

12 grudnia 2011

Streszczenie

W skrypcie przedstawimy zestawy zada« do odbywaj¡cego si¦ co dwa tygo-

dnie laboratorium komputerowego do semestralnego wykªadu z Matematyki

Obliczeniowej. Zakªadamy, »e zadania zostan¡ rozwi¡zane przy wykorzysta-

niu pakietu oblicze« numerycznych octave albo pakietu matlab.

Spis tre±ci

1 Wst¦p

2 Wst¦pne zapoznanie si¦ z octavem

3 Obliczenia w arytmetyce zmiennopozycyjnej

4 Interpolacja Lagrange'a, bazy wielomianów

5 Splajny kubiczne i liniowe. Interpolacja splajnowa

6 Wielomiany Czebyszewa

7 Kwadratury

8 Rozwi¡zywanie równa« nieliniowych

9 Ukªady równa« liniowych - rozkªady typu LU i LL'

10 LZNK. Rozkªad QR. Metoda Householdera

11 Numeryczne zadanie wªasne

12 Algorytm FFT

13 Przykªadowe projekty zaliczeniowe

Rozdziaª 1

Wst¦p

Celem zada« komputerowych przeprowadzanych w laboratorium komputero-

wym jest przetestowanie numerycznych metod omawianych w czasie wykªadu

i ¢wicze« z Matematyki Obliczeniowej.

Zadania komputerowe przedstawione w tym zbiorze polegaj¡ na imple-

mentacji metod numerycznego rozwi¡zywania zada« z wykorzystaniem pa-

kietu octave, czyli ±rodowiska oblicze« numerycznych-naukowych i przete-

stowaniu funkcji octave'a, które s¡ w stanie rozwi¡za¢ zadania omawiane w

trakcie kursu Matematyki Obliczeniowej.

Cz¦±¢ funkcji octave'a wywoªuje odpowiedni¡ bibliotek¦ numeryczn¡, w

której jest zaimplementowana odno±na - zazwyczaj zaawansowana - metoda.

Ale cz¦±¢ funkcji octave'a jest zaimplementowana wprost w samym octave'ie.

Cz¦sto nie wiemy jakiej metody numerycznej u»ywa octave. Szczególnie, »e

kolejne wersje tego pakietu u»ywaj¡ innych bibliotek, mimo »e nazwa funkcji

- np. rozwi¡zuj¡cej równanie nieliniowe - jest wci¡» taka sama.

Pakiet octave jest zarówno ±rodowiskiem oblicze« numerycznych, jak i

j¦zykiem programowania. Umo»liwia on proste rozwi¡zywanie podstawowych

zada« numerycznych jak: numeryczne obliczenie caªki, rozwi¡zanie zada«

liniowych lub nieliniowych, równa« ró»niczkowych zwyczajnych itp. Mo»na

go u»ywa¢ z linii komend, pisa¢ wªasne skrypty, czy m-pliki (czyli pliki z

implementacjami wªasnych funkcji octave'a).

Pakiet octave jest programem ogólnodost¦pnym jako wolne oprogramo-

wanie. Rozprowadzany jest na zasadach licencji GNU GPL. Warto doda¢, »e

pakiet octave jest odpowiednikiem ±rodowiska matlab, które jest programem

komercyjnym, szeroko stosowanym do oblicze« numerycznych.

Zadania z tego zbioru s¡ sformuªowane tak, »e bez kªopotów mog¡ by¢

wykonane zarówno w octave, jak i w ±rodowisku matlaba.

Zalet¡ octave'a jest to, »e jest on dost¦pny bez opªat. Octave jest dystry-

buowany w wersjach binarnych zarówno pod ró»ne dystrybucje Linuxa, jak i

w wersji binarnej pod system Windows.

Poni»ej zaª¡czamy link do stron octave'a:

1. Gªówna strona octave'a

2. Rozbudowany podr¦cznik do octave'a - w j¦zyku angielskim dost¦pny

on-line.

3. Strona z linkami do plików z octavem

4. Strona octave-forge'a - czyli rozszerze« octave'a

Podstawowym podr¦cznikiem do kursu z matematyki obliczeniowej jest

ksi¡»ka [11]. Innymi wartymi polecenia podr¦cznikami w j¦zyku polskim s¡

np. [8], [3], [17] i [1], [16]. W j¦zyku angielskim opublikowano wiele podr¦cz-

ników do matematyki obliczeniowej, inaczej nazywanej te» metodami nume-

rycznymi, czy analiz¡ numeryczn¡. Niektóre z nich zawieraj¡ zadania kompu-

terowe. Warto poleci¢ ksi¡»ki: [15], [14]. Obie zawieraj¡ wi¦cej materiaªu ni»

standardowy kurs matematyki obliczeniowej na wydziale Matematyki, Infor-

matyki i Mechaniki Uniwersytetu Warszawskiego. Innymi wartymi polecenia

podr¦cznikami i monograami dotycz¡cymi metod numerycznych s¡: [5], [6],

[2], [12], [13], [9], [18], [10], [7] i inne. Niektóre z tych podr¦czników obejmuj¡

tylko cz¦±¢ materiaªu z wykªadu matematyki obliczeniowej lub wykraczaj¡

poza ten materiaª.

Je±li chodzi o pakiet octave - polecamy przeczyta¢ manual, tzn. [4] do-

st¦pny w dziale dokumentacji na stronach octave'a on-line pod adresem:

http://www.gnu.org/software/octave/.

Zadania w tym zbiorze s¡ podzielone na rozdziaªy. Standardowe zaj¦cia

z laboratorium komputerowego do Matematyki Obliczeniowej polegaj¡ za-

zwyczaj na rozwi¡zaniu kilku zada« z zadanego tematu. Pozostaªe zadania

mog¡ zosta¢ ewentualnie zadane do rozwi¡zania samodzielnego.

Zaj¦cia z laboratorium do Matematyki Obliczeniowej na wydziale Mate-

matyki, Informatyki i Mechaniki Uniwersytetu Warszawskiego odbywaj¡ si¦

co dwa tygodnie, wi¦c ilo±¢ rozdziaªów w tym skrypcie przewy»sza standar-

dow¡ liczb¦ zaj¦¢. Prowadz¡cy zaj¦cia mo»e wybra¢, które rozdziaªy omówi¢

na zaj¦ciach. Zadania w poszczególnych rozdziaªach s¡ tak skomponowane,

»e standardowy scenariusz zaj¦¢ dotycz¡cych tematu danego rozdziaªu powi-

nien obejmowa¢ kilka pierwszych zada«. Prowadz¡cy zaj¦cia mo»e dokona¢

samodzielnego wyboru zada«, pomijaj¡c niektóre z nich.

Scenariusze i zadania

komputerowe

Rozdziaª 2

Wst¦pne zapoznanie si¦ z octavem

W przypadku przeznaczenia dwóch laboratoriów scenariusz ka»dych zaj¦¢ po-

lega na rozwi¡zaniu mo»liwie du»ej ilo±ci kolejnych zada« z tego rozdziaªu.

W przypadku jednych zaj¦¢ nale»y dokona¢ wyboru rozwi¡zuj¡c tylko kilka-

najwa»niejszych pozycji z listy zada«.

Zadania w tym rozdziale ilustruj¡ podstawowe operacje, struktury i wªa-

sno±ci octave'a. Przedstawimy w zadaniach octave jako kalkulator naukowy.

Omówimy te» operator dwukropek sªu»¡cy np. tworzeniu indeksów. Prze-

testujemy, jak tworzy¢ macierze, wektory; jak zapisywa¢ zmienne do plików

(czyta¢ z plików) w formatach: tekstowym i binarnym. Sprawdzimy two-

rzenie macierzy z podmacierzy, wycinanie podmacierzy i inne podstawowe

operacje na macierzach - mno»enie, dodawanie, transponowanie, funkcje ma-

tematyczne od macierzy, normy wektorów/macierzy.

Zadania obejmuj¡ równie» tworzenie wykresów funkcji matematycznych

przy pomocy funkcji octave'a plot(). Sprawdzimy te» tworzenie i u»ywanie

skryptów i funkcji (m-pliki) w octavie, oraz podstawowe instrukcje warun-

kowe i p¦tle:

•

if else endif;

•

switch case endswitch

•

while( ) do endwhile;

•

do .. until( );

•

for .. endfor.

Zadania obejm¡ te» wska¹niki do funkcji (function handle) i operator @ -

zwracaj¡cy wska¹nik do funkcji.

Zadanie 1 Octave jako kalkulator. Otwórz sesj¦ octave'a. Zapoznaj si¦

z pomoc¡ do funkcji sqrt() oraz sin(). Policz w octave'ie, ile

wynosi pierwiastek z 5 oraz policz warto±¢ funkcji sin na tym

pierwiastku.

Zadanie 2 Operacje macierzowe. Operator dwukropek.

•

Utwórz wektor z liczbami od 1 do 20 oraz wektor ze wszyst-

kimi liczbami parzystymi od −6 do 4.

•

Utwórz dowolne macierze 3x4 A i 3x5 B, a nast¦pnie macierz

3x8 C, której pierwsze 3 kolumny to A, a kolejne to B.

•

Z macierzy C 'wytnij' podmacierz D skªadaj¡c¡ si¦ z 1 gªów-

nego minora tzn. 3x3 od C(1,1) do C(3,3).

•

Zamie« kolejno±¢ kolumn D.

•

Zamie« kolejno±¢ wierszy D.

•

Wytnij dolnotrójk¡tn¡, a potem górnotrójk¡tn¡ cz¦±¢ ma-

cierzy D

•

Wstaw D z powrotem do C jako gªówny minor.

•

Policz sin(D) = (sin(D

)

od D.

•

Zapisz D do pliku (binarnego i ASCII) - zamie« element

D(1,1) na -100 i wczytaj now¡ macierz do octave'a.

Zadanie 3 Funkcje matematyczne od macierzy. Normy macierzy i

wektorów
Policz dyskretn¡ norm¦ maksimum od (sin(x))

na [0, 1] (wektorowo-

czyli bez u»ycia p¦tli).

Zadanie 4 Wykresy funkcji matematycznych.

•

Narysuj wykres funkcji sin(x) na odcinku [0, 4].

•

Narysuj wykres funkcji sin(x) na odcinku [0, 4] - wykres po-

winien by¢ podpisany, narysowany za pomoc¡ gwiazdek w

kolorze czerwonym.

•

Narysuj w jednym oknie podpisane wykresy funkcji sin(x) i
log(x)

na odcinku [1, 6] podpisane odpowiednio.

Zadanie 5 Znajd¹ przybli»one maksimum i minimum funkcji

f (x) = x ∗ (3 + 2 ∗ cos(x))

na odcinku [−1, 5] bez u»ycia p¦tli, oraz przybli»enia punktów

ekstremalnych.
Powtórz to zadanie dla jakiego± wielomianu stopnia dwa i trzy

np. x

+ x

− x − 4

Zadanie 6 Funkcje w octave'ie, m-pliki, czyli tzw. pliki funkcyjne.

Utwórz funkcj¦ w m-pliku obliczaj¡c¡ dla zadanego x warto±¢

funkcji

f (x) = (sin(x))

Zadanie 7 Zmienne globalne

Utwórz m-plik z funkcj¡ octave'a obliczaj¡c¡ warto±¢ funkcji ma-

tematycznej z parametrem: sin(a ∗ x) - parametr przeka» jako

zmienn¡ globaln¡.

Zadanie 8 Implementacja wektorowa funkcji w octave'ie.

Utwórz m-plik z funkcj¡ obliczaj¡cy warto±¢ funkcji f(x) = 1 +
(cos(x))

dla argumentu b¦d¡cego macierz¡, tzn. je±li X = (x

i,j

)

i,j

macierz wymiaru M × N, to funkcja powinna zwróci¢ Y macierz

wymiaru M × N tak¡, »e

Y = (y

i,j

)

i,j

= f (x

i,j

Narysuj wykres f na odcinku [−1, 5] z wykorzystaniem tylko jed-

nego wywoªania tak zaimplementowanej funkcji octave'a.

Zadanie 9 Funkcje anonimowe.

Utwórz funkcj¦ uchwyt do funkcji anonimowej, która dla danego

argumentu x zwraca warto±¢ równ¡ (sin(x))

Zadanie 10 P¦tle, instrukcje warunkowe, instrukcja printf ().

•

Zapoznaj si¦ z pomoc¡ octave'a do p¦tli while, p¦tli for,

instrukcji warunkowej if oraz funkcji drukuj¡cej napisy na

ekranie printf().

•

Przy pomocy p¦tli

for ( k = . . . )

endfor

while ( warunek stopu )

endwhile

oblicz sum¦

= 1 + ... + N

dla N = 100.

•

U»yj instrukcji warunkowej

i f ( )

endif

by sprawdzi¢, czy otrzymane S

zostaªo obliczone popraw-

nie, tzn. czy otrzymali±my 0.5 ∗ N ∗ (N + 1)

•

wyprowad¹ na ekran komunikat u»ywaj¡c

printf ( )

Rozdziaª 3

Obliczenia w arytmetyce

zmiennopozycyjnej

Zadania z tego rozdziaªu powinny wykaza¢ pewne charakterystyczne wªasno-

±ci arytmetyki zmiennopozycyjnej.

W octavie wykorzystywane s¡ domy±lnie liczby zmiennopozycyjne o po-

dwójnej precyzji, jakkolwiek w najnowszych wersjach octave'a mo»na równie»

sztucznie wymusi¢ u»ywanie zmiennych w precyzji pojedynczej przy pomocy

funkcji single (a) zwracaj¡cej zmienn¡ pojedynczej precyzji z t¡ sam¡ war-

to±ci¡.

Jedn¡ z wa»nych wªasno±ci arytmetyki zmiennopozycyjnej, która wynika

z jej konstrukcji, jest to, »e odejmowanie dwóch warto±ci o tym samym znaku

o maªej ró»nicy mo»e skutkowa¢ du»¡ utrat¡ dokªadno±ci wzgl¦dnej.

Zadanie 1 Funkcje single (a) i eps.

Wywoªaj pomoc do tych funkcji

w octave'ie. Sprawd¹, czy prawd¡ jest, »e 1 + eps obliczone w

arytmetyce podwójnej precyzji w octave'ie jest wi¦ksze od jeden.

Przyjmij, »e a = single(eps) i sprawd¹, czy ponownie 1 + a jest

wi¦ksze od jeden.

Zadanie 2 Epsilon maszynowy w arytmetyce podwójnej precyzji

Wyznacz samodzielnie epsilon maszynowy - czyli najmniejsz¡

liczb¦ w arytmetyce zmiennopozycyjnej tak¡, »e po dodaniu jej

do jeden otrzymujemy liczb¦ wi¦ksz¡ od jeden. B¦dziemy szukali

liczby postaci 2

−t

dla t - ilo±ci bitów mantysy. Porównaj z eps

komend¡ octave'a. Zadanie mo»na te» wykona¢ w C/C++. Czy

otrzymane wyniki s¡ takie same jak te otrzymane w octavie (dla

liczb typu double)?

Zadanie 3 Epsilon maszynowy w arytmetyce pojedynczej precyzji

Powtórz poprzednie zadanie dla arytmetyki pojedynczej precyzji.

Funkcja octave'a single (x) tworzy zmienne takiego typu. Wyko-

rzystuj¡c t¦ funkcj¦ ponownie wyznacz epsilon maszynowy jako

liczb¦ postaci 2

−t

, ale dla liczb w pojedynczej precyzji.

Zadanie 4 Narysuj wykres funkcji f(x) = (x + a) − a na [0, 1] dla ró»nych

warto±ci a = 10

dla k = 1, 2, .., 20. Tutaj wa»ne jest aby oblicza¢

warto±¢ f(x) dokªadnie ze wzorów: b = (x − a), f(x) = b − a,

cho¢ matematycznie f(x) = x.

Zadanie 5 Policz

f (x) = x −

√

1 + x ∗ x

algorytmem wprost wynikaj¡cym z tego wzoru, a nast¦pnie z

wykorzystaniem równowa»nego wzoru

f (x) =

−1

x +

√

1 + x ∗ x

tzn. prosz¦ zastosowa¢:
Algorytm 1

a =

√

1 + x

= x − a

oraz
Algorytm 2

a =

√

1 + x

−1

x + a

dla x = 10

i k = 4, . . . , 10. Czy wida¢ ró»nic¦ w wyniku?

Powtórz zadanie w arytmetyce pojedynczej precyzji, tzn. z wy-

korzystaniem funkcji octave'a single (x).

Zadanie 6 Wykres wielomianu na dwa sposoby

Oblicz warto±ci wielomianu

(x − 2)

= x

− ..... + 16

na siatce równomiernej 1000 punktowej na [2 − a, 2 + a] dla a =
10

−3

za pomoc¡ dwóch algorytmów:

Algorytm 1

a = (x − 2),

(x) = a

;

Algorytm 2

(x) = x ∗ x ∗ x ∗ x − . . . + 16.

Matematycznie f

≡ f

, ale wyniki obliczone w arytmetyce zmien-

nopozycyjnej mog¡ si¦ ró»ni¢.
Narysuj wykresy obu funkcji i policz bª¡d kf

(x) − f

(x)k

∞

, czyli

max

(x(k) − f

(x(k))|

. Tu x(k) = 2 − a + k ∗ h dla h =

a/500

. Wektor x mo»na utworzy¢ w octavie przy pomocy funkcji

octave'a linspace().

Zadanie 7 Powtórz poprzednie zadanie dla arytmetyki pojedynczej precyzji,

tzn. powtórz obliczenia wielomianu (x − 2)

= x

− ..... + 16

odcinku [2 − a, 2 + a] dla a = 10

−3

, 10

−2

, 10

−4

dla zmiennych

w pojedynczej precyzji uzyskanych za pomoc¡ funkcji octave'a:

single (x). Narysuj wykresy wielomianu obliczanego obydwoma

algorytmami.

Zadanie 8 Przybli»enie exp(x) z rozwini¦cia w szereg P

∞
k=0

/k!

. Za odpo-

wiedni¡ aproksymacje exp(x) bierzemy najpierw sto pierwszych

elementów szeregu czyli przybli»amy exp(x) przez

(x) =

k=0

dla N = 100, a potem przybli»amy przez tysi¡c elementów sze-

regu, czyli ustalamy N = 1000.
Sprawd¹ bª¡d wzgl¦dny |F

(x) − exp(x)|/| exp(x)|

dla x od −100

do 100 (np. dla liczb ró»ni¡cych si¦ o dziesi¦¢, czyli −100+k ∗10

dla k = 0, . . . , 20) dla obu warto±ci N.
Czy bª¦dy dla liczb ujemnych i dodatnich s¡ tego samego rz¦du?
Jak zmodykowa¢ powy»sz¡ metod¦ przybli»onego obliczania funk-

cji eksponencjalnej exp(x) dla x << 0 tak aby bª¡d wzgl¦dny byª

na tym samym poziomie co dla x > 0?

Zadanie 9 Policz caªki I

/(5+x)dx n = 0, .., 20

dwoma algorytmami

ze wzoru

+ 5 ∗ I

n−1

= 1/n

Pierwszy algorytm przyjmuje I

= log(6/5)

i oblicza z powy»-

szego wzoru kolejne

= 1/n − 5 ∗ I

n−1

dla n = 1, 2, 3, . . ..
Drugi algorytm wykorzystuje fakt, »e

(n + 1)6

≤ I

≤

(n + 1)5

(3.1)

W tym algorytmie przyjmujemy za I

jak¡kolwiek warto±¢ z tego

przedziaªu np. I

= 1/180

i iterujemy w tyª, tzn. dla n =

30, 29, . . . , 20, . . . , 0

obliczamy

n−1

(1/n − I

Porównaj wyniki obu algorytmów dla 0 ≤ n ≤ 20 oraz sprawd¹

czy wyniki otrzymane w octave'ie speªniaj¡ oszacowanie (3.1) dla
n = 1, . . . , 20

Dlaczego jeden z algorytmów dziaªa zdecydowanie lepiej w aryt-

metyce zmiennopozycyjnej?
Jako dodatkowe zadanie teoretyczne pozostawimy uzasadnienie

wzoru rekurencyjnego i oszacowania (3.1) wykorzystywanych w

algorytmach.

Zadanie 10 Zastosuj algorytm bisekcji (algorytm poªowienia odcinka) dla funk-

cji (x − 2)

liczonej z wzoru na rozwini¦cie dwumianu x

− ... − 8

startuj¡c z a = 2 − 10

−3

a b = 2 + 10

−3

. Jako warunek zako«-

czenia dziaªania algorytmu przyjmujemy, »e bª¡d jest mniejszy

od 10

−20

(za przybli»enie rozwi¡zania przyjmujemy ±rodek da-

nego odcinka w metodzie bisekcji, czyli warunkiem zako«czenia

iteracji jest to, »e dªugo±ci odcinka, w którym jest rozwi¡zanie

powinna by¢ mniejsza od 2 ∗ 10

−20

Czy algorytm zwraca przybli»enie liczby dwa jedynego pierwiastka

tego wielomianu?
Narysuj wykresy tej funkcji obliczone z obu wzorów. tzn. (x−2)

i x

− ... − 8

na odcinkach 2 + [−h, 2 ∗ h] dla h = 10

−3

, 10

−4

, 10

−5

Czy z wykresów wynika, »e ten wielomian ma tylko jedno miejsce

zerowe w otoczeniu dwa?

Rozdziaª 4

Interpolacja Lagrange'a, bazy

wielomianów

W tym rozdziale zajmiemy si¦ interpolacj¡ wielomianow¡. Zadanie interpo-

lacji wielomianowej polega na znalezieniu wielomianu stopnia nie wi¦kszego

od n, speªniaj¡cego n + 1 warunków interpolacyjnych.

W octave'ie istnieje funkcja znajduj¡ca wspóªczynniki wielomianu inter-

polacyjnego Lagrange'a. Jest to funkcja polyt(), czyli funkcja obliczaj¡ca

wspóªczynniki wielomianu interpolacyjnego w bazie pot¦gowej dla zadanych

warto±ci i w¦zªów. Z kolei funkcja polyval() jest funkcj¡ obliczaj¡c¡ warto±¢

wielomianu zadanego poprzez wspóªczynniki w bazie pot¦gowej w jednym lub

równocze±nie wielu punktach - czyli wektorowo. Co wa»ne, obie funkcje s¡ ze

sob¡ zgodne. Trzeba uwa»a¢ na kolejno±¢ wspóªczynników; octave numeruje

wspóªczynniki w bazie pot¦gowej

)

n
j=0

przestrzeni wielomianów stopnia nie wi¦kszego od n w odwrotnej kolejno±ci

tzn.

, x

n−1

, . . . , 1)

czyli np. wektor wspóªczynników (3, 2, 1) odpowiada wielomianowi 3x

+2x+

. Oczywi±cie stopnie« wielomianu, a dokªadniej: dla jakiego n rozpatrujemy

baz¦ pot¦gow¡ dla tego wektora, jest dªugo±ci¡ wektora wspóªczynników po-

mniejszon¡ o jeden.

Zadanie 1 Zapoznaj si¦ z pomoc¡ octave'a do funkcji octave'a polyval().

Oblicz warto±¢ wielomianu x

− 1

dla x = −1, 0, 1.

Zadanie 2 Korzystaj¡c z funkcji polyval() narysuj wykres wielomianu x

x − 2

bez wykorzystania p¦tli - czyli wektorowo.

Zadanie 3 Test funkcji polyt().

Zapoznaj si¦ z pomoc¡ octave'a do funkcji octave'a polyt().
Wykorzystuj¡c funkcj¦ polyt() znajd¹ wielomian interpolacyjny
L

dla funkcji F (x) = sin(x) dla w¦zªów −1, 0, 1. Policz warto-

±ci ró»nicy F − L

w w¦zªach oraz korzystaj¡c z funkcji plot()

narysuj wykresy F i L

na jednym rysunku.

Wykonaj powtórnie to zadanie ale dla w¦zªów −1, 0, 1, 10.
Oblicz warto±ci wielomianu bez u»ycia p¦tli z wykorzystaniem

funkcji polyval().

Zadanie 4 Interpolacja Lagrange'a - zbie»no±¢ ci¡gu wielomianów interpo-

lacyjnych dla w¦zªów równoodlegªych i w¦zªów Czebyszewa:

•

Wykorzystuj¡c funkcj¦ polyt() znajd¹ wielomiany inter-

polacyjne L

dla funkcji f = sin() dla N w¦zªów równo-

odlegªych na [0, 2 ∗ π] dla N = 4, 8, 16, 32, 64.

•

Oblicz dyskretn¡ norm¦ maksimum ró»nicy f − L

siatce tysi¡ca równoodlegªych punktów na tym odcinku, tzn.
e

= max |f (x

) − L

f (x

, gdzie x

to punkty siatki. Po-

licz stosunek e

dla N < 64.

Czy bª¦dy malej¡ do zera? Jak zachowuje si¦ e

Siatk¦ tysi¡ca równoodlegªych punktów na odcinku [a, b]

najpro±ciej utworzy¢ z wykorzystaniem funkcji octave'a:

linspace(a,b,1000).

•

Narysuj na ekranie wykresy sin(x) i tych wielomianów dla

ró»nych N - u»ywaj¡c funkcji polyval() i plot().

Zadanie 5 Powtórz zadanie 4 dla tej samej funkcji i tego samego odcinka

dla w¦zªów Czebyszewa. W¦zªy Czebyszewa to pierwiastki wie-

lomianu:

n+1

(t) = cos((n + 1)arccos(t))

na [−1, 1] odpowiednio przesuni¦te i przeskalowane.
Przetestuj, czy bª¦dy e

dla w¦zªów Czebyszewa s¡ mniejsze ni»

dla w¦zªów równoodlegªych.

Zadanie 6 Napisz funkcj¦ znajduj¡c¡ wspóªczynniki w bazie pot¦gowej wie-

lomianu interpolacyjnego zadanego stopnia dla w¦zªów równo-

odlegªych oraz w¦zªów Czebyszewa dla danej funkcji, odcinka
[a, b]

, tzn. napisz funkcj¦ octave'a (w m-pliku):

function [LN, eN]= LagrInterp (FCN, a , b ,N, type=0)

która dla parametrów:

•

zadanego wska¹nika funkcyjnego FCN (ang. function han-

dle) do funkcji jednego argumentu: function y=f(x),

• a, b

- ko«ców odcinka [a, b],

• N

- stopnia wielomianu interpolacyjnego

•

type - typu w¦zªów 0 - równoodlegªych, 1 - Czebyszewa

zwróci wektor LN - wspóªczynniki L

wielomianu interpolu-

j¡cego funkcj¦ w tych w¦zªach oraz eN przybli»enie dyskretnej

normy maksimum ró»nicy L

f − f

na tym odcinku.

Przybli»enie normy maksimum liczymy na dyskretnej siatce za-

wieraj¡cej tysi¡c punktów.

Zadanie 7 Powtórz zadanie 4 dla obu typów w¦zªów, tzn. powtórz znajdo-

wanie wielomianów interpolacyjnych na w¦zªach równoodlegªych

i w¦zªach Czebyszewa dla funkcji

f (x) = log(1 + x)

na odcinkach

• [0, 1]

• [0, 10]

Czy dla tej funkcji i obu odcinków bª¦dy w normach maksimum

malej¡ wraz ze wzrostem N? Porównaj wyniki otrzymane w

tym zadaniu w obliczeniach z oszacowaniami teoretycznymi bª¦du

interpolacji Lagrange'a.

Zadanie 8 Interpolacja Lagrange'a - przykªad Rungego. Powtórz zadanie 4

dla obu typów w¦zªów, tzn. znajdowanie wielomianów interpo-

lacyjnych na w¦zªach równoodlegªych i w¦zªach Czebyszewa, ale

dla funkcji:

f (x) = 1/(1 + x ∗ x)

na [−5, 5].
Czy obliczone wyniki wskazuj¡ na to, »e obliczone ci¡gi wielomia-

nów interpolacyjnych zbiegaj¡ do f jednostajnie dla obu typów

w¦zªów?

Zadanie 9 Napisz funkcj¦ octave'a obliczaj¡c¡ warto±¢ wielomianu zadanego

w bazie pot¦gowej tzn.

w(x) =

k=0

algorytmem Hornera. Parametrami funkcji b¦d¡

•

wektor wspóªczynników a

•

macierz warto±ci x.

• N

- stopnie« wielomianu (to mo»e by¢ parametr opcjonalny,

domy±lnie przyjmuj¡cy warto±¢ równ¡ dªugo±ci wektora a

minus jeden)

Funkcja ma zwróci¢ warto±ci wielomianu dla warto±ci w x.
Przetestuj funkcj¦ dla wielomianów 1 + x

oraz 1 − 2x + x

dla

w¦zªów równoodlegªych na [−1, 2], tzn. policz warto±ci wielomia-

nów dla kilku warto±ci oraz narysuj wykresy tych wielomianów z

wykorzystaniem tej funkcji.

Zadanie 10 Napisz funkcj¦ octave'a znajduj¡c¡ dla danego wielomianu stop-

nia n:

w(x) =

k=0

oraz danej liczby q wspóªczynniki (b

)

n
k=0

wielomianu

p(x) =

n−1

k=0

oraz warto±¢ r takie, »e

w(x) = (x − q) ∗ p(x) + r,

obliczone z wykorzystaniem algorytmu Hornera.

Zadanie 11 Napisz funkcj¦ octave'a znajduj¡c¡ dla danego wielomianu stop-

nia n:

w(x) =

k=0

oraz liczby q warto±¢ w(q) i pochodnej w

(q)

, obliczone z wyko-

rzystaniem algorytmu Hornera.

Zadanie 12 Algorytm Hornera w bazie Newtona. Ró»nice dzielone.

Zaprogramuj w octavie funkcj¦ ze zmodykowanym algorytmem

Hornera zwracaj¡c¡ warto±¢ wielomianu zadanego w bazie New-

tona dla danych w¦zªów. Parametrami b¦d¡

• x

punkt, w którym obliczamy wielomian (ewentualnie ta-

blica punktów, ale wtedy funkcja te» musi zwróci¢ wektor z

warto±ciami wielomianu w tych punktach),

• N

- stopnie« wielomianu,

•

wektor dªugo±ci N + 1 ze wspóªczynnikami wielomianu w

bazie Newtona.

Przetestuj na kilku prostych przykªadach: dla w¦zªów −1, 0, 1 i

wielomian w(x) = x

, który w bazie Newtona zwi¡zanej z tymi

w¦zªami ma nast¦puj¡c¡ posta¢: x

= (x + 1)x − (x + 1) + 1

Zadanie 13 Napisz funkcj¦ octave'a, która dla danego wielomianu w(x), któ-

rego wspóªczynniki w bazie pot¦gowej znamy, oblicza wspóªczyn-

niki tego wielomianu w bazie Newtona dla zadanych w¦zªów po-

danych w wektorze y. Tzn. parametrami funkcji b¦d¡:

•

wektor a = (a

)

taki, »e

w(x) =

k=0

• y = (y

)

wektor wspóªczynników bazy Newtona.

Funkcja powinna zwróci¢ wektor wspóªczynników b

takich, »e

w(x) =

k=0

gdzie

(x) = Π

k−1
j=0

(x − y

Zadanie 14 Napisz funkcj¦, która dla danego wielomianu w(x), którego wspóª-

czynniki w bazie Newtona (w

)

n
k=0

znamy, oblicza wspóªczynniki

tego wielomianu w bazie pot¦gowej (1, x, x

, . . . , x

)

. Tzn. pa-

rametrami funkcji jest wektor wspóªczynników b = (b

)

takich,

»e

w(x) =

k=0

i wektor y = (y

)

w¦zªów bazy Newtona (w

)

n
k=0

dla

(x) = Π

k−1
j=0

(x − y

Funkcja ma zwróci¢ wektor wspóªczynników a = (a

)

takich, »e

w(x) =

k=0

Zadanie 15 Sprawd¹ eksperymentalnie ile wynosi dla ró»nych warto±ci N =

4, 8, 16, 32, 64, . . .

przybli»enie:

r,N

k=0

∞,[−1,1]

dla {l

}

N
k=0

bazy Lagrange'a dla w¦zªów równoodlegªych na [−1, 1],

tzn. dla x

= −1 + k ∗ h

dla h = 2/N.

Norm¦ maksimum funkcji kfk

∞,[a,b]

liczymy w sposób przybli»ony

obliczaj¡c dyskretn¡ norm¦ maksimum na max{1000, 100 ∗ N}

równoodlegªych punktach z odcinka [a, b].

Zadanie 16 Sprawd¹ eksperymentalnie ile wynosi dla ró»nych N np.

N = 4, 8, 16, 32, 64, . . .

przybli»enie:

c,N

k=0

∞,[−1,1]

dla {l

}

N
k=0

bazy Lagrange'a dla w¦zªów Czebyszewa na [−1, 1],

tzn. dla x

zer wielomianu T

N +1

(x) = cos((N + 1)arccos(x))

Norm¦ maksimum funkcji kfk

∞,[a,b]

liczymy w sposób przybli»ony

obliczaj¡c dyskretn¡ norm¦ maksimum na max{1000, 100 ∗ N}

równoodlegªych punktach z odcinka [a, b].

Zadanie 17 Powtórz dwa poprzednie zadania dla w¦zªów równoodlegªych i

w¦zªów Czebyszewa na odcinku [0, 10] i dla odpowiedniej normy

maksimum na tym odcinku.

Zadanie 18 Policz iloraz

f k

∞,[−5,5]

N
k=0

∞,[−5,5]

dla f = 1/(1+x

)

i L

wielomianu interpoluj¡cego f w w¦zªach

równoodlegªych na [−5, 5] dla N = 10, 20, 40, 80. Tutaj {l

}

N
k=0

baza Lagrange'a dla tych w¦zªów.
Norm¦ maksimum funkcji kfk

∞,[a,b]

liczymy w sposób przybli»ony

obliczaj¡c dyskretn¡ norm¦ maksimum na max{1000, 100 ∗ N}

równoodlegªych punktach z odcinka [a, b].

Zadanie 19 Powtórz poprzednie zadanie dla tych samych funkcji i odcinka

dla w¦zªów Czebyszewa zamiast w¦zªów równoodlegªych.

Rozdziaª 5

Splajny kubiczne i liniowe.

Interpolacja splajnowa

W tym rozdziale zajmiemy si¦ interpolacj¡ splajnow¡, czyli interpolowaniem

danej funkcji za pomoc¡ splajnów - inaczej funkcji gi¦tych.

Skupimy si¦ na splajnach kubicznych, czyli funkcjach, które s¡ klasy C

na odcinku [a, b] i dla danego podziaªu tego odcinka:

a = x

< x

. . . < x

= b

na pododcinki. Te funkcje obci¦te do ka»dego pododcinka [x

, x

i+1

]

s¡ wielo-

mianami kubicznymi.

Zadanie interpolacji splajnami kubicznymi polega na znalezieniu spaljnu

kubicznego s speªniaj¡cego:

s(x

) = y

s(x

) = y

...

s(x

) = y

dla zadanych warto±ci y

. Okazuje si¦, »e tak postawione zadanie nie jest

jednoznaczne; trzeba doda¢ dwa dodatkowe warunki na s. Zazwyczaj s¡ to

odpowiednie warunki brzegowe, tzn. zwi¡zane z warto±ciami s, pierwszych

lub drugich pochodnych s w ko«cach odcinka.

Zadanie 1 Funkcje octave'a spline() and ppval.

Zapoznaj si¦ z pomoc¡ do tych funkcji (help spline i help ppval).
Wykorzystuj¡c te funkcje narysuj wykres splajnu kubicznego s

na podziale równomiernym odcinka [−3, 3] z w¦zªami {x

= k}

dla k = −3, −2, . . . , 3 przyjmuj¡cego warto±ci s

) = (−1)

tych w¦zªach.
Nast¦pnie znajd¹ wspóªczynniki splajnu kubicznego s

na tym

samym podziale odcinka i przyjmuj¡cego te same warto±ci w w¦-

zªach co s

, ale który dodatkowo przyjmuje warto±ci pochodnych

w ko«cowych w¦zªach równe zero, tzn. wywoªaj funkcje spline()

podaj¡c dwie warto±ci wi¦cej.
Nast¦pnie narysuj wykresy splajnów s

i s

na tym samym ry-

sunku.
Czy otrzymali±my te same splajny?
Policz przybli»on¡ norm¦ maksimum ró»nicy s

− s

na odcinku

[−3, 3]

Zadanie 2 Splajn kubiczny bazowy.

Dla danych w¦zªów równoodlegªych {k}

k=−5,−4,...,5

na [−5, 5] na-

rysuj wykres splajnu kubicznego typu not-a-knot (czyli splajnu,

którego wspóªczynniki zwróci funkcja spline() przy najprost-

szym wywoªaniu przez podanie wektora w¦zªów i wektora war-

to±ci w tych w¦zªach, por. help spline) takiego, »e s(0) = 1 i
s(k) = 0

dla w¦zªów k 6= 0.

Okre±l na podstawie wykresu no±nik tego splajnu.

Zadanie 3 Splajn kubiczny o minimalnym no±niku.

Dla danych w¦zªów równoodlegªych {k}

k=−5,−4,...,5

na [−5, 5] na-

rysuj wykres splajnu kubicznego takiego, »e s(−1) = s(1) =
1, s(0) = 4

i s(k) = 0 dla w¦zªów k 6∈ {−1, 0, 1} oraz ma po-

chodne równe zero w w¦zªach skrajnych, tzn. : −5 i 5. Czy poza
[−2, 2]

ten splajn jest równy zero?

Policz przybli»one normy maksimum na [−5, −2] i [2, 5] dla tego

splajnu.

Zadanie 4 Testowanie eksperymentalne rz¦du zbie»no±ci splajnu interpola-

cyjnego kubicznego z hermitowskimi warunkami brzegowymi.
Korzystaj¡c z funkcji octave'a spline() znajd¹ wspóªczynniki in-

terpolacyjnego splajnu kubicznego hermitowskiego S

na N w¦-

zªach równoodlegªych dla funkcji f(x) = sin(x) na odcinku [−π, 2∗
π]

dla N = 2

dla N

= 5

i k = 1, 2, 3, 4, 5.

Nast¦pnie

•

narysuj wykresy funkcji f(x) i splajnów S

dla ró»nych N.

•

oblicz dyskretn¡ norm¦ maksimum na siatce równomiernej

zªo»onej z tysi¡ca punktów na tym odcinku, tzn. e

max

| sin(x

) − S

dla x

punktów siatki.

•

policz równocze±nie wspóªczynnik

. Czy prawd¡ jest, »e

≈ 2

dla jakiego± p caªkowitego np. p = 4, 8 lub 16?

Zadanie 5 Testowanie eksperymentalne rz¦du zbie»no±ci splajnu interpola-

cyjnego kubicznego bez warunków brzegowych (splajn typu not

a knot).
Powtórz zadanie 4, ale dla splajnów interpolacyjnych otrzyma-

nych przez spline() bez podawania warto±ci pochodnych w skraj-

nych w¦zªach. Czy wspóªczynniki

s¡ te same? Tzn. czy

szybko±¢ zbie»no±ci k sin(x) − S

∞

jest taka sama?

Zadanie 6 Testowanie eksperymentalne rz¦du zbie»no±ci splajnu interpola-

cyjnego kubicznego naturalnego (warunek brzegowy - zerowanie

si¦ drugich pochodnych w ko«cach odcinka).
Powtórz zadanie 4, ale dla splajnów interpolacyjnych natural-

nych. Tu trzeba wykorzysta¢ funkcj¦ z octave-forge (czyli rozsze-

rzenia pakietu octave)

pp=csape(x,y,'variational')

- ostatni argument okre±la to, »e splajn b¦dzie naturalny.
Podajemy link do strony www z pomoc¡ do funkcji csape():

http://octave.sourceforge.net/splines/function/csape.html

Zadanie 7 Przykªad Rungego, czyli f(x) = 1/(1 + x ∗ x) i odcinek [−5, 5], a

zbie»no±¢ interpolacji splajnami kubicznymi.
Przetestuj jak w poprzednich zadaniach, czy splajny interpola-

cyjne kubiczne z podanymi warunkami na pochodne w ko«cach

odcinka zbiegaj¡ w normie supremum do f, tzn. korzystaj¡c z

funkcji octave'a spline() znajd¹ wspóªczynniki splajnu interpo-

lacyjnego kubicznego S

na N w¦zªach równoodlegªych dla f na

odcinku [−5, 5] dla N = 2

dla N

= 5

i k = 1, 2, 3, 4, 5 oraz

narysuj wykresy f i tych splajnów dla ró»nych N.

Nast¦pnie oblicz dyskretn¡ norm¦ maksimum na siatce zªo»onej

z tysi¡ca punktów na tym odcinku, tzn. e

= max | sin(x

) −

sin(x

dla x

= −5 + k ∗ 0.01

z k = 0, . . . , 1000.

Policz równocze±nie wspóªczynnik

. Czy

≈ 2

dla jakiego±

caªkowitego?

Zadanie 8 Funkcja octave'a mkpp() . Zapoznaj si¦ z t¡ funkcj¡ (help mkpp()).

Utwórz przy pomocy mkpp() splajn kubiczny s na podziale −1 ≤
0 ≤ 1

odcinka [−1, 1] taki, »e s jest wielomianem trzeciego stop-

nia na caªym odcinku [−1, 1] np. s(x) = x

lub (x + 1)

Zadanie 9 Funkcja octave'a umkpp().

Zapoznaj si¦ z t¡ funkcj¡ (help umkpp()).
Utwórz przy pomocy spline() splajn kubiczny s na podziale
−1 ≤ 0 ≤ 1

odcinka [−1, 1] taki, »e s interpoluje wielomian

trzeciego stopnia na caªym odcinku [−1, 1] np. f(x) = (x + 1)

Nast¦pnie sprawd¹ wspóªczynniki s w bazie {(x − x

)

}

j=0,1,2,3

za pomoc¡ umkpp() na obu przedziaªach, tzn. dla x

= −1

przedziale [−1, 0] i x

= 0

na [0, 1].

Zadanie 10 Znajd¹ za pomoc¡ funkcji octave'a umkpp() wspóªczynniki splajnu

z zadania 3 na wszystkich pododcinkach [k, k + 1] w bazie {(x −
k)

}

j=0,...,3

dla k = −2, . . . , 1 czyli tam, gdzie ten B-splajn ma

no±nik.
Porównaj z wynikami otrzymanymi teoretycznie.

Zadanie 11 Interpolacja splajnami liniowymi.

Dla danego równomiernego podziaªu odcinka [−π, 2∗π] na N po-

dodcinków utwórz za pomoc¡ mkpp() struktur¦ splajnu liniowego
s

interpoluj¡cego funkcj¦ sin(x) w w¦zªach dla N = 3, 6, 9, 18.

•

Narysuj wykresy funkcji sin(x) oraz tych splajnów liniowych

na jednym wykresie

•

Policz przybli»one normy maksimum (na 1000 punktach z

tego odcinka) bª¦du e

= k sin −s

∞

•

Policz wspóªczynnik

. Czy wida¢, »e

≈ 2

dla jakiego± p caªkowitego np. p = 2, 4 lub 8?

W tym zadaniu mo»na wykorzysta¢ funkcj¦ z nast¦pnego zada-

nia, tj. zadania 12.

Zadanie 12 Napisz w octavie funkcj¦ function pp=linspline(x,y), która dla

• x

wektora N + 1 ró»nych w¦zªów uszeregowanych (a = x

. . . < x

= b

)

• y

- wektora N + 1 warto±ci funkcji y = f(x)

zwróci w strukturze pp wspóªczynniki splajnu liniowego s dla

podziaªu zadanego w¦zªami x

Struktur¦ nale»y utworzy¢ funkcj¡ octave'a mkpp() w taki spo-

sób, aby mo»na byªo obliczy¢ warto±¢ tego splajnu w punkcie

(tablicy punktów) za pomoc¡ funkcji octave'a ppval().

Zadanie 13 Testowanie rz¦du zbie»no±ci interpolacji splajnami liniowymi w

zale»no±ci od gªadko±ci funkcji.
Dla funkcji

(x) = (x)

j
+

x < 0

x > 0

j = 1, 2, 3

oraz dla podziaªu odcinka [a, b] z a = −π i b = 3 na w¦zªy

równoodlegªe

= −π + k ∗ h}

dla h = (b − a)/N i N = 4, 8, 16, 32, 64.
Przetestuj rz¡d zbie»no±ci splajnu liniowego interpolacyjnego.
Bardziej szczegóªowo:

•

przy pomocy funkcji octave'a z zadania 12 znajd¹ wspóª-

czynniki odpowiedniego splajnu liniowego s

•

policz przybli»on¡ norm¦ maksimum bª¦du (na 1000 równo-

miernych punktach), tzn. przybli»enie

= ks

− f

∞,[a,b]

•

policz wspóªczynniki

. Czy wida¢, »e

≈ 2

dla jakiego± p caªkowitego np. p = 2, 4 lub 8? Czy wida¢

ró»nic¦ dla ró»nych j?

Rozdziaª 6

Wielomiany Czebyszewa

Wielomiany Czebyszewa deniujemy rekurencyjnie: T

−1

≡ 0

, T

≡ 1

oraz

n+1

(x) = 2xT

(x) − T

n−1

(x)

(6.1)

albo ze wzoru:

(x) = cos(n ∗ arccos(x))

x ∈ [−1, 1].

(6.2)

W tym rozdziale przetestujemy podstawowe wªasno±ci tych wielomianów.

Zadanie 1 Napisz rekurencyjn¡ funkcj¦ octave'a obliczaj¡c¡ warto±¢ wie-

lomianu zadanego poprzez wspóªczynniki w bazie Czebyszewa

wprost ze wzoru rekurencyjnego (6.1), tzn. parametrami funkcji

b¦d¡:

• x

argument dla którego chcemy obliczy¢ warto±¢ wielomianu

• a

- wektor zawieraj¡cy wspóªczynniki a

, . . . , a

takie, »e

w(x) =

n
k=0

(x)

• n

- stopie« wielomianu (parametr opcjonalny, mo»na t¦ war-

to±¢ uzyska¢ z wektora a)

Funkcja octave'a ma zwróci¢ warto±¢ w(x).
Przetestuj t¦ funkcj¦ dla w(x) = T

(x)

, tzn. narysuj dwa wy-

kresy T

(x)

raz korzystaj¡c z tej funkcji oraz drugi raz korzysta-

j¡c wprost ze wzoru (6.2).

Zadanie 2 Napisz nierekurencyjn¡ funkcj¦ octave'a

function [Y]=Czebyszew (X, a , n)

obliczaj¡c¡ warto±¢ wielomianu zadanego poprzez wspóªczynniki

w bazie Czebyszewa z wzoru (6.1), tzn. parametry wej±ciowe

funkcji to:

•

macierz X = (x

)

, wymiaru k × l

•

wektor wspóªczynników

a = [a

, . . . , a

]

takich, »e

w(x) =

k=0

(x).

• n

stopie« wielomianu; ten parametr mo»e by¢ opcjonalny,

je±li go nie podamy to funkcja powinna przyj¡¢, »e n to

wymiar wektora wspóªczynników minus jeden.

Funkcja octave'a powinna zwróci¢ macierz Y = (y

)

, wymiaru

k × l

z warto±ciami y

= w(x

)

Funkcja ma korzysta¢ z wzoru rekurencyjnego (6.1), ale jako

funkcja octave'a ma by¢ nierekurencyjna. Koszt algorytmu

powinien wynosi¢ k ∗ l ∗ O(n).
Narysuj dwa wykres T

(x)

- jeden korzystaj¡c z tej funkcji oraz

drugi - ze wzoru (6.2).

Zadanie 3 Sprawd¹ obliczeniowo, czy powy»sze wzory na wielomiany Cze-

byszewa: wzór deniuj¡cy wielomiany Czebyszewa rekurencyjnie

(6.1) i (6.2), s¡ zgodne na odcinku [−1, 1]. Tzn. dla du»ej ilo±ci

np. 10000 losowych punktów na [−1, 1] policz warto±ci T

wzoru

(x) = cos(32 ∗ arccos(x))

oraz ze wzoru rekurencyjnego (6.1).
Policz maksimum warto±ci absolutnych z ró»nicy mi¦dzy warto-

±ciami wielomianu obliczonymi obiema metodami. Porównaj czas

obliczania warto±ci tego wielomianu obydwoma wzorami.
Powtórz testy ale dla wielomianów innych stopni, np. T

, T

Zadanie 4 Narysuj korzystaj¡c z wzoru

(x) = cos(n ∗ arccos(x))

wykres wielomianów T

dla n = 0, 1, 2, 3, 4, 5. Policz dyskretn¡

norm¦ maksimum na siatce. Czy wynosi jeden? Policz na wy-

kresie zera oraz ekstrema ka»dego z wielomianów.

Zadanie 5 Wyznacz zera i ekstrema T

i T

ze wzoru T

(x) = cos(n ∗

arccos(x))

. Sprawd¹ warto±ci tych wielomianów w jego zerach i

ekstremach.

Zadanie 6 Sprawdzenie wªasno±ci optymalno±ci zer wielomianu Czebyszewa

jako w¦zªów interpolacji Lagrange'a.
Chcemy eksperymentalnie sprawdzi¢, czy wielomian 2

−n

N +1

(x)

ma minimaln¡ norm¦ supremum na [−1, 1] w±ród wszystkich wie-

lomianów postaci Π

N
k=0

(x − x

)

Policz dla N + 1 ró»nych losowych w¦zªów x

, k = 0, . . . , N z

odcinka [−1, 1] dla N = 2, 4, 8, 16, 32 i kilkuset ró»nych losowych

zestawów przybli»on¡ norm¦ supremum wielomianu Π

n
k=0

(x−x

)

Przeprowad¹ testy równie» dla w¦zªów Czebyszewa (czyli zer T

N +1

)

oraz w¦zªów równoodlegªych dla tego samego N.

Zadanie 7 Sprawdzenie wªasno±ci ekstremalnej wielomianów Czebyszewa.

Z teorii wiadomo, »e wielomian 2

−n

N +1

(x)

ma minimaln¡ norm¦

supremum na [−1, 1] w±ród wszystkich wielomianów postaci

w(x) = x

n+1

k=0

Chcemy sprawdzi¢ eksperymentalnie, czy ta wªasno±¢ si¦ potwier-

dzi.
Policz dla N + 1 losowych wspóªczynników a

dla k = 0, . . . , N

przybli»on¡ norm¦ supremum wielomianu:

w(x) = x

n+1

k=0

Przetestuj dla N = 2, 4, 8, 16, 32 i co najmniej kilku tysi¦cy ró»-

nych losowych zestawów w¦zªów. W¦zªy losujemy funkcj¡ randn()

zwracaj¡c¡ liczby losowe o rozkªadzie normalnym.

Zadanie 8 Sprawdzenie wªasno±ci ekstremalnej wielomianów Czebyszewa przyj-

muj¡cych ustalon¡ warto±¢ poza [−1, 1].

Chcemy sprawdzi¢, czy dla dowolnego punktu a = 2 i b 6= 1

wielomian

∗

(x) = T

(x)/T

(2)

ma minimaln¡ norm¦ supremum na [−1, 1] w±ród wszystkich wie-

lomianów w stopnia nie wi¦kszego od n+1 przyjmuj¡cych warto±¢
b

w punkcie a.

Policz dla n losowych wspóªczynników a

z k = 1, . . . , n dla

n = 2, 4, 8, 16, 32

i co najmniej kilku tysi¦cy ró»nych losowych

zestawów w¦zªów przybli»on¡ norm¦ supremum wielomianu

w(x) = 1 +

k=1

(x − 2)

i porównaj z norm¡ supremum w

∗

na [−1, 1] równ¡ 1/T

(2)

(dla-

czego tyle ona wynosi?).
W¦zªy losujemy funkcj¡ randn() zwracaj¡c¡ liczby losowe o roz-

kªadzie normalnym.

Zadanie 9 Funkcja octave'a sªu»¡ca przybli»onemu obliczaniu caªek po od-

cinkach (zadanie pomocnicze)
Zapoznaj si¦ z funkcj¡ octave'a quad().
Policz caªk¦ przy pomocy funkcji quad() z sin(x) na [−π, 2 ∗ π].

Zadanie 10 (trudne) Iloczyn skalarny typu L

2
w

(a, b)

Napisz funkcj¦ octave'a:

function [ nl2 ]=IlSkL2w (FCN,GCN, a , b ,FCNW)

komendy octave ' a

nl2 = . . . .

endfunction

która oblicza iloczyn skalarny typu L

2
w

(a, b)

z wag¡ w(x) tzn.:

(f, g)

(a,b)

f (x)g(x)w(x) dx

dla danej funkcji wagowej w(x) i funkcji f, g.
Parametry funkcji to:

•

wska¹niki do funkcji F CN, GCN do dwóch funkcji octave'a

function y=f ( x )

y = . . . . . .

endfunction

function y=g ( x )

y = . . . . . .

endfunction

obliczaj¡cych odpowiednio warto±ci funkcji f(x) i g(x) dla
x ∈ [a, b]

• a

, b - ko«ce odcinka caªkowania [a, b],

• F CN W

to wska¹nik do funkcji wagowej

function y=w( x )

y = . . . . . .

endfunction

Funkcja powinna zwróci¢ przybli»enie iloczynu skalarnego obli-

czone za pomoc¡ funkcji octave'a quad().
W przypadku wywoªania funkcji IlSkL2() tylko z czterema pierw-

szymi argumentami (tzn. bez podania wska¹nika do wagi) funk-

cja powinna zwróci¢ iloczyn skalarny z wag¡ w ≡ 1.

Zadanie 11 Ortogonalno±¢ wielomianów Czebyszewa w L

√

1−x2

(−1, 1)

Sprawd¹ eksperymentalnie w octavie, np. za pomoc¡ funkcji

quad() lub wªasnej funkcji z poprzedniego zadania, czy wielo-

miany Czebyszewa tworz¡ ukªad ortogonalny w L

√

1−x2

(−1, 1)

tzn. czy prawd¡ jest, »e

−1

(x)T

(x)

√

1 − x

dx =







m 6= n,

m = n > 0,

m = n = 0.

Sprawd¹ powy»sze zale»no±ci dla ró»nych m, n ró»nej wielko±ci

np. dla m = 0, 1, 2, 3, 100, 1000 i n = 0, 3, 10, 500, 1004.

Rozdziaª 7

Kwadratury

Zadania z tego rozdziaªu sªu»¡ przetestowaniu najprostszych kwadratur nu-

merycznych, czyli metod przybli»onego obliczania caªek po odcinkach.

W tym rozdziale zapoznamy si¦ tak»e z funkcj¡ octave'a quad() sªu»¡c¡

przybli»onemu obliczaniu caªek jednowymiarowych postaci

f (t) dt

− ∞ ≤ a < b ≤ +∞,

Mo»emy wi¦c znajdowa¢ przybli»enia caªek po caªej prostej rzeczywistej czy

póªprostych.

Zadanie 1 Zapoznaj si¦ z pomoc¡ do funkcji octave'a quad().

Policz przy pomocy quad() caªk¦ R

f (t) dt

dla nast¦puj¡cych

funkcji i odcinków:

• x

na [−1, 1]

• x

na [0, 1]

• sin(x)

na [0, π]

• cos(100 ∗ x)

na [0, π]

• cos(100 ∗ x) ∗ cos(1000 ∗ x)

na [0, π]

Czy wyniki s¡ zgodne z teori¡?

Zadanie 2 Zªo»ona kwadratura trapezów.

(a) Zaprogramuj w octavie funkcj¦

function c=kwadtrapez(FCN,a,b,n)

obliczaj¡c¡ zªo»on¡ kwadratur¦ trapezów:

f = h

0.5[f (a) + f (b)] +

n−1

k=1

f (a + k ∗ h)

dla h = (b − a)/n.
Parametry funkcji:

• F CN

- wska¹nik do funkcji octave'a function y=f(x)

obliczaj¡cej warto±¢ funkcji podcaªkowej

• a

- lewy koniec odcinka

• b

- prawy koniec odcinka

• n

- ilo±¢ oblicze« warto±ci funkcji podcaªkowej w kwa-

draturze trapezów minus jeden (tak, jak we wzorze po-

wy»ej)

Funkcja zwraca przybli»on¡ warto±¢ caªki obliczon¡ za po-

moc¡ powy»szego wzoru, tzn. zªo»onej kwadratury trape-

zów.
Funkcja powinna dziaªa¢ równie» je±li j¡ wywoªamy tylko z

trzema parametrami. Wtedy n powinno domy±lnie przyj¡¢

warto±¢ sto.

(b) Przetestuj funkcj¦ octave'a z poprzedniego podpunktu dla

= 2

k = 0, 1, 2, . . .

z ustalonym N

= 5

licz¡c:

dla bª¦du

= |

f dt − T

f |

dla nast¦puj¡cych funkcji, dla których warto±ci caªek znamy:

• x

na [0, 1],

• sin(x)

na [0, π] i N

= 5

- funkcja analityczna,

• sin(100 ∗ x)

na [0, π] i N

= 5

- funkcja analityczna,

silnie oscyluj¡ca - du»e warto±ci drugiej pochodnej,

• x

j+0.5

na [0, 1] dla j = 0, 1, 2 czyli funkcji w C

([0, 1])

o nieograniczonej w otoczeniu zera j + 1 pochodnej.

Porównaj wyniki obliczane kwadratur¡ trapezów z wyni-

kami funkcji quad(). Mo»na sprawdzi¢ dla jakiego n bª¡d

obliczony zªo»on¡ kwadratur¡ trapezów jest na poziomie

bª¦du funkcji octave quad() - i porówna¢ ilo±¢ wywoªa«

funkcji f przez obie procedury.

Zadanie 3 Czy wielomiany Czebyszewa tworz¡ ukªad funkcji ortogonalnych

w L

na [−1, 1] z wag¡

√

1−x

Policz za pomoc¡ funkcji octave quad():

−1

√

1 − x ∗ x

(x)T

(x) dx

dla m = 2 i n = 3. Czy wynik jest zgodny z teori¡?

Zadanie 4 Kwadratura Gaussa-Czebyszewa dla

I(f ) =

−1

√

1 − x

f (x) dx

czyli

n+1

f =

n + 1

k=0

f (x

) =

n + 1

k=0

cos

0.5 ∗ π + kπ

n + 1

dla x

zer T

n+1

- n + 1 wielomianu Czebyszewa.

Zaimplementuj funkcj¦ function c=gaussczeb(FCN,n) oblicza-

j¡c¡ warto±¢ kwadratury GC

. Parametry - wska¹nik do funkcji

function y=f(x) i ilo±¢ punktów kwadratury.
Przetestuj jej dziaªanie:

•

policz przybli»enia caªek I(T

)

. Czy GC

)

przybli»a π/2

dla k > 0 dla n = 20, 40, 80?

•

policz przybli»enie caªki z wielomianów ró»nych stopni dla
n = 4, 10, 20

- czy kwadratura jest dokªadna dla wielomia-

nów Czebyszewa stopnia 0 < k < 2n (tzn. czy zwraca zero)?

•

dla funkcji matematycznej f(x) = exp(arccos(x)), której

caªk¦ z t¡ wag¡ mo»emy obliczy¢ dokªadnie.
Porównaj warto±¢ kwadratury z wynikiem dokªadnym dla
n = 10, 20, 40, 80

, tzn. policz

= |GC

(f ) − I(f )|

przyjmuj¡c za I(f) warto±¢ dokªadn¡.

•

analogicznie do poprzedniego podpunktu tego zadania, po-

równaj wynik obliczony za pomoc¡ kwadratury Gaussa -

Czebyszewa, tzn. GC

(f )

, dla funkcji f(x) = exp(arccos(x))

z wynikiem obliczonym funkcj¡ quad(). Wypisz na ekran
e

= |GC

(f ) − c|

dla c warto±ci uzyskanej za pomoc¡ funk-

cji octave'a quad().

•

Przetestuj rz¡d zbie»no±ci tej kwadratury dla funkcji f(x) =
exp(arccos(x))

, jak w przypadku kwadratury trapezów, tzn.

policz bª¦dy e

, analogicznie do poprzednich podpunktów i

wypisz e

Zadanie 5 Zªo»ona kwadratura prostok¡tów.

(a) Zaprogramuj w octave funkcj¦

function c=kwadprost(FCN,a,b,n),
która zwraca warto±¢ zªo»onej kwadratury prostok¡tów dla

funkcji f na odcinku [a, b] na n równoodlegªych punktach:

f = h

k=1

f (a + (k − 0.5) ∗ h)),

h = (b − a)/n.

Parametry funkcji:

• F CN

- wska¹nik do funkcji octave'a postaci

function y=f ( x )

y = . . . . . .

endfunction

obliczaj¡cej warto±¢ funkcji, której caªk¦ chcemy obli-

czy¢, tj. f(x), w danym punkcie x.

• a, b

- ko«ce odcinka,

• n

- ilo±¢ w¦zªów wykorzystywanych przez kwadratur¦.

Funkcja zwraca przybli»on¡ warto±¢ caªki obliczon¡ za po-

moc¡ powy»szego wzoru, tzn. zªo»onej kwadratury prosto-

k¡tów.
Funkcja powinna dziaªa¢ równie» je±li j¡ wywoªamy bez ostat-

niego parametru. Wtedy n powinno domy±lnie przyj¡¢ usta-

lon¡ warto±¢, np. dwie±cie.

(b) Testy dla funkcji, których caªki znamy.

Przetestuj t¦ funkcj¦ octave'a kwadprost() dla N

= 2

k = 0, 1, 2, . . .

z ustalonym N

= 5

licz¡c:

= |

f dt − P

f |

dla caªek z nast¦puj¡cych funkcji po odpowiednich odcin-

kach:

• x

na [0, 1]

• sin(x)

na [0, π] i N

= 5

- funkcja analityczna,

• sin(100∗x)

na [0, π] i N

= 5

- funkcja analityczna silnie

oscyluj¡ca (du»e warto±ci drugiej pochodnej),

• x

j+0.5

na [0, 1] dla j = 0, 1, 2 czyli funkcji w C

([0, 1])

o nieograniczonej w otoczeniu zera j + 1 pochodnej.

Porównaj wyniki obliczane kwadratur¡ prostok¡tów z wyni-

kami funkcji quad() .

Zadanie 6 Porównaj wyniki obliczane kwadratur¡ prostok¡tów z wynikami

funkcji obliczaj¡cej caªk¦ za pomoc¡ zªo»onej kwadratury trape-

zów.
Porównaj bª¦dy dla warto±ci N i f(x) z poprzedniego zadania -

czyli caªek, których warto±¢ teoretycznie znamy - z wynikami dla

obu kwadratur:

•

dla tej samej ilo±ci wywoªa« funkcji f, tzn. porównaj P

z T

N −1

•

porównaj obie kwadratury dla tego samego h, tzn. T

, aby oceni¢ bª¡d w terminach parametru h = (b−a)/N.

Zadanie 7 Kwadratura Romberga.

(a) Zaprogramuj w octave funkcj¦

function c=Romberg(FCN,a,b,p,N0)
z kwadratur¡ Romberga R

p,N 0

obliczaj¡c¡ przybli»enie caªki

f (t) dt

Parametry funkcji:

• F CN

wska¹nik do funkcji octave'a function y=f(x)

obliczaj¡cej warto±¢ funkcji podcaªkowej

• a

lewy koniec odcinka

• b

prawy koniec odcinka

• p

ilo±¢ poziomów w kwadraturze Romberga

• N 0

startowa ilo±¢ punktów w kwadraturze Romberga

Funkcja zwraca przybli»on¡ warto±¢ caªki obliczon¡ za po-

moc¡ kwadratury Romberga.
Funkcja powinna dziaªa¢ równie» je±li j¡ wywoªamy tylko z

trzema albo czterema parametrami. Parametr p powinien

wtedy domy±lnie przyj¡¢ warto±¢ 5, a N0 warto±¢ sto.

(b) Przetestuj kwadratur¦ Romberga dla caªki R

f dt

dla na-

st¦puj¡cych funkcji, odcinków i warto±ci parametrów N0 i
p

• sin(x)

dla [0, π]

• sin(x)

dla [0, 100 ∗ π]

• sin(10 ∗ x)

dla [0, π]

• x

j+0.5

na [0, 1] dla j = 0, 1, 2, . . . , 10,

dla N0 = 100 i dla p = 2, 3, 4, 5, 6.
Czy dla

√

i rosn¡cego p bª¡d maleje?

Czy stosunek E

maleje tak samo dla wszystkich funk-

cji?
Tutaj u»yli±my oznaczenia na bª¡d E

= |

f dt − R

p,N 0

Rozdziaª 8

Rozwi¡zywanie równa«

nieliniowych

W poni»szych zadaniach przetestujemy jak dziaªaj¡ cztery metody rozwi¡-

zywania równa« nieliniowych skalarnych, tzn. metoda bisekcji, metoda New-

tona, metoda siecznych i metoda iteracji prostych oraz dwie metody rozwi¡-

zywania ukªadów równa« nieliniowych, czyli wielowymiarowa metoda New-

tona i metoda iteracji prostych w najprostszej formie.

Zapoznamy si¦ równie» z funkcjami octave'a fzero() oraz fsolve () sªu-

»¡cymi do rozwi¡zywania równa« nieliniowych skalarnych i ukªadów równa«

nieliniowych.

Zadanie 1 Zaimplementuj metod¦ bisekcji w skrypcie - dla rozwi¡zania rów-

nania cos(x) = 0 na odcinku [0, 2]. Przetestuj, czy funkcja znaj-

dzie dobre przybli»enie π/2.

Zadanie 2 Przetestuj metod¦ bisekcji z poprzedniego zadania do rozwi¡zania

innych problemów:

− 2 = 0

exp(x) = 2

cos(10 ∗ x) = 0

startuj¡c z odcinka [0, 110].

Zadanie 3 Napisz funkcj¦ octave'a:

function [ y , i t e r , kod]= b i s e k c j a (FCN, a , b , . . .

tola , maxit ) ,

której parametrami b¦d¡

•

wska¹nik do funkcji FCN (inaczej: function handle),

• a, b

- ko«ce przedziaªu,

•

tola - »¡dana tolerancja bezwzgl¦dna - bª¡d powinien by¢

mniejszy od tola,

•

maxit - maksymalna ilo±¢ iteracji.

Funkcja ma zwróci¢

•

y - przybli»one rozwi¡zanie,

•

iter - ilo±¢ iteracji,

•

kod - kod wyniku:

0 metoda zbiegªa,
1 - metoda zatrzymaªa si¦ z powodu przekroczenia mak-

symalnej ilo±ci iteracji,

2 - warto±ci w ko«cach odcinka funkcji maj¡ ten sam

znak.

Funkcja ma dziaªa¢ równie» je±li podamy tylko trzy lub cztery

pierwsze argumenty - wtedy maksymalna ilo±¢ iteracji domy±lnie

powinna wynosi¢ sto, a tolerancja 10

−5

Zadanie 4 Zapoznaj si¦ z pomoc¡ dla funkcji octave'a fzero() - rozwi¡zu-

j¡c¡ skalarne równania nieliniowe, przetestuj j¡ na kilku prostych

przykªadach skalarnych np.

• cos(x) = 0

• x

− 2 = 0

• x

− 2 = 0

• x − sin(x) = 0

i innych prostych równaniach nieliniowych.
Prosz¦ przetestowa¢ t¦ funkcj¦ dla ró»nych warto±ci pocz¡tko-

wego przedziaªu x

= [a, b]

Zadanie 5 Zapoznaj si¦ z pomoc¡ dla funkcji octave'a fsolve () - rozwi¡zu-

j¡c¡ ukªady równa« nieliniowych. Przetestuj j¡ na kilku prostych

przykªadach skalarnych, np. na tych z poprzedniego zadania.

Zadanie 6 Zaimplementuj metod¦ Newtona w octavie i przetestuj jej zbie»-

no±¢ dla nast¦puj¡cych funkcji:

f1(x)= x ∗ x − 2 z x

= 2

f2(x)= x ∗ x ∗ x − 27 z x

= 27

f3(x)= exp(x) − 2 z x

= 10, −10, 1e3

f4(x)= sin(x) dla x

= 2

f5(x)= cos(x) z x

= 2

f6(x)= (x − 2)

x x

= 3

dla k = 2, 4, 8, 16,

f7(x)= x ∗ x − 2 z x

= 10

, sprawd¹, czy metoda Newtona

zbiega, a je±li tak - to czy zbiega ona kwadratowo,

f8(x)= 1/x − a dla danych a = 0.5, 2, 4, 100 (oczywi±cie imple-

mentuj¡c bez dzielenia) jak dobra¢ x

Dla wszystkich tych funkcji znamy rozwi¡zania, wi¦c mo»na wy-

±wietla¢ równocze±nie na ekranie bª¡d

= x

− r

(tutaj r - znane rozwi¡zanie) i bada¢ eksperymentalnie rz¡d zbie»-

no±ci, tzn. drukowa¢ na ekranie stosunki bª¦du bie»¡cego do po-

przedniego w odpowiedniej pot¦dze:

|/|e

n−1

dla p = 1, 2, 3.
Prosz¦ powtórzy¢ testy z innymi ró»nymi warto±ciami startowymi
x

Zadanie 7 Powtórz zadanie 6 zast¦puj¡c metod¦ Newtona przez metod¦

siecznych.
Przetestuj, czy

/(e

n−1

n−2

)

asymptotycznie zbiega do staªej, i czy

|/|e

n−1

dla p = (1 +

√

5)/2

zbiega do staªej np. dla x ∗ x − 2 z x

= 2

lub innych równa« z poprzedniego zadania.
Za x

mo»emy przyj¡¢ x

+ 10

−7

Zadanie 8 Sprawd¹, czy metoda iteracji prostych

= cos(x

n−1

)

zbiega do x

∗

= cos(x

∗

)

Zbadaj eksperymentalnie, czy zbie»no±¢ jest liniowa, tzn. czy

∗

− x

∗

− x

n−1

zbiega do staªej. Za dobre przybli»enie x

∗

mo»na przyj¡¢ rozwi¡-

zanie równania obliczone za pomoc¡ wywoªania funkcji octave'a:

fzero().

Zadanie 9 Zaimplementuj przybli»on¡ metod¦ Newtona, w której pochodn¡

przybli»amy ilorazem ró»nicowym, tzn.

n+1

= x

−

f (x

)

f (x

+h)−f (x

)

= x

−

f (x

) ∗ h

f (x

+ h) − f (x

)

dla ustalonego h.
Przetestuj ró»ne h, np. h = 10

−4

, 10

−7

, 10

−10

itp.

Porówna¢ zbie»no±¢ z dokªadn¡ metod¡ Newtona (szczególnie

ostatnie iteracje) dla funkcji z zadania 6.

Zadanie 10 Rozwikªywanie funkcji:

Dla funkcji y(x) zadanej w sposób uwikªany równaniem

g(x, y) = 2x

+ 3y

− 3 = 0

znajd¹ przybli»one warto±ci y

= y(x

)

dla x

= k ∗ h

dla k =

−N, . . . , N

i h = 1/N, tzn. speªniaj¡ce g(x

, y

) = 0

Testuj dla ró»nych warto±ci N, np. N = 10, 20, 40, 80, ....
oblicz y

rozwi¡zuj¡c równanie

g(x

, y

) = 0.

korzystaj¡c z odpowiedniej funkcji octave'a lub wªasnej imple-

mentacji metody Newtona.
Jak w kolejnych krokach dobiera¢ przybli»enie startowe w meto-

dzie Newtona?

Zadanie 11 Odwracanie funkcji

Rozpatrzmy dan¡ funkcj¦, np.

f (x) = sin(x) + 2 ∗ x.

Znajd¹ warto±ci funkcji odwrotnej f

−1

na odcinku [0, 5] na siatce

k ∗ h

dla k = 0, .., 100. Narysuj wykres funkcji f

−1

Sam wykres mo»na narysowa¢ du»o pro±ciej bez wyliczania war-

to±ci f

−1

. Jak to zrobi¢ w octave'ie?

Zadanie 12 Wielowymiarowa metoda Newtona

Zaimplementuj wielowymiarow¡ metod¦ Newtona w wersji z do-

kªadnym Jakobianem i w wersji, gdy Jakobian przybli»amy ró»-

nicami dzielonymi z parametrem h.
Zastosuj t¦ metod¦ do rozwi¡zania ukªadu

, x

) = x

+ 2x

= 1

, x

) = 3 ∗ x

2
1

+ x

2
2

= 1

dla ró»nych przybli»e« pocz¡tkowych.
W przypadku zbie»no±ci policz

kx − x

n−1

dla p = 1, 2, ∞ oraz q = 1, 2, 3. Czy ten iloraz zbiega do staªej

dla jakiego± z tych q?
Funkcj¦ F (x) := (f

(x), f

(x))

w octavie mo»na zdeniowa¢ na-

st¦puj¡co:

function y=F( x )

#x− wektor pionowy

y =[[1 2]∗ x −1;3.∗ x (1)∗ x(1)+x (2)∗ x (2) −1];

endfunction

Zadanie 13 Wielowymiarowa metoda iteracji prostych

Zastosuj metod¦ iteracji prostych

= x

− aF (x)

do ukªadu z zadania 12, tzn. przyjmijmy, »e F (x) := (f

(x), f

(x))

a parametr a 6= 0 np. a = 1, −1, 10

−1

itp.

Czy metoda zawsze zbiega, a je±li tak - to jak szybko? W przy-

padku zbie»no±ci policz

kx − x

n−1

dla p = 1, 2, ∞.

Zadanie 14 Do ukªadu równa« z zadania 12 zastosuj funkcj¦ fsolve () - po-

równaj z implementacj¡ metody Newtona w zadaniu 12, która

byªa zastosowana do rozwi¡zania tego ukªadu.
Czy obie metody zbiegaj¡ do tego samego rozwi¡zania dla tych

samych warto±ci wektorów startowych?
Sprawd¹ - za pomoc¡ funkcji tic i toc - czas potrzebny do roz-

wi¡zania tego ukªadu równa« za pomoc¡ obu metod.

Rozdziaª 9

Ukªady równa« liniowych -

rozkªady typu LU i LL'

W tym rozdziale zapoznamy si¦ z metodami sªu»¡cych do rozwi¡zywania

ukªadów równa« liniowych przy pomocy uzyskiwaniu odpowiednich rozkªa-

dów macierzy typu LU i LL

oraz obliczaniu macierzy odwrotnej.

Zapoznamy si¦ z odpowiednimi funkcjami octave'a sªu»¡cymi znajdowa-

niu odpowiednich rozkªadów, czy rozwi¡zywaniu ukªadów równa« liniowych

z pomoc¡ tych rozkªadów.

Podstawowe funkcje i operatory octave'a zwi¡zane z rozkªadem LU to:

•

[L,U,P]=lu(A) - funkcja zwracaj¡ca czynniki rozkªad LU nieosobliwej

macierzy A, czyli P A = LU,

•

[R]=chol(A) - funkcja zwracaj¡ca czynnik rozkªadu Choleskiego ma-

cierzy symetrycznej dodatnio okre±lonej A = A

> 0

, czyli A = R

(zazwyczaj w literaturze podaje si¦ ten rozkªad w terminach czynnika
L = R

•

operator x=A\b - operator rozwi¡zuj¡cy ukªad równa« Ax = b dla A

macierzy nieosobliwej i b wektora prawej strony,

•

inv(A) - funkcja zwracaj¡ca macierz odwrotn¡ do A macierzy nieoso-

bliwej,

•

cond(A) - funkcja zwracaj¡ca wspóªczynnik uwarunkowania macierzy
A

, czyli kAk

∗ kA

−1

Zadanie 1 Operator octave'a \ sªu»¡cy m.in. do rozwi¡zywaniu ukªadów

równa« liniowych w octavie.

Przetestuj operator octave'a \ dla ukªadu równa« z macierz¡
A = [1, 2; 3, 4]

i x = [1; 3] z f=A∗x. Czy y=A\f jest rozwi¡zaniem

tego ukªadu?
Policz normy residualn¡ kAy − fk

i norm¦ bª¦du kx − xk

dla

p = 1, 2

Powtórz testy dla macierzy osobliwej [1, 2; 3, 6] i prawie osobliwej
[1, 2; 3, 6 − 10

−6

]

. Co si¦ wtedy dzieje?

Zadanie 2 Przetestuj solver octave'a, tzn operator backslash dla dwóch

prostych ukªadów równa« liniowych z macierzami nieosobliwymi:
A

= [1, 1; 1, 0.98]

i b = [2; 1], A

= [1, 1; 1, 0.99]

i wektorem pra-

wej strony ukªadu równa« liniowych b = [2; 1] jak w poprzed-

nim zadaniu. Policz w normie drugiej ró»nic¦ rozwi¡za«, norm¦

Frobeniusa ró»nicy A

− A

oraz uwarunkowania tych macierzy

korzystaj¡c z funkcji octave'a cond().
Powtórz zadanie dla macierzy: A

= [1, 1; 1, 1 − 2]

i b = [2; 1],

= [1, 1; 1, 1 − ]

dla = 10

−p

z p = 3, 4, . . . , 16.

Zadanie 3 Funkcja inv().

Zapoznaj si¦ z pomoc¡ do funkcji: inv(). Przetestuj, obliczaj¡c

macierz odwrotn¡ do A = [1, 1; 1, −1]. Czy macierz B obliczona

za pomoc¡ tej funkcji rzeczywi±cie jest macierz¡ odwrotn¡?
Policz normy pierwsz¡ i Frobeniusa kA∗B −Ik oraz kB ∗A−Ik.
Zastosuj funkcje do rozwi¡zywania ukªadu równa« liniowych Ax =
f

dla znanego x (liczymy f = A ∗ x). Policz y jako iloczyn B z f

i policz bª¦dy residualny kAy − fk i kx − yk w normie pierwszej

i drugiej.
Powtórz to zadanie dla macierzy N ×N losowej (funkcja octave'a

rand() zwraca macierz losow¡) dla n = 10, 50, 250, 1250 ze zna-

nym rozwi¡zaniem - oszacuj czas przy pomocy funkcji tic i toc.

Porównaj czas i bª¦dy w normie pierwszej dla rozwi¡zania uzy-

skanego przy pomocy operatora octave'a \.

Zadanie 4 Funkcje lu() i chol().

•

Zapoznaj si¦ z pomoc¡ do funkcji: lu() i chol().

•

Dla macierzy A = [1, 1; 1, −1] znajd¹ jej rozkªady: P A =
LU

i rozkªad Choleskiego A = L

T
1

•

Sprawd¹ te rozkªady licz¡c normy macierzowe pierwsz¡ i

Frobeniusa bª¦dów np. P A − LU i A − L

T
1

•

Zastosuj te rozkªady do znalezienia rozwi¡zania ukªadu rów-

na« liniowych Ax = f ze znanym rozwi¡zaniem np. x =
[1; 1]

i f = [2; 0].

•

Policz normy pierwsz¡ i drug¡ wektorowe pomi¦dzy x, a

wynikiem algorytmu w, polegaj¡cym na zastosowaniu od-

powiedniego rozkªadu, oraz takie same normy residualne,

tzn. normy Aw − f.

Zadanie 5 W octavie przetestuj eliminacje Gaussa z cz¦±ciowym wyborem

i bez wyboru dla macierzy A = [e, 1; 1, 1] z e = eps/10 (funk-

cja octave'a eps zwraca epsilon maszynowy) i wektorem prawej

strony f = [1; 0]

. Trzeba tu zaprogramowa¢ samodzielnie elimi-

nacj¦ Gaussa bez wyboru elementu gªównego, tzn. bez permuta-

cji ani wierszy, ani kolumn dla macierzy 2 × 2.

Zadanie 6 Testy solvera octave'a dla macierzy Hilberta H(N)

Polecenie octave'a hilb() tworzy macierz Hilberta.

•

Stwórz macierz H(N) dla N = 10, 20, 30,

•

Dla znanego rozwi¡zania, np. staªego równego jeden na ka»-

dej pozycji, czyli sol

= 1

, policz H(N) ∗ sol = f,

•

Rozwi¡» w octavie ukªad z H(N) i f,

•

Policz normy (1,2 itd) kH(N)x − fk i kx − solk dla x roz-

wi¡zania wyliczonego przez octave,

•

Policz uwarunkowanie macierzy Hilberta dla powy»szych N.

Zadanie 7 Powtórz zadanie 6 w arytmetyce pojedynczej precyzji.

Nale»y tu sztucznie wymusi¢, za pomoc¡ funkcji octave'a single (),

aby zmienne byªy zmiennymi pojedynczej precyzji.

Zadanie 8 Przetestuj funkcj¦ invhilb() tworz¡c¡ macierz odwrotn¡ do ma-

cierzy Hilberta (por. zadanie 6).

•

Policz normy: macierzowa norm¦ Frobeniusa i norm¦ macie-

rzow¡ pierwsz¡ ró»nicy macierzy E = H(N)∗iH(N)−I dla
N = 10, 20, 30

, gdzie iH(N) macierz odwrotna do macierzy

Hilberta obliczona przez invhilb().

•

Wykorzystaj iH(N) do rozwi¡zania ukªadu równa« z ma-

cierz¡ Hilberta H(N), tzn.:

(a) stwórz macierze H(N) oraz iH(N) dla N = 10, 20, 30

korzystaj¡c z funkcji hilb() i invhilb(),

(b) dla znanego rozwi¡zania sol, policz H(N) ∗ sol = f,

x = iH(N ) ∗ f

(d) policz normy kH(N)x−fk

i kx−solk

dla p = 1, 2, ∞.

Zadanie 9 Stwórz w octavie macierz trójdiagonaln¡ A = (a

)

n
i,j=1

wymiaru

n × n

z 2 na diagonali i −1 na pod- i nad diagonal¡, tzn.

i,j







i = j

−1

|i − j| = 1

|i − j| > 1

i, j = 1, . . . , n.

przy pomocy funkcji octave'a diag(), jak i wprost w p¦tli. Po-

równaj czas u»ywaj¡c tic i toc dla n = 10, 100, 1000.
Policz uwarunkowanie macierzy dla ró»nych N.
Policz macierz odwrotn¡ przy pomocy inv() (czy jest trójdiago-

nalna?).

Zadanie 10 Sprawd¹ z wykorzystaniem funkcji octave'a chol(), czy macierz

z zadania 9 jest symetryczna i dodatnio okre±lona.

Zadanie 11 Zaimplementuj rozkªad LU dla macierzy trójdiagonalnej N ×

bez wyboru elementu gªównego, tzn. stwórz wªasn¡ funkcj¦

octave'a:
function [x,d,l]=lu3diag(a,b,c, f ,N)
Parametry funkcji:

• a, b, c

- przek¡tna, podprzek¡tna i nadprzek¡tna macierzy A,

• f

- wektor prawej strony,

• N

- wymiar zadania - dªugo±¢ przek¡tnej a.

Funkcja powinna zwróci¢ x - rozwi¡zanie Ax = f oraz czynniki

rozkªadu macierzy A = LU, czyli diagonal¦ macierzy górnotrój-

k¡tnej U w wektorze d i poddiagonal¦ L - macierzy dolnotrójk¡t-

nej, trójdiagonalnej z jedynkami na diagonali, czyli wektor l.
Naddiagonala U równa si¦ nadiagonali A, tzn. wektorowi c.

Zadanie 12 Zastosuj funkcj¦ z zadania 11 do uzyskania rozkªadu LU macierzy

z zadania 9 dla N = 4, 10, 15. Porównaj z wynikiem uzyskanym

przy pomocy funkcji octave'a lu().

Zadanie 13 Zastosuj funkcj¦ z zadania 11 do rozwi¡zania ukªadu równa« z

macierz¡ z zadania 9 dla N = 10, 100, 1000, czy 10000:

•

Rozwi¡» ukªad dla znanego rozwi¡zania, np. x = (1, . . . , 1)

•

Porównaj czas rozwi¡zywania wªasnym algorytmem i algo-

rytmem octave'a, czyli operatorem A\f,

•

Policz bª¡d rezydualny i bª¡d rzeczywisty w normie drugiej,

tzn. kx − ˜xk

i kf − A ∗ ˜xk

, gdzie f = A ∗ x, a ˜x to

rozwi¡zanie obliczone przez octave'a lub obliczone wªasnym

algorytmem.

Zadanie 14 Zaimplementuj rozkªad LU dla macierzy trójdiagonalnej z cz¦-

±ciowym wyborem elementu gªównego.
Nast¦pnie testuj ten rozkªad jak w zadaniu 13.
Zastosuj ten rozkªad do rozwi¡zania ukªadu równa« z macierz¡ z

poprzedniego zadania A i z macierz¡ A − 1.5 ∗ I.

Zadanie 15 Zaprogramuj rozkªad Choleskiego typu dla macierzy A = A

> 0

trójdiagonalnej, tzn. policz L dolnotrójk¡tn¡ z jedn¡ poddia-

gonal¡ (czyli reprezentowan¡ przez dwa wektory z przek¡tn¡ i

podprzek¡tn¡) tak¡, »e A = L

Zastosuj do rozwi¡zania ukªadu równa« z zadania 9.
Testuj analogicznie do zadania 13.
Porównaj macierz L uzyskan¡ w ten sposób z macierz¡ uzyskan¡

przez chol() zastosowan¡ do tej samej macierzy.

Rozdziaª 10

LZNK. Rozkªad QR. Metoda

Householdera

W tym rozdziale zajmiemy si¦ liniowym zadaniem najmniejszych kwadratów

(LZNK).

Dla danej macierzy A wymiaru M × N i wektora b wymiaru M chcemy

znale¹¢ wektor x wymiaru N taki, »e

kAx − bk

= min

kAy − bk

Je±li A jest macierz¡ kolumnami regularn¡ (rz¡d A jest maksymalny równy
N

), to zadanie to ma jednoznaczne rozwi¡zanie i nazywamy je regularnym

LZNK (RLZNK).

Podstawowym algorytmem sªu»¡cym jego rozwi¡zaniu jest metoda Ho-

useholdera, czyli znalezienia rozkªadu QR macierzy A, gdzie Q to macierz

ortogonalna - iloczyn macierzy Householdera, a R to macierz górnotrójk¡tna.

W tym rozdziale przetestujemy podstawowy operator octave'a sªu»acy do

rozwi¡zywania dowolnego LZNK, tzn. operator \. Zauwa»my, »e je±li M = N

i A jest macierz¡ kolumnami regularn¡, to A jest nieosobliwa i to regularne

LZNK jest równowa»ne rozwi¡zaniu ukªadu równa« liniowych Ax = b.

Przetestujemy w zadaniach funkcj¦ octave'a qr() sªu»¡c¡ znajdowaniu

rozkªadu QR macierzy, kilka wªasno±ci macierzy (przeksztaªce«) Househol-

dera i rozwi¡»emy kilka konkretnych LZNK.

Zadanie 1 Operator octave'a \ sªu»¡cy m.in. do rozwi¡zywaniu ukªadów

równa« liniowych i LZNK w octavie.

•

Przetestuj operator octave'a \ rozwi¡zuj¡c RLZNK dla ma-

cierzy A ze znanym konkretnym rozwi¡zaniem x:

A = [1, 1; 1, −1; 1, 3],

x = [1; 2]

przyjmuj¡c, »e f = Ax.
Czy y=A\f jest rozwi¡zaniem tego ukªadu?
Policz normy residualn¡ kAy − fk

i norm¦ bª¦du kx − yk

•

Przetestuj ten operator dla nieregularnego LZNK dla ma-

cierzy A

z x = [1; 1; 1] i f = Ax.

Czy y=A'\f jest rozwi¡zaniem tego ukªadu?
Policz normy residualn¡ kA

y − f k

i norm¦ bª¦du kx−yk

Zadanie 2 Rozkªad QR w octave. Funkcje qr().

(a) Zapoznaj si¦ z pomoc¡ do funkcji: qr().

(b) Dla macierzy A = [1, 1; 1, −1; 1, 3] znajd¹ jej rozkªad A =

z pomoc¡ funkcji qr().

Frobeniusa QQ

− I

i Q

Q − I

(d) Sprawd¹ ten rozkªad licz¡c normy macierzowe: norm¦ drug¡

i norm¦ Frobeniusa bª¦du A − QR.

(e) Zastosuj ten rozkªad do znalezienia rozwi¡zania LZNK Ax =

ze znanym rozwi¡zaniem, np. x = [1; 0] i f = [1; 1; 1].

(f) Policz norm¦ drug¡ wektorow¡ pomi¦dzy x, a wynikiem al-

gorytmu w, polegaj¡cym na zastosowaniu odpowiedniego

rozkªadu oraz takie same normy residualne, tzn. normy dru-

gie Aw − f oraz Rw − Q

Zadanie 3 Ukªad równa« normalnych, a rozkªad QR.

Rozpatrzmy macierz A

2n,k

- pod-macierz wymiaru 2n × k ma-

cierzy Vandermonde'a A

2n,2n

dla 2n w¦zªów równoodlegªych na

[0, 1]

LZNK z A

2n,k

z wektorem prawej strony f równym pierwszej

kolumnie tej macierzy (rozwi¡zanie to pierwszy wersor) rozwi¡»

trzema sposobami:

(a) u»ywaj¡c operator \,

(b) u»ywaj¡c rozkªad QR uzyskanym funkcj¡ qr(),

Bx = g

dla

B = A

T
2n,k

2n,k

g = A

T
2n,k

tzn. tworzymy macierz ukªadu równa« normalnych B, wek-

tor prawej strony g ukªadu równa« normalnych, a nast¦pnie

rozwi¡zujemy ukªad równa« normalnych operatorem \ .

Macierz Vandermonde'a mo»na w octave'ie utworzy¢ za pomoc¡

funkcji vander().
Przeprowad¹ testy dla N = 10, 20, 40, 80 i k = 2, 4, n. Porównaj

•

czas oblicze«

•

bª¡d - kx − yk

•

bª¡d residualny kAx − fk

dla x rozwi¡zania dokªadnego LZNK, f wektora prawej strony

LZNK, y przybli»enia rozwi¡zania uzyskanego dan¡ metod¡.

Zadanie 4 Rozkªad QR a operator \ przy rozwi¡zywaniu ukªadów równa«

liniowych,
Rozpatrzmy macierz A

n,n

Vandermonde'a dla n w¦zªów równo-

odlegªych na [0, 1].
Ukªad równa« liniowych z wektorem prawej strony równym pierw-

szej kolumnie tej macierzy (rozwi¡zanie to pierwszy wersor) roz-

wi¡» dwoma sposobami:

(a) operatorem \,

(b) rozkªadem QR uzyskanym funkcj¡ qr(),

Przetestuj dla N = 10, 20, 40, 80. Porównaj

•

czas oblicze«,

•

bª¡d kx − yk

•

bª¡d rezydualny: kAx − fk

dla x rozwi¡zania dokªadnego tego ukªadu równa«, f wektora

prawej strony ukªadu i y przybli»enia rozwi¡zania uzyskanego

dan¡ metod¡.

Zadanie 5 Krzywa najlepiej pasuj¡ca do danych punktów.

Zastosuj octave'a do znalezienia wspóªczynników a, b krzywej naj-

lepiej pasuj¡cej do zadanych punktów: (x

, y

)

, tzn. znajd¹ takie

a, b

, »e

|ax

2
k

+ by

− 1|

= min

c,d

|cx

2
k

+ dy

− 1|

Za (x

, y

)

przyjmiemy zaburzone punkty z danej elipsy

1 − 4 ∗ x

2
k

+ z

gdzie z

to zaburzenie wylosowane z [0, 10

−2

]

a x

= 1/k

lub

= −1 + h ∗ k

dla h = 2/N k = 1, ..., N.

Czy obliczone a i b jest bliskie 4 i 1?
W jednym oknie zaznacz punkty (x

, y

)

plusami oraz narysuj

fragmenty wykresów obu elips: pierwszej - dla a = 4, b = 1 i

drugiej elipsy - dopasowanej do zaburzonych punktów.
Powtórz obliczenia dla ró»nych zaburze« z

Zadanie 6 Zaprogramuj funkcj¦ octave'a

function y=H(x ,w, nw)

która dla danych wektorów ~x i ~w tego samego wymiaru N i ska-

laru nw = kwk

zwróci wektor y = H

dla

= I − 2 ∗

w ~

czyli przeksztaªcenia (macierzy) Householdera.
Skalar mo»e by¢ parametrem opcjonalnym. Je±li funkcja b¦dzie

wywoªana z dwoma tylko parametrami, to norm¦ w mo»na obli-

czy¢ w tej funkcji.
Przetestuj t¦ funkcj¦ dla losowych wektorów ~x i ~w i sprawd¹, czy

= k~

x) = x.

Zadanie 7 Napisz ogólniejsz¡ wersj¦ funkcji z poprzedniego zadania, tzn.:

funkcj¦:

function Y=Hm(X,w, nw) ,

gdzie X macierz N × M i wtedy zwracany wynik to macierz
Y = H

∗ X

. Pozostaªe dwa parametry funkcji pozostan¡ bez

zmian.
Czy mo»na zaimplementowa¢ tak¡ funkcj¦ bez u»ycia p¦tli?
Sprawd¹, wykorzystuj¡c t¦ funkcj¦, czy mno»enie przez macierz

Householder nie zmienia norm macierzowych drugiej i Frobe-

niusa, tzn. czy:

kAk

= kH

∗ Ak

= kA ∗ H

kAk

= kH

∗ Ak

= kA ∗ H

dla losowej macierzy A i H

macierzy Householdera dla losowego

wektora w 6= 0.

Zadanie 8 Znajd¹ wektor Householdera ~w taki, »e odpowiednie przeksztaª-

cenie Householdera przeprowadza dany wektor ~u 6= 0 na kierunek

drugiego danego niezerowego wektora l ∗ ~v 6= 0 dla l =

kvk

kuk

Przetestuj dla dowolnych dwóch ró»nych wektorów o tej samej

dªugo±ci, czy rzeczywi±cie H

u = ~

Zadanie 9 Zastosuj metod¦ Householdera do rozwi¡zania zadania znalezie-

nia prostej y = ax+b najlepiej przybli»aj¡cej N punktów (x

, y

) =

(k, 1 + 2 ∗ k +

)

dla k = 1, . . . , N gdzie (

, . . . ,

)

to losowy

wektor za zakresu [−, ], tzn.:

|ax

+ b − y

= min

c,d

|cx

+ d − y

Nale»y testowa¢ dla warto±ci = [1, 10

−1

, 10

−2

, 10

−3

]

Funkcja rand(n) generuje wektor losowy o rozkªadzie jednostaj-

nym na [0, 1] w octavie.
Porównaj z wynikami otrzymanymi za pomoc¡ standardowej funk-

cji octave'a, tzn. \ , oraz przy wykorzystaniu funkcji octave'a

qr(A).

Zadanie 10 Zaprogramuj metod¦ Householdera rozwi¡zywania ukªadu rów-

na« liniowych A~x = ~b dla A macierzy trójdiagonalnej N × N,

tzn. napisz funkcj¦ octave'a:
function [x]=hous3diag(a,b,c,f,N)
Parametry funkcji:

• a, b, c

przek¡tna, pod-przek¡tna i nad-przek¡tna macierzy A,

• f

- wektor prawej strony,

• N

- wymiar zadania - dªugo±¢ przek¡tnej a.

Funkcja zwraca x rozwi¡zanie Ax = f.
Przetestuj dziaªanie funkcji analogicznie do zadania 4 dla macie-

rzy trójdiagonalnej o staªych diagonalach, np. takiej, »e elementy

na gªównej diagonali s¡ równe dwa, a elementy na pod- i nad-

diagonalach s¡ równe minus jeden dla N = 10

z p = 1, 2, . . . , 9..

Za wektor prawej strony f mo»emy przyj¡¢ pierwsz¡ kolumn¦

macierzy A.

Rozdziaª 11

Numeryczne zadanie wªasne

W tym rozdziale zajmiemy si¦ symetrycznym zadaniem wªasnym, tzn. za-

daniem znajdowania warto±ci i/lub wektorów wªasnych dla macierzy syme-

trycznej A = A

. W zadaniach zbadamy, jak dziaªa podstawowa funkcja

octave'a rozwi¡zuj¡ca zadanie wªasne, tzn. funkcja octave'a eig() oraz zba-

damy zbie»no±¢ ró»nych wersji metody pot¦gowej.

Metoda pot¦gowa jest zdeniowana nast¦puj¡co: dla dowolnego wektora

= x

o normie drugiej równej jeden (np. losowego), deniujemy ci¡g

iteracji metody pot¦gowej nast¦puj¡co:

= Ax

k−1

(iteracja)

kˆ

(normowanie)

= x

T
k

(iloraz

Rayleigha)

k > 0

(11.1)

Zgodnie z teori¡, o ile x

nie jest ortogonalny do wektora wªasnego (pod-

przestrzeni wªasnej) dla najwi¦kszej warto±ci wªasnej, to x

zbiegnie do tego

wektora wªasnego, a r

- iloraz Rayleigha - do tej warto±ci wªasnej.

Zadanie 1 Zapoznaj si¦ z pomoc¡ do funkcji octave'a eig().

Za jej pomoc¡ oblicz warto±ci wªasne macierzy





0 1 1
1 0 1
1 1 0





Sprawd¹, czy dla λ warto±ci wªasnej A prawd¡ jest, »e wyznacznik
A − λ ∗ I

wynosi zero.

Wyznacz macierz V tak¡, »e jej kolumny to wektory wªasne A.
Przetestuj w normie macierzowej Frobeniusa, czy kAV − ΛV k

wynosi zero. Tutaj Λ - to macierz diagonalna z warto±ciami wªa-

snymi A.

Zadanie 2 Zastosuj metod¦ pot¦gow¡ (11.1) do zadania wªasnego z macie-

rz¡ A z poprzedniego zadania. Za warunek stopu przyjmiemy, »e

iloraz Rayleigha r

praktycznie przestanie si¦ zmienia¢, co mo»e

sugerowa¢, »e jest bliski granicy, czyli odpowiedniej warto±ci wªa-

snej, tzn. je±li

k+1

− r

| < tol

to zatrzymujemy iteracj¦. Przyjmijmy, »e tol = 10

−12

Przetestuj, czy otrzymane r

rzeczywi±cie jest dobrym przybli»e-

niem najwi¦kszej co do moduªu warto±ci wªasnej A - czyli dwa, a
x

zbiegªo do wektora wªasnego dla tej warto±ci wªasnej o normie

drugiej równej jeden, tzn. do v

= (

√

)

, lub −v

Przetestuj algorytm bior¡c za x

(a) wektory losowe - otrzymane za pomoc¡ funkcji losowych

rand() lub randn().

(b) wektor wªasny dla drugiej warto±ci wªasnej, tzn. dla minus

jeden, czyli unormowany wektor v

= (−2, 1, 1)

= (

√

)

, np. wek-

tor otrzymany przez zaburzenie ostatniej skªadowej wektora
v

: (−2, 1, 1.01)

Czy we wszystkich tych przypadkach metoda zbiegnie do v

lub

−v

Czy wybór x

wpªywa na ilo±¢ iteracji metody pot¦gowej?

Zadanie 3 Napisz funkcj¦ octave'a z implementacj¡ metody pot¦gowej (11.1),

tzn. napisz m-plik z funkcj¡:

function [ x , r , i t ]=power (A, tol , x0 )

która za parametry ma

• A

- macierz wymiary N × N, której wektora i warto±ci wªa-

snej szukamy

• tol

- warunek stopu, domy±lnie przyjmuj¡cy warto±¢ 10

−12

• x

-wektor startowy niezerowy, domy±lnie losowy,

i która zwraca

•

przybli»enie - unormowanego w normie drugiej - wektora

wªasnego x dla najwi¦kszej warto±ci wªasnej co do moduªu

macierzy A,

• r

- przybli»enie tej warto±ci wªasnej,

• it

- ilo±¢ iteracji metody pot¦gowej.

Warunek stopu przyjmiemy jak w poprzednim zadaniu, tzn. »e

iloraz Rayleigha r

przestanie si¦ zmienia¢. Je±li

k+1

− r

| < tol

to zatrzymujemy iteracj¦.

Zadanie 4 Utwórz macierz symetryczn¡ (ale nie diagonaln¡) N × N o war-

to±ciach wªasnych 1, . . . , N. Sprawd¹ dla N = 10, 100, 1000 czy

funkcja octave'a eig() znajdzie warto±ci wªasne tej macierzy.

Zadanie 5 Utwórz macierze symetryczne A

(niediagonalne) wymiaru 11 ×

o warto±ciach wªasnych 1, . . . , 10, a dla a = 11, 101, 1001 i ze

znanymi wektorami wªasnymi.
Sprawd¹, ile iteracji b¦dzie potrzebowaªa metoda pot¦gowa dla

losowego wektora startowego do znalezienia najwi¦kszej warto±ci

wªasnej.
Czy potwierdzaj¡ si¦ wyniki teoretyczne, »e szybko±¢ zbie»no±ci

zale»y od stosunku najwi¦kszej co do moduªu warto±ci wªasnej

do drugiej co do moduªu warto±ci wªasnej, czyli w tym zadaniu
a/10

Zadanie 6 Dla macierzy z zadania 5 dla a = 11, 101 ze znanym wektorem

wªasnym v

dla warto±ci wªasnej a, przetestuj zbie»no±¢ x

obli-

czanego przy pomocy metody pot¦gowej do v

. W szczególno±ci

sprawd¹, czy

− x

= O(|λ

/λ

)

dla λ

= 10

i λ

= a

Zadanie 7 Dla macierzy z zadania 5 dla a = 11, 101 dla warto±ci wªasnej

, przetestuj zbie»no±¢ r

ilorazu Rayleigha obliczanego przy po-

mocy metody pot¦gowej do λ

. W szczególno±ci sprawd¹, czy

|λ

− r

| = O(|λ

/λ

)

dla λ

= 10

i λ

= a

Zadanie 8 Dla macierzy z zadania 5 dla a = 11, 101 ze znanym wektorem

wªasnym v

dla warto±ci wªasnej 11 przetestuj zbie»no±¢ metody

pot¦gowej dla ró»nych wektorów startowych:

•

dla losowego wektora x

ortogonalnego do v

(jak go mo»na

utworzy¢?)

•

dla wektora prawie ortogonalnego - za x

bierzemy w+10

−6

∗

, gdzie w

= 0

i z to wektor losowy o normie drugiej

równej jeden.

•

dla dowolnego losowego wektora x

Prosz¦ powtórzy¢ testy kilka razy dla ka»dego podpunktu - czy

zdarzyªo si¦, »e wyniki byªy istotnie ró»ne w zale»no±ci od wylo-

sowanego wektora?

Zadanie 9 Zaimplementuj funkcj¦ analogiczn¡ do tej z zadania 3 z imple-

mentacj¡ odwrotnej metody pot¦gowej, tzn. metody pot¦gowej

zastosowanej do macierzy odwrotnej A

−1

zamiast A:

function [ x , r , i t ]= invpower (A, tol , x0 )

Wszystkie parametry i zwracane warto±ci b¦d¡ takie same jak w

funkcji power() z zadania 3.

Zadanie 10 Utwórz macierz z zadania 5 dla a = 11, 101, 1001. Sprawd¹, ile

iteracji b¦dzie potrzebowaªa odwrotna metoda pot¦gowa dla lo-

sowego wektora startowego do znalezienia przybli»enia najmniej-

szej warto±ci wªasnej przy ustalonym warunku stopu |r

−r

k+1

| ≤

−12

. Tutaj r

- to odpowiedni iloraz Rayleigha.

Zadanie 11 Zaimplementuj funkcj¦ analogiczn¡ do tej z zadania 3 z imple-

mentacj¡ odwrotnej metody pot¦gowej z przesuni¦ciem (shift),

tzn. metody pot¦gowej zastosowanej do macierzy (A − b ∗ I)

−1

dla parametru b:

function [ x , r , i t ]= invpowwiel (A, b , tol , x0 )

Tutaj b b¦dzie dodatkowym parametrem tej funkcji, a pozostaªe

parametry i zwracane warto±ci pozostaj¡ takie same jak w funkcji

z zadania 3.

Zadanie 12 Utwórz macierz z zadania 5 dla a = 11, 101, 1001.

Sprawd¹, ile iteracji b¦dzie potrzebowaªa odwrotna metoda pot¦-

gowa z przesuni¦ciem dla losowego wektora startowego i ró»nych

warto±ci parametru b = 1.3, 1.5, 1.8, 2.1, 2.01 do znalezienia ró»-

nych warto±ci wªasnych macierzy A

Zadanie 13 Zaimplementuj funkcj¦ analogiczn¡ do tej z zadania 3 z imple-

mentacj¡ metody ilorazów Rayleigha:

function [ x , r , i t ]= raylegh (A, b , tol , x0 )

Metoda polega na tym, »e dla danego wektora startowego x

drugiej normie równej jeden, obliczamy iloraz Rayleigha r

T
0

i iterujemy:

(A−r

k−1

−1

k−1

k(A−r

k−1

−1

k−1

= x

T
k

k > 0

(11.2)

Oczywi±cie w implementacji prosz¦ nie oblicza¢ macierzy odwrot-

nej do A − r

∗ I

, tylko rozwi¡zywa¢ ukªad równa« liniowych z t¡

macierz¡.
Wszystkie parametry i zwracane warto±ci b¦d¡ takie same jak w

funkcji power() z zadania 3.

Zadanie 14 Przetestuj metod¦ ilorazów Rayleigha na przykªadzie macierzy

symetrycznej z zadania 5 dla a = 11, 101, 1001 bior¡c losowy

wektor startowy. W kolejnych testach prosz¦ wzi¡¢ wektory star-

towe bliskie wektorom wªasnym dla warto±ci wªasnych 2 i 7.

Zadanie 15 Zaimplementuj funkcj¦ octave'a sprowadzaj¡c¡ macierz syme-

tryczn¡ A wymiaru N ×N do macierzy podobnej trójdiagonalnej:
Q

BQ = A

lub B = QAQ

za pomoc¡ odpowiednich przeksztaª-

ce« Householdera:

function [B,Q]= sym2tridiag (A, typ )

która za parametry ma

• A

- symetryczn¡ macierz wymiary N × N

• tol

- typ zwracanej macierzy Q, typ=0 oznacza, »e Q jest

zwracana jako macierz N × N - wymno»one macierze Ho-

useholdera (domy±lna opcja), typ=1 oznacza, »e w macierzy
Q

zwracamy kolumny odpowiednich macierzy Householdera

i która zwraca

• B

- trójdiagonaln¡ macierz symetryczn¡ podobn¡ do A w

formacie rzadkim octave'a (sparse),

• Q

macierz b¦d¡c¡ iloczynem przeksztaªce« Householdera

lub macierz, w której kolumnach s¡ odpowiednie wektory

Householdera w zale»no±ci od warto±ci parametru typ.

Zadanie 16 Przetestuj funkcj¦ z poprzedniego zadania dla jakiej± macierzy

symetrycznej A (ani diagonalnej, ani trójdiagonalnej, ani trójk¡t-

nej) o znanych warto±ciach wªasnych, np. dla macierzy z zadania

5 dla a = 11, 101, 1001.
Policz norm¦ Frobeniusa kA − Q

BQk

i kQAQ

− Bk

Sprawd¹, czy funkcja octave'a eig() dla obu macierzy zwróci te

same warto±ci wªasne, i czy odpowiednie wektory wªasne speª-

niaj¡:

= w

= v

gdzie w

to wektor wªasny macierzy B dla warto±ci wªasnej λ

a v

- to wektor wªasny dla tej samej warto±ci wªasnej λ

dla

macierzy A.

Rozdziaª 12

Algorytm FFT

W tym rozdziale przedstawimy zadania laboratoryjne, w których przetestu-

jemy dyskretn¡ transformat¦ Fouriera (DFT) oraz szybki algorytm jej obli-

czania, czyli algorytm szybkiej trasformaty Fouriera; FFT (ang. Fast Fourier

Transform)

W opisie matematycznym przyjmujemy, »e mamy do czynienia z wekto-

rami okresowymi o okresie dªugo±ci N indeksowanymi od zera i x

= x

k+N

dla k ujemnych.

Przez F

oznaczymy wynik dziaªania operatora DFT na wektorze x, a

przez F

−1

oznaczymy wynik dziaªania operatora odwrotnego do DFT na

wektorze x.

Zadanie 1 Funkcja octave'a z algorytmem FFT t (), która oblicza warto±¢

DFT i funkcja octave'a z algorytmem odwrotnym FFT it (),

która oblicza warto±¢ transformaty odwrotnej do DFT.
Sprawd¹, czy rzeczywi±cie operacje wykonywane przez t () i

it () s¡ do siebie odwrotne, tzn. dla ró»nych losowych wektorów

o dªugo±ci N = 4, 8, 16, 32 policz y = F

za pomoc¡ t (), a

nast¦pnie z = F

−1

za pomoc¡ it ().

Policz normy drugie ró»nicy z − x.

Zadanie 2 Skalowanie w funkcjach octave'a t () i it ().

W literaturze zdarza si¦ denicja DFT bez 1/N przed macierz¡.

Wtedy macierz odwrotn¡ do DFT nale»y przemno»y¢ przez N

albo obie przez 1/

√

Sprawd¹, jaka jest staªa przed DFT obliczanym przez funkcj¦

t (), tzn. wyznacz C

takie, »e

= C

N −1

k=0

exp(−2π ∗ j ∗ k/N )

dla F

obliczanego przez t¦ funkcj¦ octave'a .

Zadanie 3 Sprawd¹, czy DFT obliczona przy pomocy funkcji t () speªnia

= C

kxk

ze staª¡ C

dla stu losowych wektorów x, oraz czy analogicznie

odwrotna DFT obliczona it () speªnia:

−1

= C

−1

kxk

Zadanie 4 Algorytm szybkiego mno»enia wielomianów.

Napisz funkcj¦ octave'a

function w=polymult (p , q )

szybko mno»¡c¡ dwa wielomiany stopnia nie wi¦kszego od N,

dla których znamy ich wspóªczynniki w bazie pot¦gowej. Funk-

cja powinna korzysta¢ z DFT i odwrotnej DFT, czyli dyskretnej

transformaty Fouriera i odwrotnej dyskretnej transformaty Fo-

uriera.
Parametry funkcji powinny by¢:

•

wektor p = (p

)

ze wspóªczynnikami wielomianu P (x) =

•

wektor q = (q

)

ze wspóªczynnikami wielomianu Q(x) =

Funkcja powinna zwróci¢ wspóªczynniki wielomianu W (x) = P (x)∗
Q(x) =

Sprawd¹ dziaªanie funkcji dla wielomianów P (x) = 1+x i Q(x) =
2 + x

, czyli (P ∗ Q)(x) = 2 + 3x + x

, a potem dla wielomianów

du»ych stopni, czyli np. dla Q(x) = x

i P (X) = x

+ 1

funkcja

powinna zwróci¢ wspóªczynniki wielomianu (P ∗ Q)(x) = x

100

Porównaj czas oblicze« tej funkcji z wynikiem funkcji octave'a

conv() dla wielomianów ró»nych stopni z losowymi wspóªczyn-

nikami.

Zadanie 5 DFT a dyskretny splot dwóch wektorów.

Napisz dwie funkcje obliczaj¡ce splot dwóch wektorów z = x ? y

•

wprost z denicji, tzn. z

N −1
j=0

k−j

•

z wykorzystaniem DFT - tu musimy skorzysta¢ z tego, »e

(x ? y))

= α

∗ (F

k = 0, . . . , N − 1

dla pewnej staªej α

. Trzeba t¦ staª¡ wyznaczy¢ teoretycz-

nie albo eksperymentalnie.
Przetestuj obie funkcje dla losowych wektorów i ró»nych
N = 4, 8, 64, 1024

. Porównaj z wynikami funkcji octave'a

conv(). Nale»y zapozna¢ si¦ z pomoc¡ do tej funkcji.

Zadanie 6 Filtry a DFT.

W niektórych zastosowaniach wspóªczynniki wektora x mog¡ od-

powiada¢ próbkom sygnaªu (np. d¹wi¦ku). Stosuje si¦ wtedy

ltry w postaci splotu x z zadanym wektorem F peªni¡cym rol¦

ltru: x ? F .
Rozwa»my ltr postaci F =

1
4

∗ (1, 1, 1, 1, 0, . . . , 0)

(pierwsze

cztery wspóªrz¦dne s¡ równe

1
4

, pozostaªe - zero).

Zastosuj ten ltr do wektora imituj¡cego próbki sygnaªu z szu-

mem: z = x + y, gdzie x

= sin(k ∗ h)

dla h = 2 ∗ π/N - imituje

sygnaª bez szumów, a y jest wektorem losowym o warto±ciach

w [−0.1, 0.1] imituj¡cym losowe zaburzenie, czyli tzw. szumy.

Przetestuj ltr dla du»ych N, np. N = 1024 lub 2048, tzn.:

•

Narysuj wykres z przed zastosowaniem ltra i po zastosowa-

niu, tzn. wykresy z i z?F . Mo»na narysowa¢ maªy fragment

wykresu, np dla 100 ≤ k ≤ 120. Jaki efekt optyczny wida¢

na wykresach?

•

Policz maksimum z moduªu ró»nicy kolejnych elementów

obu wektorów, tzn. max

k>0

− z

k−1

i max

k>0

|(z ? F )

−

(z ? F )

k−1

. Im ta warto±¢ dla danego wektora z próbkami

sygnaªu jest mniejsza, tym sygnaª jest gªadszy, czyli mo»emy

uzna¢ go za sygnaª z mniejsz¡ ilo±ci¡ szumów.

Projekty zaliczeniowe

Rozdziaª 13

Przykªadowe projekty

zaliczeniowe

W tej cz¦±ci skryptu przedstawimy przykªady projektów na zaliczenia zaj¦¢

z laboratorium komputerowego z matematyki obliczeniowej. Projekty mo»na

potraktowa¢ jako trudniejsze zadania laboratoryjne.

Niektóre projekty s¡ proste je±li chodzi o stron¦ programistyczn¡. W

przypadku wi¦kszo±ci takich projektów gªówn¡ cz¦±ci¡ jest przetestowanie

danej metody. Inne projekty mog¡ polega¢ na zaimplementowaniu bardziej

skomplikowanej metody z wykorzystaniem ró»nych funkcji octave'a.

Interpolacja wielomianowa. Algorytm ró»nic dzie-

lonych.

1. Zaprogramuj w octave funkcje obliczaj¡c¡ warto±ci wielomianu zada-

nego w bazie Newtona dla danych w¦zªów zmodykowanym algoryt-

mem Hornera dla danej tablicy punktów. Parametrami maj¡ by¢:

• x

tablica punktów - wektor dªugo±ci N z w¦zªami bazy Newtona,

• a

wektor dªugo±ci N + 1 ze wspóªczynnikami wielomianu w bazie

Newtona,

• N

- stopie« wielomianu (mo»na opcjonalnie obliczy¢ N z wektora

wspóªczynników).

Funkcja zwraca y tablice warto±ci wielomianu w punktach x.

2. Zaprogramuj funkcj¦ obliczaj¡c¡ ró»nice dzielone dla danych N + 1

w¦zªów i warto±ci funkcji w tych w¦zªach algorytmem ró»nic dzielonych

(jak z tabelki ró»nic dzielonych). Parametry funkcji to:

• x

wektor dªugo±ci N + 1 z ró»nymi w¦zªami

• y

wektor dªugo±ci N + 1 z warto±ciami funkcji f w w¦zªach.

Funkcja powinna zwróci¢ wektor rd z ró»nicami dzielonymi rd(k) =
f [x

, . . . , x

]

dla k = 0, . . . , N.

Czy mo»na to zrobi¢ za pomoc¡ tylko jednej p¦tli w octavie?

3. Testy: Interpolacja funkcji f(x) = sin(x) i g(x) = sin(4∗x) na [−π, 2π]

dla w¦zªów równoodlegªych i Czebyszewa.

•

Porównanie algorytmu ró»nic dzielonych z algorytmem z funkcji

octave'a polyt.
Znajd¹ wielomiany interpolacyjne w w¦zªach Czebyszewa i równo-

odlegªych stopnia N dla N = 4, 8, 16, 32, 64 w bazie Newtona za

pomoc¡ funkcji z projektu i za pomoc¡ funkcji polyt octave'a

w bazie pot¦gowej. Oblicz dyskretn¡ norm¦ maksimum (na 1000

punktach) mi¦dzy wielomianem w bazie Newtona uzyskanym wªa-

sn¡ funkcj¡, a jej wielomianem interpolacyjnym uzyskanym polyt.

Czy ró»nice s¡ pomijalne?

•

Bª¡d interpolacji

Oblicz dyskretn¡ norm¦ maksimum (na 1000 punktach) mi¦dzy

funkcj¡, a jej wielomianem interpolacyjnym otrzymanym za po-

moc¡ funkcji z projektu dla N = 4, 8, 16, 32, 64. Narysuj wykresy

funkcji i wykres jej wielomianu interpolacyjnego (na jednym ry-

sunku).

Warto±¢ wielomianu w bazie Newtona obliczamy algorytmem Hornera

z pierwszego podpunktu.

Równania nieliniowe. Metoda Steensena

Zaprogramuj funkcj¦ octave'a z metod¡ Steensena zdeniowan¡ wzorem

n+1

= x

−

f (x

)

f (x

+ f (x

)) − f (x

)

- która ma znale¹¢ przybli»enie f(x

∗

) = 0

Napisz funkcj¦:

function [ x , i t ]= s t e f (FCN, x0 )

w m-pliku steff.m z parametrami:

•

FCN - 'wska¹nik' do funkcji (function handle) obliczaj¡cej warto±¢ f(x)

dla danego argumentu x,

• x0

- liczba, przybli»enie startowe.

Funkcja powinna zwraca¢

• x

obliczone przybli»enie pierwiastka,

• it

- ilo±¢ iteracji.

Funkcja powinna si¦ zatrzyma¢, drukuj¡c komunikat o braku zbie»no±ci

na ekranie komputera, gdy ilo±¢ iteracji przekroczy 100.

W przypadku przekroczenia ilo±ci iteracji funkcja ma zwróci¢ komunikat

o tym na ekran.

Przetestuj metod¦ Steensena z wykorzystaniem funkcji ste () na przy-

kªadach równa« rozpatrywanych w zadaniach w 8.

W szczególno±ci nale»y sprawdzi¢, czy metoda zbiega kwadratowo lokal-

nie, o ile f

∗

) 6= 0

dla konkretnej funkcji f(x) = x

− 27

Równania nieliniowe. Metoda Halleya

Zaprogramuj funkcj¦ octave'a z metod¡ Halleya, która jest metod¡ New-

tona zastosowan¡ do funkcji

g(x) :=

f (x)

p|f

(x)|

która ma znale¹¢ przybli»enie f(x

∗

) = 0

Wyprowad¹ wzór na kolejn¡ iteracj¦ metody Halleya, w którym b¦dziemy

potrzebowa¢ odwoªania tylko do warto±ci f, f

i f

Napisz funkcj¦ z metod¡ Halleya

function [ x , i t ]= h a l l e y (FCN,DFCN, D2FCN, x0 )

w m-pliku halley.m z parametrami:

•

FCN - 'wska¹nik' do funkcji (function handle) obliczaj¡cej warto±¢ f(x)

dla danego argumentu x,

•

DFCN - 'wska¹nik' do funkcji (function handle) obliczaj¡cej warto±¢

pochodnej f

(x)

dla danego argumentu x,

•

D2FCN - 'wska¹nik' do funkcji (function handle) obliczaj¡cej warto±¢

drugiej pochodnej f

(x)

dla danego argumentu x,

• x0

- liczba, przybli»enie startowe.

Funkcja powinna zwraca¢

• x

obliczone przybli»enie pierwiastka,

• it

- ilo±¢ iteracji.

Jako warunek stopu nale»y przyj¡¢ warunek speªnienie, którego± z wa-

runków |f(x

)| < 10

−7

lub |x

n+1

− x

| < 10

−10

Funkcja powinna si¦ zatrzyma¢, drukuj¡c komunikat o braku zbie»no±ci

na ekranie komputera, gdy ilo±¢ iteracji przekroczy 100.

W przypadku przekroczenia ilo±ci iteracji funkcja ma zwróci¢ komunikat

o tym na ekran.

Przetestuj metod¦ Halleya z wykorzystaniem funkcji halley() na przykªa-

dach równa« rozpatrywanych w zadaniach w 8.

W szczególno±ci nale»y sprawdzi¢, czy metoda zbiega kubicznie lokalnie,

o ile f

∗

) 6= 0

dla f(x) = x

− 4

dla x

∗

= 2

Tak, tego Halleya od komety Halleya.

Równania nieliniowe. Rozwikªywanie funkcji.

Rozpatrzmy dan¡ krzyw¡ okre±lon¡ równaniem f(x, y) = 0, gdzie funkcja
f (x, y)

jest okre±lona na prostok¡cie [a, b] × [c, d].

Chcemy znale¹¢ przybli»one warto±ci funkcji y(x) zadanej w sposób uwi-

kªany przez

f (x, y(x)) = 0

w punktach x

= a + k ∗ h

dla h = (b − a)/N, znaj¡c dobre przybli»enie

y(x

) = y(a) = y0

. Tutaj N - to ustalona liczba naturalna.

•

Napisz funkcj¦ (w m-pliku) octave'a, która dla danych parametrów

FCN wska¹nika do funkcji dwóch argumentów f(x, y),
ustalonych a, b - ko«ców odcinka
y0 b¦d¡cego przybli»eniem y

N ilo±ci punktów, w których chcemy znale¹¢ przybli»enie y

≈

y(x

)

- ten parametr mo»e by¢ opcjonalny. Je±li nie zostanie po-

dany to powinien przyjmowa¢ domy±ln¡ warto±¢ sto.

Funkcja powinna zwróci¢ jako wektor y

= y(x

)

dla k = 0, . . . , N.

W ka»dym kroku trzeba rozwi¡za¢, u»ywaj¡c funkcji octave'a fsolve ()),

równanie nieliniowe:

g(y

) = f (x

, y

) = 0

bior¡c za startowe przybli»enie y

k−1

(obliczone w poprzednim kroku dla

k − 1

Je±li si¦ oka»e, »e fsolve() nie potra rozwikªa¢ funkcji powinien na

ekranie pojawi¢ si¦ komunikat o bª¦dzie.

•

Testowa¢ na dwóch przykªadach:

Fragment elipsy:

f (x, y) = 2 ∗ x ∗ x + y ∗ y − 4.

Interesuje nas tu y(x) na [−1.5, 1.5] na siatce 101 punktowej. Za
y0

mo»emy wzi¡¢ y0 = 1

g(x, y) = x

+ y

− 4

na [0, 1] × [0, 1].

Metoda Householdera. LZNK

Zaªó»my, »e w wyniku do±wiadczenia otrzymujemy m punktów (x

, y

)

, które

powinny le»e¢ na jednej prostej, ale w wyniku bª¦dów pomiaru le»¡ tylko

blisko prostej. Chcemy wyznaczy¢ prost¡ y = ax + b, która le»y najbli»ej

tych punktów w sensie najmniejszych kwadratów, tzn. takie a, b, »e suma

k=1

|ax

+ b − y

jest minimalna. To zadanie mo»emy przedstawi¢ jako odpowiednie regularne

LZNK. Je±li istniej¡ dwa punkty o ró»nych odci¦tych, to LZNK ma jedno-

znaczne rozwi¡zanie.

Celem projektu jest zaprogramowanie funkcji rozwi¡zuj¡cej problem zna-

lezienia wspóªczynników prostej y = ax + b najlepiej przybli»aj¡cej dane

punkty (x

, y

) k = 1, . . . , m

za pomoc¡ rozkª¡du QR macierzy przy pomocy

metody Householdera. Tzn. szukamy (a, b) takich, »e

k=1

|ax

+ b − y

= min

a,ˆ

k=1

|ˆ

+ ˆ

b − y

Czyli: w funkcji rozwi¡zujemy LZNK A ∗ [a; b] ≈ ~y z macierz¡

A = [~

x, ~1]

dla wektorów

~1 =








1
1

...








x =








...








i wektora prawej strony

y =








...








Jako parametry wej±ciowe funkcji traktujemy wektory ~x, ~y dªugo±ci m,

funkcja powinna zwraca¢:

•

wektor [a; b] z rozwi¡zaniem tego LZNK,

•

macierz A z LZNK

•

macierz górnotrójk¡tn¡ R wymiaru 2 × 2 z rozkªadu QR macierzy A

metod¡ Householdera, tzn. A = Q ∗ [R; 0],

•

dwukolumnow¡ macierz B = [~h

, ~h

]

wymiaru m × 2 - w której od-

powiednie kolumny to wektory Householdera ~h

k = 1, 2

dla macierzy

Householdera H

takich, »e H

∗ H

= Q

Je±li kolumny macierzy A oka»¡ si¦ zale»ne liniowo - funkcja ma zwróci¢

odpowiedni komunikat na ekran.

Testy:

1. Przetestuj dla dwóch ró»nych punktów z ró»nymi x

- czy znajdzie

prost¡ przechodz¡c¡ przez te punkty.

2. Przetestuj dla punktów le»¡cych na prostej: tzn. wygeneruj kilka ró»-

nych losowych punktów x

i przyjmij y

= 2 ∗ x

− 10

. Nast¦pnie

sprawd¹, czy funkcja zwróci a = 2, b = −10.

3. Przetestuj dla punktów le»¡cych blisko danej prostej: tzn. przyjmij np.

= k/m

dla k = 1, . . . , m dla ró»nych m > 1 z y

= x

+ 2 +

dla

losowego z przedziaªu [−10

−3

, 10

−3

]

(funkcja octave'a rand() zwraca losowe punkty z przedziaªu [0, 1]).

4. Sprawd¹, czy rzeczywi±cie dla otrzymanych wektorów Householdera i

macierzy R, A zachodzi

A = H

∗ H

∗ [R; 0],

np. dla przykªadów z poprzednich podpunktów. W tym celu mo»na

stworzy¢ macierze H

i [R; 0] i te trzy macierze wymno»y¢.

Bibliograa

[1] Åke Björck i Germund Dahlquist. Metody numeryczne. Pa«-

stwowe Wydawnictwo Naukowe (PWN), Warszawa, 1983.

Przetªumaczone ze szwedzkiego przez Stefana Paszkowskiego.

[2] James W. Demmel. Applied numerical linear algebra. So-

ciety for Industrial and Applied Mathematics (SIAM), Phila-

delphia, PA, 1997.

[3] Maksymilian Dryja, Janina Jankowska, i Michaª Jankowski.

Metody numeryczne, tom 2. Wydawnictwo Naukowo Tech-

niczne (WNT), Warszawa, 1982.

[4] John W. Eaton. GNU Octave Manual. Network Theory Li-

mited, 2002.

[5] Gene Golub i James M. Ortega. Scientic computing. Aca-

demic Press Inc., Boston, MA, 1993. An introduction with

parallel computing.

[6] Gene H. Golub i James M. Ortega. Scientic computing

and dierential equations. Academic Press Inc., Boston, MA,

1992. An introduction to numerical methods.

[7] Wolfgang Hackbusch. Iterative solution of large sparse sys-

tems of equations, tom 95, Applied Mathematical Sciences.

Springer-Verlag, New York, 1994. Translated and revised

from the 1991 German original.

[8] Janina Jankowska i Michaª Jankowski. Metody numeryczne,

tom 1. Wydawnictwo Naukowo Techniczne (WNT), War-

szawa, 1981.

[9] C. T. Kelley. Solving nonlinear equations with Newton's me-

thod. Fundamentals of Algorithms. Society for Industrial and

Applied Mathematics (SIAM), Philadelphia, PA, 2003.

[10] A. Kieªbasi«ski i H. Schwetlick. Numeryczna algebra liniowa:

wprowadzenie do oblicze« zautomatyzowanych. Wydawnictwo

Naukowo Techniczne (WNT), Warszawa, 1992.

[11] David Kincaid i Ward Cheney. Analiza numeryczna. Wydaw-

nictwo Naukowo-Techniczne (WNT), 2006.

[12] J. M. Ortega i W. C. Rheinboldt. Iterative solution of non-

linear equations in several variables, tom 30, Classics in Ap-

plied Mathematics. Society for Industrial and Applied Ma-

thematics (SIAM), Philadelphia, PA, 2000. Reprint of the

1970 original.

[13] James M. Ortega. Numerical analysis, tom 3, Classics in

Applied Mathematics. Society for Industrial and Applied Ma-

thematics (SIAM), Philadelphia, PA, second edition, 1990. A

second course.

[14] Alo Quarteroni, Riccardo Sacco, i Fausto Saleri. Nume-

rical mathematics, tom 37, Texts in Applied Mathematics.

Springer-Verlag, New York, 2000.

[15] Alo Quarteroni i Fausto Saleri. Scientic computing with

MATLAB, tom 2, Texts in Computational Science and En-

gineering. Springer-Verlag, Berlin, 2003.

[16] Anthony Ralston. Wstep do analizy numerycznej. Pa«stwowe

Wydawnictwo Naukowe (PWN), Warszawa, 1983.

[17] J. Stoer i R. Bulirsch. Wstep do metod numerycznych. Tom

drugi. Pa«stwowe Wydawnictwo Naukowe (PWN), War-

szawa, 1980. Przetªumaczone z niemieckiego przez Maªgo-

rzate Mikulska.

[18] Lloyd N. Trefethen i David Bau, III. Numerical linear alge-

bra. Society for Industrial and Applied Mathematics (SIAM),

Philadelphia, PA, 1997.

Document Outline

Wstep
Wstepne zapoznanie sie z octavem
Obliczenia w arytmetyce zmiennopozycyjnej
Interpolacja Lagrange'a, bazy wielomianów
Splajny kubiczne i liniowe. Interpolacja splajnowa
Wielomiany Czebyszewa
Kwadratury
Rozwiazywanie równan nieliniowych
Układy równan liniowych - rozkłady typu LU i LL'
LZNK. Rozkład QR. Metoda Householdera
Numeryczne zadanie własne
Algorytm FFT
Przykładowe projekty zaliczeniowe

Marcinkowski L Zadania do wykładu z matematyki obliczeniowej

Document Outline