NumMod 1 pset 4

Cél: Az $f(x)=0$ egyenlet megoldása, ahol adott, szép tulajdonságú (pl. folytonosan deriválható) függvény, azaz keressük az függvény gyökeit/zérushelyeit. Ekvivalensen $g(x)=b$ megoldása az $f=g-b$ függvényre alkamazott gyökkereséssel.

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 é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+2x^2+10x-20$

egyenlet gyöke $x=\mathrm{1.368808107...}$ Ellenőrizzük le az intervallumfelezéses módszerrel, hogy jól számolt-e!
Megoldás: Legyen a toleranciánk $10^{-10}$ és indítsuk az iterációt az $[1,2]$ intervallumon!

Megoldás:
halving.py

def f(x):
	return x*x*x + 2*x*x + 10*x - 20
 
def halving(f, a, b, tolerance):
	assert f(a) * f(b) < 0
	if f(a) > 0:
		f = -f
 
	while a - b > tolerance:
		half = (a + b) / 2
		if f(half) < 0:
			a = half
		else:
			b = half
 
	return a, b
 
print(halving(f, 1, 2, 1e-10))

$ python3 halving.py
> (1.3688081077998504, 1.368808107858058)

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\le q<1$ együtthatóval, akkor van $x^{\ast }$ fixpontja, amire $f(x^{\ast })=x^{\ast }$; továbbá az $x_{n+1}=f(x_n)$ sorozat tetszőleges $X$-beli pontból indítva $x^{\ast }$-hoz fog tartani; továbbá $e_n=x_n-x^{\ast }$ választással látható, hogy

$\Vert e_{n+1}\Vert =\Vert x_{n+1}-x^{\ast }\Vert =\Vert f(x_n)-f(x^{\ast })\Vert \le q\Vert x_n-x^{\ast }\Vert =q\Vert e_n\Vert$
és
$\Vert e_n\Vert =\Vert (f\circ \dots \circ f)(x_0)-(f\circ \dots \circ f)(x^{\ast })\Vert \le q^n\Vert x_0-x^{\ast }\Vert =q^n\Vert e_0\Vert$

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 U$ olyan, hogy az $x+h\cdot [0,1]\subseteq U$, azaz az $x$-ből induló $x+h$ végű szakasz része $U$-nak; ekkor

$f(x+h)-f(x)={\int }_0^1{f}^{\prime }(x+th)\cdot hdt.$

Ebből kapjuk a középérték-tétel magasabbdimenziós változatát:

$\Vert f(x+h)-f(x)\Vert \le \left({\sup }_{0\le t\le 1}\Vert f^{\prime }(x+th)\Vert \right)\Vert h\Vert .$

Ha $K\subseteq U\subseteq X$ és $K$ kompakt, akkor

$M={\sup }_K\Vert f^{\prime }\Vert$

választással $x,y\in K$ esetén ha az $x$ és $y$ pontok által határolt szakasz $K$-ban van (pl. ha konvex), akkor

$\Vert f(x)-f(y)\Vert \le M\Vert x-y\Vert .$

Ha $f(K)\subseteq K$ (hívjuk ezt ráképezésnek) és $M<1$, akkor $f{|}_K:K\to K$ kontrakció. Mivel teljes tér zárt részhalmaza teljes, $K$ ezért teljes, 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 -nél kisebb.

Megoldás:
Kell egy Lipschitz fuggveny mely megadja a $\cos x=x$ egyenlet megoldasat.
Erre lehet hasznalni fixpont iteraciot, mert $\cos x$ kontrakcio $[0,\frac{\pi }{2})$ intervallumon

iteration.py

from math import cos
 
def fix_point_iteration(start, f, max_steps, tolerance):
	x = start
	step_counter = 1
	while step_counter < max_steps:
		change = f(x) - x
		x = f(x)
		if abs(change) < tolerance:
			break
 
		step_counter += 1
	return f"Ended after {steps} iteration with value = {x}"
 
print(fix_point_iteration(0, cos, 100, 1e-5))

$ python3 iteration.py
> Ended after 30 iteration with value = 0.7390822985224023

A lineáris esettel analóg módon bizonyos feltételek mellett fixpontiterációt készíthetünk az alábbi átalakítással:

$$\begin{array}{rl}f(x)& =b\\ 0& =b-f(x)\\ x& =x-f(x)+b\end{array}$$

Az iterációt ugyanolyan módon el tudjuk végezni mint a lineáris esetben, amíg $f(x)=x-(f(x)-b)$ kontrakció. Ha nem az, akkor valamely $c$ megválasztása mellett még lehet esélyünk arra, hogy az $f_c(x)=x-c(f(x)-b)$ függvény kontrakció legyen, a kérdés, hogy ezt hogy válasszuk meg.

Egy lehetséges választ kaphatunk erre a kérdésre az alábbi gradiens-módszerre vonatkozó tételben, mely analóg lesz a lineáris esetnél látottakkal, melyek szerint a Richardson iteráció tulajdonképpen egy állandó lépésközű Gradiens-ereszkedésnek felelt meg.

Gradiens-módszer

A gradiens-módszert is tudjuk alkalmazni nemlineáris feladatokra. A következő tételt alkalmazhatjuk speciális $A$ leképezések esetén. A tételt általános esetben mondjuk ki, egy $A:H\to H$ Hilbert-tér operátor esetén, az $A(u)=b$ egyenletre.

Tétel: Legyen $H$ valós Hilbert-tér (gondolhatunk $H={\mathbb{R}}^n$-re is az Euklideszi skalárszorzattal ellátva) és $A:H\to H$ legyen deriválható minden $u\in H$ pontban, azaz létezik az $A^{\prime }(u)\in B(H)$ lineáris operátor. Legyen $A^{\prime }(u)$ önadjungált minden $u\in H$ esetén. Tegyük fel, hogy léteznek $0<m\le M$ állandók, hogy

$m\Vert h{\Vert }^2\le ⟨A^{\prime }(u)h,h⟩\le M\Vert h{\Vert }^2\forall u,h\in H.$

Ekkor létezik olyan $\varphi :H\to \mathbb{R}$ funkcionál, amelyre teljesül, hogy ${\varphi }^{\prime }=A$, és a $\varphi$ funkcionál szigorúan konvex, továbbá, mint a lineáris esetben láthattuk:
$A(u^{\ast })=b$ pontosan akkor áll fenn, ha $(\varphi (u^{\ast })-⟨b,u^{\ast }⟩{)}^{\prime }=0$, azaz ha $u^{\ast }$ minimumhelye az $u\mapsto \varphi (u)-⟨b,u⟩$ leképezésnek. Ekkor tetszőleges $u_0$ esetén az

$u_{n+1}=u_n-\tau (A(u_n)-b)$

sorozat konvergál az $u^{\ast }$-hoz, $\tau =\frac{2}{m+M}$ mellett.

Megjegyzés: Feladattól függően nem állandó lépésköz is választható, ezt azonban most nem tárgyaljuk.

Ez a tétel valamilyen értelemben analóg a lineáris esetben tanultakkal: a feltételek miatt létezik az egyenletben adott monoton $A$ függvény “primitív függvénye”, ami konvex is, így az erre felírt minimalizási feladat az eredeti feladatunk egyértelmű megoldása lesz.

3. Feladat

Gondoljuk meg mit jelentenek a tétel feltételei $H={\mathbb{R}}^n$ esetén, illetve azt is hogy a képletben milyen dimenziójú objektumok szerepelnek. Gondoljuk meg, hogy a tétel visszaadja a lineáris esetnél látottakat.

Megoldás:

4. Feladat

Legyen $H={\mathbb{R}}^n$. Ellenőrizzük, hogy ha teljesülnek az előbbi tétel $A$-ra vonatkozó feltételei, akkor $F(u)=u-\tau (A(u)-b)$ kontrakció. Mutassuk meg, hogy bármely $0<\tau <\frac{2}{M}$ választás jó.

Megoldás:

$$F(u)=u-\tau (A(u)-b)$$

Kell:

$$||F(u)-F(v)||\le L\cdot ||u-v||,l\le 1$$

tudjuk:

$$||F(u)-F(v)||\le \underset{u\in {\mathbb{R}}^n}{\sup }||F^{\prime }(u)|{|}_{2,2}\cdot ||u-v|{|}_2$$ $$F^{\prime }(u)=I-\tau \cdot A^{\prime }(u)$$

$F^{\prime }(u)$ sajatertekei: $1-\tau \lambda (A^{\prime }(u))$, ahol $\lambda (A^{\prime }(u))$ az $A^{\prime }(u)$ sajaterteke.

$$|1-\tau \cdot \lambda (A^{\prime }(u))|=\max \{1-\tau \cdot \lambda (A^{\prime }(u)),\tau \cdot \lambda (A^{\prime }(u))-1\}\le \max \{1-\tau \cdot {\lambda }_{\min },\tau \cdot {\lambda }_{\max }-1\}$$ $$\le \max \{1-\tau \cdot m,\tau \cdot M-1\}$$

Kell, hogy ez $<1$ legyen, ehez eleg hogy

$$\begin{array}{rl}1-\tau \cdot m& <1⟺0<\tau \cdot m\\ \tau \cdot M-1& <1⟺\tau <\frac{2}{M}\end{array}$$

Tehat belattuk a feladat allitasat, hogy $0<\tau <\frac{2}{M}$ valasztas jo.

5. Feladat

Legyen $A:{\mathbb{R}}^n\to {\mathbb{R}}^n$ olyan függvény amelyre léteznek valamely $B,C\in {\mathbb{R}}^{n\times n}$ invertálható mátrixok, amelyekre $B\cdot A^{\prime }(x)\cdot C$ SZPD mátrix minden $x\in K\subset {\mathbb{R}}^n$ esetén, egy alkalmas $K$ esetén, ahol $A(K)\subset K$ (vagyis $A$ ráképezés $K$-n). Ennek segítségével hogyan tudjuk megoldani az $A(x)=b$ feladatot iterációval?

Próbáljuk megoldani az $\{\begin{array}{ll}-x(6x^2+2y^2)& =2\\ y(3x^2+y^2)& =-2\end{array}$

feladatot.

Megoldás:

6. Feladat

Keressük meg az $f(x)=x^3+2x^2+10x-20$ egyenlet gyökét, ha valahonnan sejtjük, hogy 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.

Megoldás:

7. Feladat

Tekintsük az alábbi függvényt, ahol $c\in {\mathbb{R}}^{+}$:

$f(x)=\frac{1}{x}-\frac{x}{c}$

Keressük a függvény $x^{\ast }>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?

Megoldás:

Kiegészítés lineáris egyenletrendszerek megoldására: Konjugált gradiens módszer

Emlék előadásról: A gradiens módszer helyett, gyakran a konjugált gradiens módszert használjuk, ha lineáris egyenletendszert szeretnénk megoldani. A módszer gyorsabban konvergál és elméletben véges lépés után visszakapjuk a pontos megoldást. A módszer lépései:

1. Adott a kezdeti megoldás $x_0$ és a kezdeti gradiens $p_0=r_0=A{x}_0-b$, ahol $A$ a mátrix és $b$ a vektor a lineáris egyenletrendszerben.
2. Ha $x_n$ és $p_n$ ismert

$$x_{n+1}=x_n-{\alpha }_n{p}_n,{\alpha }_n=\frac{⟨r_n,p_n⟩}{⟨A{p}_n,p_n⟩},r_{n+1}=A{x}_{n+1}-b$$ $$p_{n+1}=r_{n+1}-{\beta }_n{p}_n,{\beta }_n=\frac{⟨A{r}_{n+1},p_n⟩}{⟨A{p}_n,p_n⟩}$$

8. Feladat

Írjunk programot, amely megold egy egyenletrendszert a konjugált gradiens módszerrel, amelynek bemenetei $A$, $b$ és az iterácók száma, majd oldjuk meg az alábbi egyenletet. Próbáljuk ki úgy is, hogy $1,2$, illetve $3$ iterációs lépést teszünk.

$$\{\begin{array}{l}4x+2y=7,\\ 2x+3y=10.\end{array}$$

Megoldás:

related: NumMod 1