• Nie Znaleziono Wyników

Zadania i scenariusze zaj¦¢ z laboratorium komputerowego do wykªadu z Matematyki Obliczeniowej

N/A
N/A
Protected

Academic year: 2022

Share "Zadania i scenariusze zaj¦¢ z laboratorium komputerowego do wykªadu z Matematyki Obliczeniowej"

Copied!
76
0
0

Pełen tekst

(1)

Zadania i scenariusze zaj¦¢ z laboratorium komputerowego do wykªadu z Matematyki

Obliczeniowej

Leszek Marcinkowski

12 grudnia 2011

(2)

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.

(3)

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

(4)

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

(5)

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.

(6)

Scenariusze i zadania

komputerowe

(7)

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.

(8)

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))

(9)

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

(10)

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 ( )

(11)

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)?

(12)

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;

(13)

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 Pk=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

(14)

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?

(15)

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.

(16)

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):

(17)

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?

(18)

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.

(19)

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

(20)

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]

(21)

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.

(22)

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}

(23)

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

(24)

• 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.

(25)

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?

(26)

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?

(27)

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)

(28)

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))

(29)

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].

(30)

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:

(31)

• 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 L21

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 L21

1−x2

(−1, 1), tzn. czy prawd¡ jest, »e

(32)

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.

(33)

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)

(34)

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

(35)

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¡.

(36)

• 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.

(37)

(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

(38)

• 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|.

(39)

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 ) ,

(40)

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.

(41)

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.

(42)

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?

(43)

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)

(44)

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.

(45)

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.

Cytaty

Outline

Powiązane dokumenty

Proszę napisać funkcje, która mnoży dwie liczby typu unsigned (albo unsigned long), a wynik zapamiętuje w dwóch innych zmiennych, z których jedna przechowuje najmłodsze bity

W przypadku funkcji generującej liczby losowe metodą Boxa-Mullera, można zadbać o to, aby przy nieparzystym wywołaniu była wykonywana cała pro- cedura począwszy od losowania liczb x

Napisać program, który czyta ciąg liczb ze standardowego wejścia aż do wystąpie- nia znaku końca pliku (kombinacja klawiszy Ctrl–D na klawiaturze), a następnie oblicza

B¦dziemy starali si¦ zbada¢ jak mo»emy uzyska¢ jeden w¦zeª z dwóch i jak mo»na rozªo»y¢ w¦zeª na dwa prostsze w¦zªy. B¦dzie nam wygodniej patrze¢ na w¦zeª jako le»¡cy

Przegląd funkcji elementarnych – zadanie domowe do samodzielnego (!!) przygotowania- powtórzenie wiadomości ze szkoły średniej... Określić dziedzinę funkcji f: R→R,

Montrealu jest właśnie południe, u mnie w Paryżu jest godzina 18 po południu tego samego dnia?. Jaka była średnia prędkość tego samochodu na

Powtórz zadanie dla innych zadanych przez prowadz¡cego lub samodzielnie wybranych

i informatyki na potrzeby gospodarki - Wiking Projekt jest wspóªnansowany z Europejskiego Funduszu Spoªecznego w ramach programu operacyjnego KAPITAŠ LUDZKI Poddziaªanie