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¦ R_a^bf (t) dt dla nast¦puj¡cych funkcji i odcinków:

• x² na [−1, 1]

• x⁹ na [0, 1]

• sin(x) na [0, π]

• cos(100 ∗ x) na [0, π]

• cos(100 ∗ x) ∗ cos(1000 ∗ x) na [0, π]

Czy wyniki s¡ zgodne z teori¡?

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

(a) Zaprogramuj w octavie funkcj¦

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

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

T_nf = 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 N_k= 2^kN₀ k = 0, 1, 2, . . .z ustalonym N0 = 5 licz¡c:

E_N E_2N dla bª¦du

E_N = | Z b

a

f dt − T_Nf |

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

• x² 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,

• x^j+0.5 na [0, 1] dla j = 0, 1, 2 czyli funkcji w C^j([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 L² na [−1, 1] z wag¡ ^√_1−x¹ 2.

Policz za pomoc¡ funkcji octave quad():

Z 1

−1

√ 1

1 − x ∗ xT_n(x)T_m(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 − x²f (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(Tk²). Czy GCn(T_k²)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

e_n= |GC_n(f ) − I(f )|

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

• analogicznie do poprzedniego podpunktu tego zadania, poró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 e_n= |GC_n(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/e_2n.

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 = 2^kN₀ k = 0, 1, 2, . . .z ustalonym N0 = 5 licz¡c:

EN

E_2N, EN = | Z b

a

f dt − PNf |

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

• x² 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),

• x^j+0.5 na [0, 1] dla j = 0, 1, 2 czyli funkcji w C^j([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 P_Nf, 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 R_a^bf dt dla na-st¦puj¡cych funkcji, odcinków i warto±ci parametrów N0 i p:

• sin(x) dla [0, π]

• sin(x) dla [0, 100 ∗ π]

• sin(10 ∗ x) dla [0, π]

• x^j+0.5 na [0, 1] dla j = 0, 1, 2, . . . , 10, dla N0 = 100 i dla p = 2, 3, 4, 5, 6.

Czy dla √

x i rosn¡cego p bª¡d maleje?

Czy stosunek EN/E_2N maleje tak samo dla wszystkich funk-cji?

Tutaj u»yli±my oznaczenia na bª¡d EN = |Rb

a f dt − R_{p,N 0}|.

Rozdziaª 8

Rozwi¡zywanie równa«

nieliniowych

W poni»szych zadaniach przetestujemy jak dziaªaj¡ cztery metody rozwi¡-zywania równa« nieliniowych skalarnych, tzn. metoda bisekcji, metoda New-tona, metoda siecznych i metoda iteracji prostych oraz dwie metody rozwi¡-zywania ukªadów równa« nieliniowych, czyli wielowymiarowa metoda New-tona i metoda iteracji prostych w najprostszej formie.

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

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

nieliniowych.

Zadanie 1 Zaimplementuj metod¦ bisekcji w skrypcie - dla rozwi¡zania rów-nania cos(x) = 0 na odcinku [0, 2]. Przetestuj, czy funkcja znaj-dzie dobre przybli»enie π/2.

Zadanie 2 Przetestuj metod¦ bisekcji z poprzedniego zadania do rozwi¡zania innych problemów:

x² − 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⁻⁵.

Zadanie 4 Zapoznaj si¦ z pomoc¡ dla funkcji octave'a fzero() - rozwi¡zu-j¡c¡ skalarne równania nieliniowe, przetestuj j¡ na kilku prostych przykªadach skalarnych np.

• cos(x) = 0,

• x²− 2 = 0,

• x¹⁰− 2 = 0,

• x − sin(x) = 0

i innych prostych równaniach nieliniowych.

Prosz¦ przetestowa¢ t¦ funkcj¦ dla ró»nych warto±ci pocz¡tko-wego przedziaªu 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 = 10⁶, sprawd¹, czy metoda Newtona zbiega, a je±li tak - to czy zbiega ona kwadratowo, f8(x)= 1/x − a dla danych a = 0.5, 2, 4, 100 (oczywi±cie

imple-mentuj¡c bez dzielenia) jak dobra¢ x0?

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

wy-±wietla¢ równocze±nie na ekranie bª¡d e_n = x_n− r

(tutaj r - znane rozwi¡zanie) i bada¢ eksperymentalnie rz¡d zbie»-no±ci, tzn. drukowa¢ na ekranie stosunki bª¦du bie»¡cego do po-przedniego w odpowiedniej pot¦dze:

|e_n|/|e_n−1|^p dla p = 1, 2, 3.

Prosz¦ powtórzy¢ testy z innymi ró»nymi warto±ciami startowymi x₀.

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

siecznych.

Przetestuj, czy

e_n/(e_n−1e_n−2) asymptotycznie zbiega do staªej, i czy

|e_n|/|e_n−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⁻⁷.

Zadanie 8 Sprawd¹, czy metoda iteracji prostych x_n= cos(x_n−1) zbiega do x^∗ = cos(x^∗).

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

|x^∗− x_n|

|x^∗− x_n−1|

zbiega do staªej. Za dobre przybli»enie x^∗ mo»na przyj¡¢ rozwi¡-zanie równania obliczone za pomoc¡ wywoªania funkcji octave'a:

fzero().

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

przybli»amy ilorazem ró»nicowym, tzn.

x_n+1= x_n− f (x_n)

f (xn+h)−f (xn) h

= x_n− f (x_n) ∗ h f (x_n+ h) − f (x_n) dla ustalonego h.

Przetestuj ró»ne h, np. h = 10⁻⁴, 10⁻⁷, 10⁻¹⁰ itp.

Porówna¢ zbie»no±¢ z dokªadn¡ metod¡ Newtona (szczególnie ostatnie iteracje) dla funkcji z zadania 6.

Zadanie 10 Rozwikªywanie funkcji:

Dla funkcji y(x) zadanej w sposób uwikªany równaniem g(x, y) = 2x²+ 3y²− 3 = 0

znajd¹ przybli»one warto±ci yk = y(x_k) dla xk = k ∗ h dla k =

−N, . . . , N i h = 1/N, tzn. speªniaj¡ce g(xk, y_k) = 0. Testuj dla ró»nych warto±ci N, np. N = 10, 20, 40, 80, ....

oblicz yk rozwi¡zuj¡c równanie

g(x_k, y_k) = 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⁻¹ na odcinku [0, 5] na siatce k ∗ h dla k = 0, .., 100. Narysuj wykres funkcji f⁻¹.

Sam wykres mo»na narysowa¢ du»o pro±ciej bez wyliczania war-to±ci f⁻¹. 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 f₁(x₁, x₂) = x₁+ 2x₂ = 1 f₂(x₁, x₂) = 3 ∗ x²₁+ x²₂ = 1 dla ró»nych przybli»e« pocz¡tkowych.

W przypadku zbie»no±ci policz kx − x_nk_p kx − x_n−1k^qp

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

x_n= x_n− aF (x)

do ukªadu z zadania 12, tzn. przyjmijmy, »e F (x) := (f1(x), f₂(x))^T, a parametr a 6= 0 np. a = 1, −1, 10⁻¹ itp.

Czy metoda zawsze zbiega, a je±li tak - to jak szybko? W przy-padku zbie»no±ci policz

kx − xnkp

kx − x_n−1k_p dla p = 1, 2, ∞.

Zadanie 14 Do ukªadu równa« z zadania 12 zastosuj funkcj¦ fsolve () - po-równaj z implementacj¡ metody Newtona w zadaniu 12, która byªa zastosowana do rozwi¡zania tego ukªadu.

Czy obie metody zbiegaj¡ do tego samego rozwi¡zania dla tych samych warto±ci wektorów startowych?

Sprawd¹ - za pomoc¡ funkcji tic i toc - czas potrzebny do roz-wi¡zania tego ukªadu równa« za pomoc¡ obu metod.

Rozdziaª 9

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

W tym rozdziale zapoznamy si¦ z metodami sªu»¡cych do rozwi¡zywania ukªadów równa« liniowych przy pomocy uzyskiwaniu odpowiednich rozkªa-dów macierzy typu LU i LL^T oraz obliczaniu macierzy odwrotnej.

Zapoznamy si¦ z odpowiednimi funkcjami octave'a sªu»¡cymi znajdowa-niu odpowiednich rozkªadów, czy rozwi¡zywaznajdowa-niu ukªadów równa« liniowych z pomoc¡ tych rozkªadów.

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

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

• [R]=chol(A) - funkcja zwracaj¡ca czynnik rozkªadu Choleskiego ma-cierzy symetrycznej dodatnio okre±lonej A = A^T > 0, czyli A = R^TR, (zazwyczaj w literaturze podaje si¦ ten rozkªad w terminach czynnika L = R^T),

• 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⁻¹k₂.

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

Przetestuj operator octave'a \ dla ukªadu równa« z macierz¡

A = [1, 2; 3, 4]i x = [1; 3] z f=A∗x. Czy y=A\f jest rozwi¡zaniem tego ukªadu?

Policz normy residualn¡ kAy − fkp i norm¦ bª¦du kx − xkp dla p = 1, 2.

Powtórz testy dla macierzy osobliwej [1, 2; 3, 6] i prawie osobliwej [1, 2; 3, 6 − 10⁻⁶]. Co si¦ wtedy dzieje?

Zadanie 2 Przetestuj solver octave'a, tzn operator backslash dla dwóch prostych ukªadów równa« liniowych z macierzami nieosobliwymi:

A1 = [1, 1; 1, 0.98]i b = [2; 1], A2 = [1, 1; 1, 0.99] i wektorem pra-wej strony ukªadu równa« liniowych b = [2; 1] jak w poprzed-nim zadaniu. Policz w normie drugiej ró»nic¦ rozwi¡za«, norm¦

Frobeniusa ró»nicy A1 − A2 oraz uwarunkowania tych macierzy korzystaj¡c z funkcji octave'a cond().

Powtórz zadanie dla macierzy: A1 = [1, 1; 1, 1 − 2] i b = [2; 1], A₂ = [1, 1; 1, 1 − ] dla  = 10^−p z p = 3, 4, . . . , 16.

Zadanie 3 Funkcja inv().

Zapoznaj si¦ z pomoc¡ do funkcji: inv(). Przetestuj, obliczaj¡c macierz odwrotn¡ do A = [1, 1; 1, −1]. Czy macierz B obliczona za pomoc¡ tej funkcji rzeczywi±cie jest macierz¡ odwrotn¡?

Policz normy pierwsz¡ i Frobeniusa kA∗B −Ik oraz kB ∗A−Ik.

Zastosuj funkcje do rozwi¡zywania ukªadu równa« liniowych Ax = f dla znanego x (liczymy f = A ∗ x). Policz y jako iloczyn B z f i policz bª¦dy residualny kAy − fk i kx − yk w normie pierwszej i drugiej.

Powtórz to zadanie dla macierzy N ×N losowej (funkcja octave'a rand() zwraca macierz losow¡) dla n = 10, 50, 250, 1250 ze zna-nym rozwi¡zaniem - oszacuj czas przy pomocy funkcji tic i toc.

Porównaj czas i bª¦dy w normie pierwszej dla rozwi¡zania uzy-skanego przy pomocy operatora octave'a \.

Zadanie 4 Funkcje lu() i chol().

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

• Dla macierzy A = [1, 1; 1, −1] znajd¹ jej rozkªady: P A = LU i rozkªad Choleskiego A = L^T1L₁.

• Sprawd¹ te rozkªady licz¡c normy macierzowe pierwsz¡ i Frobeniusa bª¦dów np. P A − LU i A − L^T1L₁.

• Zastosuj te rozkªady do znalezienia rozwi¡zania ukªadu rów-na« liniowych Ax = f ze znanym rozwi¡zaniem np. x = [1; 1] i f = [2; 0].

• Policz normy pierwsz¡ i drug¡ wektorowe pomi¦dzy x, a wynikiem algorytmu w, polegaj¡cym na zastosowaniu od-powiedniego rozkªadu, oraz takie same normy residualne, tzn. normy Aw − f.

Zadanie 5 W octavie przetestuj eliminacje Gaussa z cz¦±ciowym wyborem i bez wyboru dla macierzy A = [e, 1; 1, 1] z e = eps/10 (funk-cja octave'a eps zwraca epsilon maszynowy) i wektorem prawej strony f = [1; 0]^T. Trzeba tu zaprogramowa¢ samodzielnie elimi-nacj¦ Gaussa bez wyboru elementu gªównego, tzn. bez permuta-cji ani wierszy, ani kolumn dla macierzy 2 × 2.

Zadanie 6 Testy solvera octave'a dla macierzy Hilberta H(N) Polecenie octave'a hilb() tworzy macierz Hilberta.

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

• Dla znanego rozwi¡zania, np. staªego równego jeden na ka»-dej pozycji, czyli solk= 1, policz H(N) ∗ sol = f,

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

• Policz normy (1,2 itd) kH(N)x − fk i kx − solk dla x roz-wi¡zania wyliczonego przez octave,

• Policz uwarunkowanie macierzy Hilberta dla powy»szych N.

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

Nale»y tu sztucznie wymusi¢, za pomoc¡ funkcji octave'a single (), aby zmienne byªy zmiennymi pojedynczej precyzji.

Zadanie 8 Przetestuj funkcj¦ invhilb() tworz¡c¡ macierz odwrotn¡ do ma-cierzy Hilberta (por. zadanie 6).

• Policz normy: macierzowa norm¦ Frobeniusa i norm¦ macie-rzow¡ pierwsz¡ ró»nicy macierzy E = H(N)∗iH(N)−I dla N = 10, 20, 30, gdzie iH(N) macierz odwrotna do macierzy Hilberta obliczona przez invhilb().

• Wykorzystaj iH(N) do rozwi¡zania ukªadu równa« z ma-cierz¡ Hilberta H(N), tzn.:

(a) stwórz macierze H(N) oraz iH(N) dla N = 10, 20, 30 korzystaj¡c z funkcji hilb() i invhilb(),

(b) dla znanego rozwi¡zania sol, policz H(N) ∗ sol = f, (c) rozwi¡» ukªad H(N)x = f mno»¡c f przez iH(N), tzn.

x = iH(N ) ∗ f,

(d) policz normy kH(N)x−fkp i kx−solkpdla p = 1, 2, ∞.

Zadanie 9 Stwórz w octavie macierz trójdiagonaln¡ A = (aij)ⁿ_i,j=1 wymiaru n × n z 2 na diagonali i −1 na pod- i nad diagonal¡, tzn.

ai,j =







2 i = j

−1 |i − j| = 1 0 |i − j| > 1

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

przy pomocy funkcji octave'a diag(), jak i wprost w p¦tli. Po-równaj czas u»ywaj¡c tic i toc dla n = 10, 100, 1000.

Policz uwarunkowanie macierzy dla ró»nych N.

Policz macierz odwrotn¡ przy pomocy inv() (czy jest trójdiago-nalna?).

Zadanie 10 Sprawd¹ z wykorzystaniem funkcji octave'a chol(), czy macierz z zadania 9 jest symetryczna i dodatnio okre±lona.

Zadanie 11 Zaimplementuj rozkªad LU dla macierzy trójdiagonalnej N × N bez wyboru elementu gªównego, tzn. stwórz wªasn¡ funkcj¦

octave'a:

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

• a, b, c- przek¡tna, podprzek¡tna i nadprzek¡tna macierzy A,

• f - wektor prawej strony,

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

Funkcja powinna zwróci¢ x - rozwi¡zanie Ax = f oraz czynniki rozkªadu macierzy A = LU, czyli diagonal¦ macierzy górnotrój-k¡tnej U w wektorze d i poddiagonal¦ L - macierzy dolnotrójk¡t-nej, trójdiagonalnej z jedynkami na diagonali, czyli wektor l.

Naddiagonala U równa si¦ nadiagonali A, tzn. wektorowi c.

Zadanie 12 Zastosuj funkcj¦ z zadania 11 do uzyskania rozkªadu LU macierzy z zadania 9 dla N = 4, 10, 15. Porównaj z wynikiem uzyskanym przy pomocy funkcji octave'a lu().

Zadanie 13 Zastosuj funkcj¦ z zadania 11 do rozwi¡zania ukªadu równa« z macierz¡ z zadania 9 dla N = 10, 100, 1000, czy 10000:

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

• Porównaj czas rozwi¡zywania wªasnym algorytmem i algo-rytmem octave'a, czyli operatorem A\f,

• Policz bª¡d rezydualny i bª¡d rzeczywisty w normie drugiej, tzn. kx − ˜xk2 i kf − A ∗ ˜xk2, gdzie f = A ∗ x, a ˜x to rozwi¡zanie obliczone przez octave'a lub obliczone wªasnym algorytmem.

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

cz¦-±ciowym wyborem elementu gªównego.

Nast¦pnie testuj ten rozkªad jak w zadaniu 13.

Zastosuj ten rozkªad do rozwi¡zania ukªadu równa« z macierz¡ z poprzedniego zadania A i z macierz¡ A − 1.5 ∗ I.

Zadanie 15 Zaprogramuj rozkªad Choleskiego typu dla macierzy A = A^T > 0 trójdiagonalnej, tzn. policz L dolnotrójk¡tn¡ z jedn¡ poddia-gonal¡ (czyli reprezentowan¡ przez dwa wektory z przek¡tn¡ i podprzek¡tn¡) tak¡, »e A = L^TL.

Zastosuj do rozwi¡zania ukªadu równa« z zadania 9.

Testuj analogicznie do zadania 13.

Porównaj macierz L uzyskan¡ w ten sposób z macierz¡ uzyskan¡

przez chol() zastosowan¡ do tej samej macierzy.

Rozdziaª 10

LZNK. Rozkªad QR. Metoda Householdera

W tym rozdziale zajmiemy si¦ liniowym zadaniem najmniejszych kwadratów (LZNK).

Dla danej macierzy A wymiaru M × N i wektora b wymiaru M chcemy znale¹¢ wektor x wymiaru N taki, »e

kAx − bk₂ = min

y kAy − bk₂.

Je±li A jest macierz¡ kolumnami regularn¡ (rz¡d A jest maksymalny równy N), to zadanie to ma jednoznaczne rozwi¡zanie i nazywamy je regularnym LZNK (RLZNK).

Podstawowym algorytmem sªu»¡cym jego rozwi¡zaniu jest metoda Ho-useholdera, czyli znalezienia rozkªadu QR macierzy A, gdzie Q to macierz ortogonalna - iloczyn macierzy Householdera, a R to macierz górnotrójk¡tna.

W tym rozdziale przetestujemy podstawowy operator octave'a sªu»acy do rozwi¡zywania dowolnego LZNK, tzn. operator \. Zauwa»my, »e je±li M = N i A jest macierz¡ kolumnami regularn¡, to A jest nieosobliwa i to regularne LZNK jest równowa»ne rozwi¡zaniu ukªadu równa« liniowych Ax = b.

Przetestujemy w zadaniach funkcj¦ octave'a qr() sªu»¡c¡ znajdowaniu rozkªadu QR macierzy, kilka wªasno±ci macierzy (przeksztaªce«) Househol-dera i rozwi¡»emy kilka konkretnych LZNK.

Zadanie 1 Operator octave'a \ sªu»¡cy m.in. do rozwi¡zywaniu ukªadów równa« liniowych i LZNK w octavie.

• Przetestuj operator octave'a \ rozwi¡zuj¡c RLZNK dla ma-cierzy A ze znanym konkretnym rozwi¡zaniem x:

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

przyjmuj¡c, »e f = Ax.

Czy y=A\f jest rozwi¡zaniem tego ukªadu?

Policz normy residualn¡ kAy − fk2 i norm¦ bª¦du kx − yk2.

• Przetestuj ten operator dla nieregularnego LZNK dla ma-cierzy A^T z x = [1; 1; 1] i f = Ax.

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

Policz normy residualn¡ kA^Ty − f k₂ i norm¦ bª¦du kx−yk2. Zadanie 2 Rozkªad QR w octave. Funkcje qr().

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

(b) Dla macierzy A = [1, 1; 1, −1; 1, 3] znajd¹ jej rozkªad A = QR z pomoc¡ funkcji qr().

(c) Sprawd¹, czy uzyskana Q jest ortogonalna - policz normy Frobeniusa QQ^T − I i Q^TQ − I.

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

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

(e) Zastosuj ten rozkªad do znalezienia rozwi¡zania LZNK Ax = f ze znanym rozwi¡zaniem, np. x = [1; 0] i f = [1; 1; 1].

(f) Policz norm¦ drug¡ wektorow¡ pomi¦dzy x, a wynikiem al-gorytmu w, polegaj¡cym na zastosowaniu odpowiedniego rozkªadu oraz takie same normy residualne, tzn. normy dru-gie Aw − f oraz Rw − Q^Tf.

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

Rozpatrzmy macierz A2n,k - pod-macierz wymiaru 2n × k ma-cierzy Vandermonde'a A2n,2n dla 2n w¦zªów równoodlegªych na [0, 1].

LZNK z A2n,k z wektorem prawej strony f równym pierwszej kolumnie tej macierzy (rozwi¡zanie to pierwszy wersor) rozwi¡»

trzema sposobami:

(a) u»ywaj¡c operator \,

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

(c) poprzez rozwi¡zanie ukªadu równa« normalnych:

Bx = g dla

B = A^T_2n,kA_2n,k, g = A^T_2n,kf,

tzn. tworzymy macierz ukªadu równa« normalnych B, wek-tor prawej strony g ukªadu równa« normalnych, a nast¦pnie rozwi¡zujemy ukªad równa« normalnych operatorem \ . Macierz Vandermonde'a mo»na w octave'ie utworzy¢ za pomoc¡

funkcji vander().

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

• czas oblicze«

• bª¡d - kx − yk2

• bª¡d residualny kAx − fk2

dla x rozwi¡zania dokªadnego LZNK, f wektora prawej strony LZNK, y przybli»enia rozwi¡zania uzyskanego dan¡ metod¡.

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

liniowych,

Rozpatrzmy macierz An,n Vandermonde'a dla n w¦zªów równo-odlegªych na [0, 1].

Ukªad równa« liniowych z wektorem prawej strony równym pierw-szej kolumnie tej macierzy (rozwi¡zanie to pierwszy wersor) roz-wi¡» dwoma sposobami:

(a) operatorem \,

(b) rozkªadem QR uzyskanym funkcj¡ qr(), Przetestuj dla N = 10, 20, 40, 80. Porównaj

• czas oblicze«,

• bª¡d kx − yk2,

• bª¡d rezydualny: kAx − fk2,

dla x rozwi¡zania dokªadnego tego ukªadu równa«, f wektora prawej strony ukªadu i y przybli»enia rozwi¡zania uzyskanego dan¡ metod¡.

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

Zastosuj octave'a do znalezienia wspóªczynników a, b krzywej naj-lepiej pasuj¡cej do zadanych punktów: (xk, y_k), tzn. znajd¹ takie a, b, »e

X

k

|ax²_k+ by²_k− 1|² = min

c,d

X

k

|cx²_k+ dy_k²− 1|². Za (xk, y_k) przyjmiemy zaburzone punkty z danej elipsy

yk = q

1 − 4 ∗ x²_k+ zk,

gdzie zk to zaburzenie wylosowane z [0, 10⁻²] a xk = 1/k lub x_k = −1 + h ∗ k dla h = 2/N k = 1, ..., N.

Czy obliczone a i b jest bliskie 4 i 1?

W jednym oknie zaznacz punkty (xk, y_k) plusami oraz narysuj fragmenty wykresów obu elips: pierwszej - dla a = 4, b = 1 i drugiej elipsy - dopasowanej do zaburzonych punktów.

Powtórz obliczenia dla ró»nych zaburze« zk. Zadanie 6 Zaprogramuj funkcj¦ octave'a

function y=H(x ,w, nw)

która dla danych wektorów ~x i ~w tego samego wymiaru N i ska-laru nw = kwk2 zwróci wektor y = Hw~x dla

H_w = I − 2 ∗ 1 nww ~~w^T czyli przeksztaªcenia (macierzy) Householdera.

Skalar mo»e by¢ parametrem opcjonalnym. Je±li funkcja b¦dzie wywoªana z dwoma tylko parametrami, to norm¦ w mo»na obli-czy¢ w tej funkcji.

Przetestuj t¦ funkcj¦ dla losowych wektorów ~x i ~w i sprawd¹, czy kH_w~xk₂ = k~xk₂, H_w(H_wx) = x.

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

funkcj¦:

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

gdzie X macierz N × M i wtedy zwracany wynik to macierz Y = H_w ∗ X . Pozostaªe dwa parametry funkcji pozostan¡ bez zmian.

Czy mo»na zaimplementowa¢ tak¡ funkcj¦ bez u»ycia p¦tli?

Sprawd¹, wykorzystuj¡c t¦ funkcj¦, czy mno»enie przez macierz Householder nie zmienia norm macierzowych drugiej i Frobe-niusa, tzn. czy:

kAk₂ = kH_w ∗ Ak₂ = kA ∗ H_wk₂ i

kAk_F = kH_w∗ Ak_F = kA ∗ H_wk_F

dla losowej macierzy A i Hw macierzy Householdera dla losowego wektora w 6= 0.

Zadanie 8 Znajd¹ wektor Householdera ~w taki, »e odpowiednie przeksztaª-cenie Householdera przeprowadza dany wektor ~u 6= 0 na kierunek drugiego danego niezerowego wektora l ∗ ~v 6= 0 dla l = ^kvk_kuk²₂. Przetestuj dla dowolnych dwóch ró»nych wektorów o tej samej dªugo±ci, czy rzeczywi±cie Hw~u = ~v.

Zadanie 9 Zastosuj metod¦ Householdera do rozwi¡zania zadania znalezie-nia prostej y = ax+b najlepiej przybli»aj¡cej N punktów (xk, y_k) = (k, 1 + 2 ∗ k + _k) dla k = 1, . . . , N gdzie (1, . . . , _N) to losowy wektor za zakresu [−, ], tzn.:

X

k

|ax_k+ b − y_k|² = min

c,d

X

k

|cx_k+ d − y_k|².

Nale»y testowa¢ dla warto±ci  = [1, 10⁻¹, 10⁻², 10⁻³].

Funkcja rand(n) generuje wektor losowy o rozkªadzie jednostaj-nym na [0, 1] w octavie.

Porównaj z wynikami otrzymanymi za pomoc¡ standardowej funk-cji octave'a, tzn. \ , oraz przy wykorzystaniu funkfunk-cji octave'a qr(A).

Zadanie 10 Zaprogramuj metod¦ Householdera rozwi¡zywania ukªadu rów-na« liniowych A~x = ~b dla A macierzy trójdiagonalnej N × N, tzn. napisz funkcj¦ octave'a:

function [x]=hous3diag(a,b,c,f,N) Parametry funkcji:

• a, b, cprzek¡tna, pod-przek¡tna i nad-przek¡tna macierzy A,

• f - wektor prawej strony,

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

Funkcja zwraca x rozwi¡zanie Ax = f.

Przetestuj dziaªanie funkcji analogicznie do zadania 4 dla macie-rzy trójdiagonalnej o staªych diagonalach, np. takiej, »e elementy na gªównej diagonali s¡ równe dwa, a elementy na pod- i nad-diagonalach s¡ równe minus jeden dla N = 10^p z p = 1, 2, . . . , 9..

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

macierzy A.

Rozdziaª 11

Numeryczne zadanie wªasne

W tym rozdziale zajmiemy si¦ symetrycznym zadaniem wªasnym, tzn. za-daniem znajdowania warto±ci i/lub wektorów wªasnych dla macierzy syme-trycznej A = A^T. W zadaniach zbadamy, jak dziaªa podstawowa funkcja octave'a rozwi¡zuj¡ca zadanie wªasne, tzn. funkcja octave'a eig() oraz zba-damy zbie»no±¢ ró»nych wersji metody pot¦gowej.

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

x₀ = x₀ o normie drugiej równej jeden (np. losowego), deniujemy ci¡g iteracji metody pot¦gowej nast¦puj¡co:

ˆ

xk = Axk−1 (iteracja)

x_k = _kˆx^xˆk

kk2 (normowanie)

r_k = x^T_kAx_k (iloraz Rayleigha)

k > 0 (11.1) Zgodnie z teori¡, o ile x0 nie jest ortogonalny do wektora wªasnego (pod-przestrzeni wªasnej) dla najwi¦kszej warto±ci wªasnej, to x2kzbiegnie do tego wektora wªasnego, a rk - iloraz Rayleigha - do tej warto±ci wªasnej.

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

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





0 1 1 1 0 1 1 1 0



.

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

Wyznacz macierz V tak¡, »e jej kolumny to wektory wªasne A.

Przetestuj w normie macierzowej Frobeniusa, czy kAV − ΛV k2

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

W dokumencie Zadania i scenariusze zaj¦¢ z laboratorium komputerowego do wykªadu z Matematyki Obliczeniowej (Stron 33-76)