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
2 Wst¦pne zapoznanie si¦ z octavem 5
3 Obliczenia w arytmetyce zmiennopozycyjnej 9
4 Interpolacja Lagrange'a, bazy wielomianów 13
5 Splajny kubiczne i liniowe. Interpolacja splajnowa 20
6 Wielomiany Czebyszewa 25
7 Kwadratury 31
8 Rozwi¡zywanie równa« nieliniowych 37
9 Ukªady równa« liniowych - rozkªady typu LU i LL' 43
10 LZNK. Rozkªad QR. Metoda Householdera 48
11 Numeryczne zadanie wªasne 54
12 Algorytm FFT 60
13 Przykªadowe projekty zaliczeniowe 64
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(Dij)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))2na [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. x3+ x2− 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))2. 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))2dla argumentu b¦d¡cego macierz¡, tzn. je±li X = (xi,j)i,j macierz wymiaru M × N, to funkcja powinna zwróci¢ Y macierz wymiaru M × N tak¡, »e
Y = (yi,j)i,j, yi,j = f (xi,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))2.
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¦
SN = 1 + ... + N dla N = 100.
• U»yj instrukcji warunkowej i f ( )
endif
by sprawdzi¢, czy otrzymane SN 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 = 10kdla 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 + x2 w1 = x − a oraz
Algorytm 2
a =√
1 + x2 w2 = −1 x + a
dla x = 10k 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)4 = x4− ... + 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), f1(x) = a4;
Algorytm 2
f2(x) = x ∗ x ∗ x ∗ x − . . . + 16.
Matematycznie f1 ≡ f2, ale wyniki obliczone w arytmetyce zmien- nopozycyjnej mog¡ si¦ ró»ni¢.
Narysuj wykresy obu funkcji i policz bª¡d kf1(x) − f2(x)k∞, czyli maxk|f1(x(k) − f2(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)4 = x4− ... + 16 na 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=0xk/k!. Za odpo- wiedni¡ aproksymacje exp(x) bierzemy najpierw sto pierwszych elementów szeregu czyli przybli»amy exp(x) przez
FN(x) =
N
X
k=0
xk k!,
dla N = 100, a potem przybli»amy przez tysi¡c elementów sze- regu, czyli ustalamy N = 1000.
Sprawd¹ bª¡d wzgl¦dny |FN(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 In=R1
0 xn/(5+x)dx n = 0, .., 20dwoma algorytmami ze wzoru
In+ 5 ∗ In−1 = 1/n
Pierwszy algorytm przyjmuje I0 = log(6/5) i oblicza z powy»- szego wzoru kolejne
In = 1/n − 5 ∗ In−1
dla n = 1, 2, 3, . . ..
Drugi algorytm wykorzystuje fakt, »e 1
(n + 1)6 ≤ In ≤ 1
(n + 1)5 (3.1)
W tym algorytmie przyjmujemy za I30jak¡kolwiek warto±¢ z tego przedziaªu np. I30 = 1/180 i iterujemy w tyª, tzn. dla n = 30, 29, . . . , 20, . . . , 0 obliczamy
In−1= 1
5(1/n − In).
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)3 liczonej z wzoru na rozwini¦cie dwumianu x3− ... − 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)3 i x3− ... − 8na 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
(xj)nj=0
przestrzeni wielomianów stopnia nie wi¦kszego od n w odwrotnej kolejno±ci tzn.
(xn, xn−1, . . . , 1)
czyli np. wektor wspóªczynników (3, 2, 1) odpowiada wielomianowi 3x2+2x+
1. 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 x50− 1 dla x = −1, 0, 1.
Zadanie 2 Korzystaj¡c z funkcji polyval() narysuj wykres wielomianu x3+ 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 LnF dla funkcji F (x) = sin(x) dla w¦zªów −1, 0, 1. Policz warto-
±ci ró»nicy F − LnF w w¦zªach oraz korzystaj¡c z funkcji plot() narysuj wykresy F i LnF 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 LNf 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 − LNf na siatce tysi¡ca równoodlegªych punktów na tym odcinku, tzn.
eN = max |f (xk) − LNf (xk)|, gdzie xk to punkty siatki. Po- licz stosunek eN/e2N dla N < 64.
Czy bª¦dy malej¡ do zera? Jak zachowuje si¦ eN/e2N? 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:
Tn+1(t) = cos((n + 1)arccos(t)) na [−1, 1] odpowiednio przesuni¦te i przeskalowane.
Przetestuj, czy bª¦dy eN 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 LNf wielomianu interpolu- j¡cego funkcj¦ w tych w¦zªach oraz eN przybli»enie dyskretnej normy maksimum ró»nicy LNf − 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) =
n
X
k=0
akxk
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 + x2 oraz 1 − 2x + x2 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) =
n
X
k=0
akxk
oraz danej liczby q wspóªczynniki (bk)nk=0 wielomianu
p(x) =
n−1
X
k=0
bkxk 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) =
n
X
k=0
akxk
oraz liczby q warto±¢ w(q) i pochodnej w0(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) = x2, który w bazie Newtona zwi¡zanej z tymi w¦zªami ma nast¦puj¡c¡ posta¢: x2 = (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 = (ak)k taki, »e w(x) =
n
X
k=0
akxk
• y = (yk)k wektor wspóªczynników bazy Newtona.
Funkcja powinna zwróci¢ wektor wspóªczynników bk takich, »e w(x) =
n
X
k=0
bkwk, gdzie
wk(x) = Πk−1j=0(x − yj).
Zadanie 14 Napisz funkcj¦, która dla danego wielomianu w(x), którego wspóª- czynniki w bazie Newtona (wk)nk=0 znamy, oblicza wspóªczynniki tego wielomianu w bazie pot¦gowej (1, x, x2, . . . , xn). Tzn. pa- rametrami funkcji jest wektor wspóªczynników b = (bk)k takich,
»e
w(x) =
n
X
k=0
bkwk
i wektor y = (yk)k w¦zªów bazy Newtona (wk)nk=0 dla wk(x) = Πk−1j=0(x − yj).
Funkcja ma zwróci¢ wektor wspóªczynników a = (ak)k takich, »e
w(x) =
n
X
k=0
akxk.
Zadanie 15 Sprawd¹ eksperymentalnie ile wynosi dla ró»nych warto±ci N = 4, 8, 16, 32, 64, . . . przybli»enie:
Ar,N =
N
X
k=0
klkk∞,[−1,1]
dla {lk}Nk=0bazy Lagrange'a dla w¦zªów równoodlegªych na [−1, 1], tzn. dla xk = −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:
Ac,N =
N
X
k=0
klkk∞,[−1,1]
dla {lk}Nk=0 bazy Lagrange'a dla w¦zªów Czebyszewa na [−1, 1], tzn. dla xk zer wielomianu TN +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
kLNf k∞,[−5,5]
PN
k=0klkk∞,[−5,5]
dla f = 1/(1+x2)i LNf wielomianu interpoluj¡cego f w w¦zªach równoodlegªych na [−5, 5] dla N = 10, 20, 40, 80. Tutaj {lk}Nk=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 C2 na odcinku [a, b] i dla danego podziaªu tego odcinka:
a = x0 < x1. . . < xN = b
na pododcinki. Te funkcje obci¦te do ka»dego pododcinka [xi, xi+1]s¡ wielo- mianami kubicznymi.
Zadanie interpolacji splajnami kubicznymi polega na znalezieniu spaljnu kubicznego s speªniaj¡cego:
s(x0) = y0 s(x1) = y1
...
s(xN) = yN
dla zadanych warto±ci yk. 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 s1
na podziale równomiernym odcinka [−3, 3] z w¦zªami {xk = k}
dla k = −3, −2, . . . , 3 przyjmuj¡cego warto±ci s1(xk) = (−1)k w tych w¦zªach.
Nast¦pnie znajd¹ wspóªczynniki splajnu kubicznego s2 na tym samym podziale odcinka i przyjmuj¡cego te same warto±ci w w¦- zªach co s1, 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 s1 i s2 na tym samym ry- sunku.
Czy otrzymali±my te same splajny?
Policz przybli»on¡ norm¦ maksimum ró»nicy s1− s2 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 SN na N w¦- zªach równoodlegªych dla funkcji f(x) = sin(x) na odcinku [−π, 2∗
π] dla N = 2kN0 dla N0 = 5 i k = 1, 2, 3, 4, 5.
Nast¦pnie
• narysuj wykresy funkcji f(x) i splajnów SN dla ró»nych N.
• oblicz dyskretn¡ norm¦ maksimum na siatce równomiernej zªo»onej z tysi¡ca punktów na tym odcinku, tzn. eN = maxk| sin(xk) − SN(xk)| dla xk punktów siatki.
• policz równocze±nie wspóªczynnik ee2NN . Czy prawd¡ jest, »e eN
e2N ≈ 2p
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 ee2NN s¡ te same? Tzn. czy szybko±¢ zbie»no±ci k sin(x) − SNk∞ 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 SN na N w¦zªach równoodlegªych dla f na odcinku [−5, 5] dla N = 2kN0 dla N0 = 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. eN = max | sin(xk) − SNsin(xk)| dla xk = −5 + k ∗ 0.01z k = 0, . . . , 1000.
Policz równocze±nie wspóªczynnik ee2NN . Czy ee2NN ≈ 2pdla jakiego±
p 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) = x2 lub (x + 1)3.
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)3. Nast¦pnie sprawd¹ wspóªczynniki s w bazie {(x − xk)j}j=0,1,2,3 za pomoc¡ umkpp() na obu przedziaªach, tzn. dla x0 = −1 na przedziale [−1, 0] i x1 = 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}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 sn 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 eN = k sin −snk∞.
• Policz wspóªczynnik ee2NN . Czy wida¢, »e eN
e2N ≈ 2p
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
• xwektora N + 1 ró»nych w¦zªów uszeregowanych (a = x0 <
x1. . . < xN = 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 xk.
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
fj(x) = (x)j+ = 0 x < 0
xj 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
{xk = −π + 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 sNfj.
• policz przybli»on¡ norm¦ maksimum bª¦du (na 1000 równo- miernych punktach), tzn. przybli»enie
eN = ksNfj − fjk∞,[a,b].
• policz wspóªczynniki ee2NN . Czy wida¢, »e eN
e2N ≈ 2p
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, T0 ≡ 1 oraz Tn+1(x) = 2xTn(x) − Tn−1(x) (6.1) albo ze wzoru:
Tn(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¡:
• xargument dla którego chcemy obliczy¢ warto±¢ wielomianu
• a - wektor zawieraj¡cy wspóªczynniki a0, . . . , an takie, »e w(x) =Pn
k=0akTk(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) = T15(x), tzn. narysuj dwa wy- kresy T15(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 = (xij)ij, wymiaru k × l
• wektor wspóªczynników
a = [a0, . . . , an] takich, »e
w(x) =
n
X
k=0
akTk(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 = (yij)ij, wymiaru k × l z warto±ciami yij = w(xij).
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 T3(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 T32 ze wzoru
T32(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. T7, T50. Zadanie 4 Narysuj korzystaj¡c z wzoru
Tn(x) = cos(n ∗ arccos(x))
wykres wielomianów Tn 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 T10 i T15 ze wzoru Tn(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−nTN +1(x) ma minimaln¡ norm¦ supremum na [−1, 1] w±ród wszystkich wie- lomianów postaci ΠNk=0(x − xk).
Policz dla N + 1 ró»nych losowych w¦zªów xk, 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 Πnk=0(x−xk). Przeprowad¹ testy równie» dla w¦zªów Czebyszewa (czyli zer TN +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−nTN +1(x)ma minimaln¡ norm¦
supremum na [−1, 1] w±ród wszystkich wielomianów postaci
w(x) = xn+1+
n
X
k=0
akxk.
Chcemy sprawdzi¢ eksperymentalnie, czy ta wªasno±¢ si¦ potwier- dzi.
Policz dla N + 1 losowych wspóªczynników ak dla k = 0, . . . , N przybli»on¡ norm¦ supremum wielomianu:
w(x) = xn+1+
n
X
k=0
akxk.
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
w∗(x) = Tn(x)/Tn(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 ak 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 +
n
X
k=1
ak(x − 2)k
i porównaj z norm¡ supremum w∗ na [−1, 1] równ¡ 1/Tn(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 L2w(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 L2w(a, b) z wag¡ w(x) tzn.:
(f, g)L2w(a,b)= Z b
a
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 L2√1
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 L2√1
1−x2
(−1, 1), tzn. czy prawd¡ jest, »e
Z 1
−1
Tn(x)Tm(x) 1
√1 − x2 dx =
0 m 6= n,
π
2 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 Z b
a
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¦ Rabf (t) dt dla nast¦puj¡cych funkcji i odcinków:
• x2 na [−1, 1]
• x9 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:
Tnf = h 0.5[f (a) + f (b)] +
n−1
X
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 Nk= 2kN0 k = 0, 1, 2, . . .z ustalonym N0 = 5 licz¡c:
EN E2N dla bª¦du
EN = | Z b
a
f dt − TNf |
dla nast¦puj¡cych funkcji, dla których warto±ci caªek znamy:
• x2 na [0, 1],
• sin(x) na [0, π] i N0 = 5 - funkcja analityczna,
• sin(100 ∗ x) na [0, π] i N0 = 5 - funkcja analityczna, silnie oscyluj¡ca - du»e warto±ci drugiej pochodnej,
• xj+0.5 na [0, 1] dla j = 0, 1, 2 czyli funkcji w Cj([0, 1]) i 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 L2 na [−1, 1] z wag¡ √1−x1 2.
Policz za pomoc¡ funkcji octave quad():
Z 1
−1
√ 1
1 − x ∗ xTn(x)Tm(x) dx dla m = 2 i n = 3. Czy wynik jest zgodny z teori¡?
Zadanie 4 Kwadratura Gaussa-Czebyszewa dla
I(f ) = Z 1
−1
√ 1
1 − x2f (x) dx czyli
GCn+1f = π n + 1
n
X
k=0
f (xk) = π n + 1
n
X
k=0
f
cos 0.5 ∗ π + kπ n + 1
dla xk zer Tn+1 - n + 1 wielomianu Czebyszewa.
Zaimplementuj funkcj¦ function c=gaussczeb(FCN,n) oblicza- j¡c¡ warto±¢ kwadratury GCnf. Parametry - wska¹nik do funkcji function y=f(x) i ilo±¢ punktów kwadratury.
Przetestuj jej dziaªanie:
• policz przybli»enia caªek I(Tk2). Czy GCn(Tk2)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
en= |GCn(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. GCn(f ), dla funkcji f(x) = exp(arccos(x)) z wynikiem obliczonym funkcj¡ quad(). Wypisz na ekran en= |GCn(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 en, analogicznie do poprzednich podpunktów i wypisz en/e2n.
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:
Pnf = h
n
X
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 Nk = 2kN0 k = 0, 1, 2, . . .z ustalonym N0 = 5 licz¡c:
EN
E2N, EN = | Z b
a
f dt − PNf |
dla caªek z nast¦puj¡cych funkcji po odpowiednich odcin- kach:
• x2 na [0, 1]
• sin(x) na [0, π] i N0 = 5 - funkcja analityczna,
• sin(100∗x)na [0, π] i N0 = 5- funkcja analityczna silnie oscyluj¡ca (du»e warto±ci drugiej pochodnej),
• xj+0.5 na [0, 1] dla j = 0, 1, 2 czyli funkcji w Cj([0, 1]) i 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 PNf z TN −1f,
• porównaj obie kwadratury dla tego samego h, tzn. TNf z PNf, 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 Rp,N 0obliczaj¡c¡ przybli»enie caªki Rb
a 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
• pilo±¢ 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 Rabf 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, π]
• xj+0.5 na [0, 1] dla j = 0, 1, 2, . . . , 10, dla N0 = 100 i dla p = 2, 3, 4, 5, 6.
Czy dla √
x i rosn¡cego p bª¡d maleje?
Czy stosunek EN/E2N maleje tak samo dla wszystkich funk- cji?
Tutaj u»yli±my oznaczenia na bª¡d EN = |Rb
a f dt − Rp,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:
x2 − 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,
• x2− 2 = 0,
• x10− 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 x0 = [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 x0 = 2, f2(x)= x ∗ x ∗ x − 27 z x0 = 27,
f3(x)= exp(x) − 2 z x0 = 10, −10, 1e3, f4(x)= sin(x) dla x0 = 2,
f5(x)= cos(x) z x0 = 2
f6(x)= (x − 2)k x x0 = 3 dla k = 2, 4, 8, 16,
f7(x)= x ∗ x − 2 z x0 = 106, 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¢ x0?
Dla wszystkich tych funkcji znamy rozwi¡zania, wi¦c mo»na wy-
±wietla¢ równocze±nie na ekranie bª¡d en = xn− 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:
|en|/|en−1|p dla p = 1, 2, 3.
Prosz¦ powtórzy¢ testy z innymi ró»nymi warto±ciami startowymi x0.
Zadanie 7 Powtórz zadanie 6 zast¦puj¡c metod¦ Newtona przez metod¦
siecznych.
Przetestuj, czy
en/(en−1en−2) asymptotycznie zbiega do staªej, i czy
|en|/|en−1|p dla p = (1 +√
5)/2 zbiega do staªej np. dla x ∗ x − 2 z x0 = 2 lub innych równa« z poprzedniego zadania.
Za x1 mo»emy przyj¡¢ x0+ 10−7.
Zadanie 8 Sprawd¹, czy metoda iteracji prostych xn= cos(xn−1) zbiega do x∗ = cos(x∗).
Zbadaj eksperymentalnie, czy zbie»no±¢ jest liniowa, tzn. czy
|x∗− xn|
|x∗− xn−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.
xn+1= xn− f (xn)
f (xn+h)−f (xn) h
= xn− f (xn) ∗ h f (xn+ h) − f (xn) 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) = 2x2+ 3y2− 3 = 0
znajd¹ przybli»one warto±ci yk = y(xk) dla xk = k ∗ h dla k =
−N, . . . , N i h = 1/N, tzn. speªniaj¡ce g(xk, yk) = 0. Testuj dla ró»nych warto±ci N, np. N = 10, 20, 40, 80, ....
oblicz yk rozwi¡zuj¡c równanie
g(xk, yk) = 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 f1(x1, x2) = x1+ 2x2 = 1 f2(x1, x2) = 3 ∗ x21+ x22 = 1 dla ró»nych przybli»e« pocz¡tkowych.
W przypadku zbie»no±ci policz kx − xnkp kx − xn−1kqp
dla p = 1, 2, ∞ oraz q = 1, 2, 3. Czy ten iloraz zbiega do staªej dla jakiego± z tych q?
Funkcj¦ F (x) := (f1(x), f2(x))T 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
xn= xn− aF (x)
do ukªadu z zadania 12, tzn. przyjmijmy, »e F (x) := (f1(x), f2(x))T, 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 − xnkp
kx − xn−1kp 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 LLT 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 = AT > 0, czyli A = RTR, (zazwyczaj w literaturze podaje si¦ ten rozkªad w terminach czynnika L = RT),
• 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 kAk2∗ kA−1k2.
Zadanie 1 Operator octave'a \ sªu»¡cy m.in. do rozwi¡zywaniu ukªadów równa« liniowych w octavie.