%% Hetedik gyakorlat
%% Nemlineáris egyenletek: intervallumfelezés, fixpont iteráció
%
% Cél: Az f(x)=0 egyenlet megoldása, ahol f adott, szép tulajdonságú (pl. folytonosan deriválható) függvény, azaz keressük az f függvény gyökeit/zérushelyeit. Ekvivalensen g(x) = b megoldása az f = g-b függvényre alkamazott gyökkereséssel.
% Megjegyzés: A gyökkeresési feladat (ekvivalensen: egyenlet(rendszer) megoldás) hibaérzékenységét nemlineáris függvények esetén is jellemezhetjük egy kondíciószámmal. Egydimenziós esetben, tehát egy f: \mathbb{R} \rightarrow \mathbb{R} függvényre vonatkozólag, egy adott [a,b] intervallumon a zérushelykeresési feladat kondíciószáma legyen
% \kappa = \frac{1}{\inf_{x \in [a,b]}|f'(x)|}.
% Ha ez az érték nagy, akkor a feladat rosszul kondícionált (azaz az f függvény kis megváltoztatásával a megoldás is nagyon megváltozhat).
%% I. Intervallumfelezés
% Tegyük fel, hogy tudjuk, hogy egy adott [a,b] intervallumon a függvényünknek van egy gyöke; például ezt onnan tudhatjuk, hogy a végpontokban különböző előjelű a függvény, pontosabban:f(a) \cdot f(b) < 0. Tegyük fel, hogy ez fennáll és hogy f(a) > 0, f(b)<0. Ez utóbbi nem megszorítás, mert f és -f gyökei ugyanott vannak.
% Alapgondolat: Vizsgáljuk meg a függvény előjelét az \frac{a+b}{2} pontban! Ha itt pozitív, az azt jelenti, hogy a keresett zérushely az \left[a,\; \frac{a+b}{2}\right) intervallumban van, azaz itt érdemes folytatni a keresést. Ugyanakkor ha itt negatív, akkor a zérushely az \left( \frac{a+b}{2},\; b\right] intervallumban található, és itt kell tovább keresnünk. Ezt követően az előbbi gondolatot folytatjuk az új, kisebb intervallumokon, azaz vagy az \frac{a+b}{4}, vagy a \frac{3 (a+b)}{4} pontban vizsgáljuk a előjelet, és így tovább.
% Az algoritmus addig fut, amíg f értéke a vizsgált pontban nulla nem lesz (ez programok esetén a gépi hibák miatt ritkán következik be), vagy amikor a vizsgált intervallum már nagyon kicsi/rövid, de természetesen megadhatunk maximális lépésszámot is, amely után álljon le mindenképpen a program. Határértékben ez a módszer mindenképpen megtalál egy zérushelyet (ha több is van, akkor a zérushelyek közül valamelyiket).
% 1. feladat: Fibonacci egyik cikkében szerepel az alábbi állítás: az
% f(x) = x^3 + 2 x^2 + 10 x - 20
% egyenlet gyöke x = 1.368808107... Ellenőrizzük le az intervallumfelezéses módszerrel, hogy jól számolt-e!

%

%% II. Fixpont-iteráció
% Emlék: A Banach-féle fixponttétel következtében ha X Banach-tér és f: X \to X kontrakció 0\leq q < 1 együtthatóval, akkor van x^* fixpontja, amire f(x^*) = x^*; továbbá az x_{n+1} = f(x_n) sorozat tetszőleges X-beli pontból indítva x^*-hoz fog tartani; továbbá
% e_n = x_n - x^* választással látható, hogy
% \|e_{n+1}\| = \| x_{n+1} - x^* \| =  \| f(x_{n}) - f(x^*) \| \leq q \|x_n - x^*\| = q \| e_n\|
% és
% \|e_n\| = \|  (f\circ \ldots \circ f)(x_0) - (f\circ \ldots \circ f)(x^*) \| \leq q^n\|x_0 - x^*\| = q^n\|e_0\|.
% Emlék: Legyen X végesdimenziós valós vektortér, U \subseteq X nyílt halmaz, f: U \to X folytonosan differenciálható és legyen x, h \in Uolyan, hogy az x+h\cdot[0,1] \subseteq U, azaz az x-ből induló x+h végű szakasz része U-nak, akkor
% f(x+h) - f(x) = \int_{0}^1 f'(x+th) \cdot h\, dt.
% Ebből kapjuk a középérték-tétel magasabbdimenziós változatát:
% \|f(x+h) - f(x)\| \leq \left(\sup_{0\leq t \leq 1}\|f'(x+th)\| \right) \|h\|.
% Ha K \subseteq U\subseteq X és Kkompakt, akkor
% M = \sup_{K}\|f'\|
% választással x, y \in K esetén ha az x és y pontok által határolt szakasz K-ban van (pl. ha K konvex), akkor
% \| f(x) - f(y) \| \leq M \|x -y\|.
% Ha f(K) \subseteq K (hívjuk ezt ráképezésnek) és M <1, akkor f|_K:K \to Kkontrakció. Mivel teljes tér zárt részhalmaza teljes, ezért Kteljes, tehát f|_K-ra is alkalmazható a Banach-féle fixponttétel. Érdemes lehet megjegyezni, hogy K általában nem vektortér, azonban az eredeti tér normája által generált metrikával a halmaz teljes metrikus tér, tehát a Banach-féle fixponttétel valóban alkalmazható.
% 2. feladat:
% Oldjuk meg a \cos(x) = x egyenletet a valós számok halmazán. Keressünk alkalmas halmazon alkalmas kontrakciót, majd írjunk kódot amivel addig iteráljunk, míg két szomszédos lépés távolsága nem lesz 10^{-4-nél kisebb.

%

% Sajnos az élet nem mindig ad ilyen szép egyenleteket nekünk. Például \cos(6x) = x választással a -6\sin(x) deriváltra már nem adható ilyen egyszerűen jó felső becslés. Általános módszer nincs, de nézzünk meg egy feladatosztályt, amivel el tudunk bírni.
% Monoton (baloldalú) egyenletek
% Tekintsük az g(x) = b feladatot, amire
% g(x) = b \iff f(x) = g(x) - b  = 0,
% sőt c \neq 0 konstans esetén a lánc folytatható: 0 = cf(x) \iff x = x - c f(x) = F_c(x).
% Keressünk egy feltételt, ami mellett valamilyen c választással F_c-re alkalmazható a fixpont-iteráció.
% Tudjuk, hogy F'(x) = I - c f'(x), így ha A = cf'(x), akkor
% \| I- A\|_2^2 = \sup_{\|x\|_2 = 1}x^T(I - A^T)(I - A)x = \sup_{\|x\|_2 = 1}x^T (I - (A+A^T)+ A^TA)x = 1 + \sup_{\|x\|_2 = 1}x^T (-(A+A^T)+ A^TA)x.
% Emiatt
% \|F'(x)\|_2^2 < 1 \iff \sup_{\|x\|_2 = 1}x^T (-(A+A^T)+ A^TA)x < 0,
% átfogalmazva, minden egységhosszú x vektorra
% x^T( -(A+A^T) + A^TA)x  < 0.
% Tegyük fel, hogy A szimmetrikus. EkkorA = Q \Lambda Q^T diagonalizációval és y = Q^Tx helyettesítéssel (ami továbbra is lefedi az összes lehetséges egységvektort)
% -2 y^T \Lambda y + y^T \Lambda^2 y  < 0,
% azaz
% y^T \Lambda(\Lambda - 2I) y < 0.
% Ez nem jelent mást, mint hogy a \Lambda(\Lambda - 2I) mátrix sajátértékei negatívak, szóval minden \lambda sajátértékre \lambda(\lambda - 2) < 0.
% Emélkezzünk, hogy A = cf'(x) volt, tehát implicit feltettük, hogyf'(x) szimmetrikus. Legyen f' szimmetrikus egy K halmaz minden pontjában, továbbá ott a \mu sajátértékei abszolútértékben ne legyenek nagyobbak egy adott M konstansnál. Ekkor tehát minden \mu sajátértékre c\mu(c\mu -2) = c(c\mu^2 - 2\mu) < 0 kellene, hogy teljesüljön.
% Ez a kifejezés c-ben másodfokú, konvex polinom, mely pontosan a gyökei között negatív.
% A gyökök: c_1 = 0  és c_2 = \frac{2\mu}{\mu^2} = \frac{2}{\mu}.
% Láthatjuk, hogy ha van két ellentétes előjelű sajátérték, akkor ez a megközelítés nem működhet (mert c az egyik miatt negatív, a másik miatt pedig pozitív kellene, hogy legyen), tehát tegyük fel, hogy minden sajátérték azonos előjelű (azaz f' a K halmaz minden pontjában (poz. vagy neg.) definit).
% Látható, hogy c_2 a nullához legközelebb akkor van, ha \mu a lehető legnagyobb abszolútértékű (legyen ez az abszolútérték M).
% Ha a sajátértékek negatívak, akkor -\frac{2}{\mu}\leq -\frac{2}{M} < c < 0 jó választás. Ha a sajátérékek pozitívak, akkor 0 < c < \frac{2}{M} \leq \frac{2}{\mu} jó választás.
% Ezekben az esetekben a derivált hosszára vonatkozó felső becslés alapján kontrakciót kapunk az adott K halmazon. Ha F_c(K) \subseteq K és Kmegfelelő (pl. konvex), akkor F_c: K\to Kkontrakció és lesz fixpontja K-ban.
% Összefoglalva: Legyen f egy K zárt, konvex halmazon monoton (azaz f' ott legyen minden pontban (poz. vagy neg.) definit).
% Ha f'szimmetrikus K minden pontjában és a sajátértékeinek abszolútértékére van egy M felső korlát, akkor 0 < c < \frac{2}{M} választással
% ha f poz. monoton K-n, akkor az F_c(x) = x - cf(x),
% ha f neg. monoton K-n, akkor az F_c(x) = x + cf(x)
% függvény kontrakció K-n. Ha ez K-t önmagába is képezi (ráképezés), akkor K \to K kontrakció, tehát K-ban lesz fixpontja és ehhez tartani is tudunk egy K-beli elemből indulva.
% 3. feladat:
% Legyen g(x) = Ax valamilyen poz. def A mátrixra. Tegyük fel, hogy a legnagyobb sajátértéke kisebb, mint egy ismert M konstans. Oldjuk meg az g(x) = b egyenletet fixpont-iterációval. Addig iteráljunk, míg két szomszédos iterált távolsága nem lesz kisebb, mint 10^{-4}.
% Például, vehetjük az alábbi paramétereket:

[Q, R] = qr(magic(3));
A = Q'*diag([1,2,3])*Q
b = ones(3, 1);
M = 3.1;

% de persze lehet mást is.

%

% 4. feladat
% Keressük meg az f(x) = x^3 + 2 x^2 + 10 x - 20 egyenlet gyökét, ha valahonnan sejtjük, hogy az a K=[1,2] intervallumban van. Alkalmazzunk fixpont-iterációt és addig iteráljunk, míg a szomszédos iteráltak távolsága 10^{-10} alá nem csökken.

%

% Megjegyzés:  Az f(x) = 0 egyenletet természetesen nem csak a korábbiakhoz hasonlóan lehet manipulálni.
% Például, ha f:\mathbb{R} \to \mathbb{R}, f(x) = g(x) x - h(x) alakú, ahol g egy függvény, melynek nincs gyöke K=[a,b]-ban, akkor
% f(x) = 0 \iff x  = \frac{h(x)}{g(x)} = F(x) a Khalmazon.
% Ha F kontrakció K-n és F(K)\subseteq K, akkor alkalmazhatjuk rá a fixpont-tételt a fixpont (azaz f zérushelyének) meghatározásához.
% 5. feladat: Alkalmazzuk ezen megjegyzés ötletét az első feladatban szereplő
% f(x) = x^3 + 2 x^2 + 10 x - 20
% kifejezés [1,2] intervallumban lévő gyökének meghatározására.

%

% 6. feladat: Határozzuk meg azon pontot, ahol az f(x)=x^5 - x függvény értéke 2 lesz. Oldjuk meg a feladatot
% intervallumfelezéssel, majd
% fixpont-iterációval!

%

%% Konvergenciarend
% Definíció: Legyen az x_k \rightarrow x^* konvergens iteráció k-adik elemének hibája e_k (azaz e_k = x_k - x^*). Ekkor azt mondjuk, hogy az iteráció konvergenciarendje p, ha
% p = \lim_{k \rightarrow \infty} \frac{\ln ||e_k||}{\ln ||e_{k-1}||} =\lim_{k \rightarrow \infty} \log_{\|e_{k-1}\|} \|e_k\|,
% tehát határértékben egy hibavektor hosszának durván hányadik hatványa a következő hibavektor hossza.
% Emlék: Fixpont-iterációra a szokásos feltételek mellett
% \|e_{n+1}\| = \| x_{n+1} - x^* \| =  \| F(x_{n}) - F(x^*) \| \leq q \|x_n - x^*\| = q \|e_n\|,
% ahol q az F kontrakció együtthatója. Egyéb feltételek nélkül tehát ez a módszer (legalább) elsőrendű.
% Tétel: Fixpont-iteráció esetén ha F(x_k) \rightarrow x^* és azt is tudjuk, hogy
% F'(x_k) = F''(x_k) = \dots = F^{(p-1)}(x_k)=0 \qquad \text{de} \qquad F^{(p)}(x^*)\neq 0,
% akkor a módszer p-ed rendben konvergál.
% Motiváció:
% e_{k+1} = x_{k+1} - x^* = F(x_k) - F(x^*) = \sum\limits_{j=1}^p \frac1{j!}F^{(j)}(x^*)(x_{k} - x^*)^j + o(\|x_{k} - x^*\|^{p}) =  0+ \frac1{p!}F^{(p)}(x^*) e_k^p + o(\|e_k\|^{p}),
% így
% \|e_{k+1}\| \leq C_p \|e_k\|^p  + o(\|e_k\|^p),
% tehát határértékben \|e_{k+1}\| durván (legalább) p-edik hatványa az \| e_k\| értéknek, azaz p-edrendről beszélhetünk.
% 4. feladat: Tekintsük a következő függvényt: f(x) = \frac{x}{x-1}.
% Hol veszi fel a 2 értéket? Határozzuk meg ezt a pontot kézzel!
% Mit tudunk mondani a konvergenciájának rendjéről, ha fixpont-iterációval közelítenénk?
% Alkalmas intervallum választásával numerikusan ellenőrizzük is a konvergenciarendet!

%

%
% Szorgalmi: Tekintsük az alábbi függvényt, ahol c \in \mathbb{R}^+:
% f(x) = \frac{1}{x} - \frac{x}{c}
% (a) Keressük a függvény x^*>0 zérushelyét fixpont iteráció segítségével! Írjuk fel alkalmas F(x)  függvényt, és bizonyítsuk be, hogy valóban konvergál minden c>0 esetén! Hová fog tartani?

%

% (b) Mit mondhatunk az iteráció rendjéről c függvényében? Az állításunkat numerikusan is szemléltessük!

%

% Szorgalmi:
% Oldjuk meg az
% \begin{cases}
% x(x^2+y^2) &= 2\\
% y(x^2+y^2) &= -2
% \end{cases}
% egyenletrendszert (síkban haladó) fixpont-iterációval.
% Határozzuk meg a megoldást kézzel.
% Keressünk a megoldást tartalmazó egyszerű alakú, alkalmas Khalmazt, amin a megfelelő függvény monoton. Ezt úgy tudjuk ellenőrizni, hogy kiszámoljuk a deriváltját, majd a derivált sajátértékeit, melyeknek azonos előjelűeknek (és nem nullának) kell lenniük.
% írjunk fel egy alkalmas K esetén egy alkalmas c választás mellett egy F_c függvényt, ami K-t önmagába képezi (ráképezés) és azon kontrakció. A ráképezést nem muszáj ellenőrizni.
% Futtassuk a fixpont-iteráció kódját az adott függvényre egy alkalmas kezdőpontból indítva (ez ne az analitikus megoldás legyen)

%

%
%% Függvénytemető:

%