{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "Řešení Poissonovy rovnice metodou konečných diferencí (metodou sítí)\n", "==============================================================\n", "\n", "### Definice úlohy\n", "\n", "Chceme vyřešit 2D Poissonovu rovnici\n", "$$\n", "\\frac{\\partial^2 u}{\\partial x^2} + \\frac{\\partial^2 u}{\\partial y^2} = \\rho\n", "$$\n", "na obdélníkové oblasti\n", "$$\n", "\\Omega = \\left<0, x_{\\rm max}\\right> \\times \\left<0, y_{\\rm max}\\right>\n", "$$\n", "s Dirichletovou okrajovou podmínkou\n", "$$\n", "u(x, y) = u_0(x, y);\\quad (x,y)\\in \\partial\\Omega.\n", "$$\n", "Pozn.: faktor $-1/\\epsilon_0$ z Poissonovy rovnice v elektrostatice zde pro jednoduchost a bez újmy na obecnosti vynecháváme." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Diskretizace\n", "Provedeme následovnou diskretizaci prostorových souřadnic:\n", "$$\n", "x_i = hi;\\quad i = 0,\\ldots,M+1;\\quad hM = x_{\\rm max},\n", "$$\n", "$$\n", "y_j = hj;\\quad j = 0,\\ldots,N+1;\\quad hN = y_{\\rm max}.\n", "$$\n", "Velikost kroku $h$ pro jednoduchost volíme stejnou v $x$-ovém i $y$-ovém směru. Celkový počet neznámých označíme $N_g = MN$. Numerickou aproximaci hledaného řešení $u(x, y)$ v bodě $(x_i, y_j)$ označíme \n", "$$\n", "u_{ij} \\approx u(x_i, y_j)\n", "$$\n", "a obdobné značení zavedeme i pro ostatní veličiny, tedy např $\\rho_{ij} \\approx \\rho(x_i, y_j)$.\n", "Na konec ještě parciální derivace nahradíme druhými diferencemi a dostaneme soustavu rovnic typu\n", "$$\n", "\\frac{u_{i-1,j} - 2u_{ij} + u_{i+1,j}}{h^2} +\n", "\\frac{u_{i,j-1} - 2u_{ij} + u_{i,j+1}}{h^2} = \\rho_{ij}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Abychom s touto soustavou lineárních rovnic mohli dále pracovat, je vhodné ji zapsat v maticovém tvaru. Pro jednoduchost začněme v 1D. První rovnice má tvar\n", "$$\n", "u_0 - 2u_1 + u_2 = h^2\\rho_1\n", "$$\n", "a vzhledem k tomu, že $u_0$ je dané Dirichletovou okrajovou podmínku, převedeme jej na druhou stranu\n", "$$\n", "-2u_1 + u_2 = h^2\\rho_1 - u_0\n", "$$\n", "rovnice pro $i=2,\\ldots,M-1$ budou\n", "$$\n", "u_{i-1} - 2u_{i} + u_{i+1} = h^2\\rho_i\n", "$$\n", "a konečně poslední rovnice pro $i=M$ bude opět obsahovat okrajovou podmínku\n", "$$\n", "u_{M-1} - 2u_{M} = h^2\\rho_i - u_{M+1}.\n", "$$\n", "Maticový zápis těchto rovnic bude\n", "$$\n", "\\left[\\begin{matrix}{}\n", " -2 & 1 & & &\\\\\\\n", " 1 & -2& 1 & &\\\\\\\n", " & .& . & .&\\\\\\\n", " & & 1& -2 & 1\\\\\\\n", " & & & 1 & -2\n", " \\end{matrix}\\right]\n", " \\cdot\n", " \\left[\\begin{matrix}{}\n", " u_1\\\\\\\n", " u_2\\\\\\\n", " \\vdots\\\\\\\n", " u_{M-1}\\\\\\\n", " u_M\n", " \\end{matrix}\\right]\n", " =\\left[\\begin{matrix}{}\n", " h^2\\rho_1-u_0\\\\\\\n", " h^2\\rho_2\\\\\\\n", " \\vdots\\\\\\\n", " h^2\\rho_{n-1}\\\\\\\n", " h^2\\rho_n-u_{M+1}\n", " \\end{matrix}\\right],\n", "$$\n", "kde matici na levé straně nazýváme _matice druhých diferencí_ a značíme $K$. Dále označíme vektor neznámých $U$ a pravou stranu $F$, čímž dostaneme maticový tvar soustavy\n", "$$\n", "KU=F.\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Maticový zápis ve 2D\n", "\n", "### Metody řešení\n", "#### Gaussova eliminace\n", "#### Relaxační metody\n", "\n", "\n", "##### Jacobi\n", "Z diskretizovaného zápisu Poissonovy rovnice\n", "$$\n", "u_{i-1} - 2u_{i} + u_{i+1} = h^2\\rho_i\n", "$$\n", "si můžeme vyjádřit $i$-tou složku řešení\n", "$$\n", "u_i =\\frac{u_{i-1} + u_{i+1}}{2} - \\frac{h^2\\rho_i}{2}\n", "$$\n", "a interpretovat tuto rovnici jako iterativní vztah pro výpočet $n+1$ aproximace řešení\n", "$$\n", "u_i^{n+1} =\\frac{u_{i-1}^n + u_{i+1}^n}{2} - \\frac{h^2\\rho_i}{2}\\qquad\\text{(eq:Jacobi)}\n", "$$\n", "##### Gauss Seidel\n", "Pokud postupně pro $i=1,\\ldots,M$ vypočítáváme $u_i$ z Jacobiho iteračního vztahu, vidíme, že v době výpočtu $u_i^{n+1}$ již známe $u_{i-1}^{n+1}$ a můžeme tedy ve výpočtu použít tuto aktuálnější hodnotu. Náš iterační vztah tedy bude mít tvar \n", "$$\n", "u_i^{n+1} =\\frac{u_{i-1}^{n+1} + u_{i+1}^n}{2} - \\frac{h^2\\rho_i}{2}.\\qquad\\text{(eq:GS)}\n", "$$\n", "Výhodou tohoto uspořádání je, že nové hodnoty můžeme přímo přepisovat v paměti, čímž snížíme paměťové nároky, a zároveň tato metoda rychleji konverguje.\n", "\n", "\n", "Abychom zjistili rychlost konvergence výše uvedených relaxačních metod, zapíšeme si je v maticovém tvaru.\n", "Obecně můžeme relaxační metody pro řešení soustav lineárních rovnic interpretovat tak, že matici soustavy $K$ nahradíme maticí $P$, kterou lze narozdíl od $K$ snadno invertovat, a přitom je $P$ je v určitém smyslu blízké $K$. Význam této blízkosti odvodíme níže. Matice $P$ se nazývá _předpodmiňovací_ (preconditioner).\n", "\n", "Při odvození obecné iterační metody pro soustavu\n", "$$KU=F$$\n", "si nejprve přičteme na obě strany výraz $PU$ a obtížně invertovatelný součin $KU$ převedeme na pravou stranu:\n", "$$PU = (P-K)U + F.$$\n", "Pokud tento výraz interpretujeme jako iterační vztah, kde na pravé straně je přechozí aproximace $U^n$ a na levé straně je nová aproximace $U^{n+1}$, můžeme $U^{n+1}$ snadno vypočíst, neboť jsme zvolili $P$ snadno invertovatelné:\n", "$$U^{n+1}= P^{-1}(P-K)U^n + P^{-1}F = U^n- P^{-1}(KU^n-F).\\qquad\\text{(eq:iter)}$$\n", "Chceme-li vědět, zda tato posloupnost $U^n$ konverguje k hledanému řešení $U$, napišme si $U^n$ jako součet hledaného řešení a chyby $e^n$ v daném kroku, $U^n = U + e^n$,\n", "a dosazením tohoto výrazu do iteračního vztahu (eq:iter) dostaneme\n", "$$\n", "U + e^{n+1}= P^{-1}(P-K)(U + e^n) + P^{-1}F.\n", "$$\n", "Částečně roznásobíme:\n", "$$\n", "U + e^{n+1}=U - P^{-1}\\underbrace{KU}_F + P^{-1}(P-K)e^n + P^{-1}F,\n", "$$\n", "odečteme $U$ od obou stran rovnice a využijeme, že pro $U$ platí $KU=F$ z definice, čímž dostaneme\n", "$$\n", "e^{n+1}=\\underbrace{P^{-1}(P-K)}_Me^n,\n", "$$\n", "kde matici $\\mathbf M=P^{-1}(P-K)$ nazýváme iterační maticí. Rozložíme-li si chybu počátečního odhadu $e^0$ do báze vlastních vektorů matice $\\mathbf M$: $e^0 = \\sum_{k=1}^M e_k^0 z^k$ tak chybu v $n$-té iteraci můžeme psát jako\n", "$$\n", "e^n = \\sum_{k=1}^M\\lambda_k^n e_k^0 z^k.\n", "$$\n", "Velikost chyby konverguje k nule, pokud velikost všech vlastních čísel je menší než jedna. Definujeme-li spektrální poloměr $\\rho(\\mathbf M)$ jako maximální velikost vlastního čísla matice $\\mathbf M$,\n", "$\\rho(\\mathbf M) = \\max_k\\{|\\lambda_k(\\mathbf M)|\\}$, můžeme podmínku konvergence zapsat jako\n", "$$\n", "\\rho(\\mathbf M) < 1.\n", "$$\n", "Onen požadavek na \"blízkost\" matic $P$ a $K$, který jsme vyjádřili při odvozování iterační metody je kvantifikován právě spektrálním poloměrem matice $\\mathbf M$, který můžeme v určitém smyslu interpretovat jako jejich vzdálenost.\n", "\n", "K výpočtu spektrálního poloměru iteračních matic pro řešení Poissonovy rovnice se nám bude hodit znát vlastní čísla matice druhých diferencí, která odvodíme v následující vsuvce.\n", "\n", "### Vlastní čísla a vektory matice druhých diferencí\n", "TODO odvození\n", "$$\n", "z_i^k = \\sin\\frac{ik\\pi}{M+1}\\qquad\\rm(eq:eigvec)\n", "$$\n", "$$\n", "\\lambda_k = 2\\cos\\frac{k\\pi}{M+1} - 2\\qquad\\rm(eq:eigval)\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Konvergence Jacobiho metody\n", "V jacobiho metodě volíme jako preconditioner diagonální část matice druhých diferencí, kterou označíme $D$. Maticový zápis tedy má tvar\n", "$$\n", "U^{n+1} = D^{-1}(D-K)U^n + D^{-1}F.\n", "$$\n", "Když uvážíme, že násobení čtvercovou diagonální maticí s konstantními koeficienty můžeme nahradit násobením skalárem a výsledné matice rozepíšeme, dostáváme\n", "$$\n", "U^{n+1} = \\frac{1}{2}\\left[\\begin{smallmatrix}{}\n", " 0 & 1 & & &\\\\\\\n", " 1 & 0& 1 & &\\\\\\\n", " & .& . & .&\\\\\\\n", " & & 1& 0 & 1\\\\\\\n", " & & & 1 & 0\n", " \\end{smallmatrix}\\right]U^n - \\frac{1}{2}F,\n", "$$\n", "což skutečně odpovídá maticovému zápisu Jacobiho metody (eq:Jacobi). Nyní potřebujeme vypočítat vlastní čísla Jacobiho iterační matice $\\mathbf M = D^{-1}(D-K)$. Vzhledem k tomu, že vlastní vektory $K$ jsou zároveň vlastními vektory $D$, lze pro vlastní číslo $\\lambda_k^{\\mathbf M}$ matice $\\mathbf M$ odpovídající vlastnímu vektoru $z^k$ psát\n", "$$\n", "\\lambda_k^{\\mathbf M} = \\frac{1}{\\lambda_k^D}(\\lambda_k^D-\\lambda_k^K)\n", "$$\n", "a uvážíme-li dále, že všechna vlastní čísla $\\lambda_k^D$ jsou rovna $-2$, dostáváme\n", "$$\n", "\\lambda_k^{\\mathbf M} = -\\frac{1}{2}\\left(-2-(2\\cos\\frac{k\\pi}{M+1} - 2)\\right) = \\cos\\frac{k\\pi}{M+1}.\n", "$$\n", "Pro spektrální poloměr potom platí\n", "$$\n", "\\rho(\\mathbf M) = \\max_{k=1,\\ldots,M}\\{|\\lambda_k^{\\mathbf M}|\\} = \\lambda_1^{\\mathbf M} = -\\lambda_M^{\\mathbf M} =\n", "\\cos\\frac{\\pi}{M+1} < 1,\n", "$$\n", "čímž jsme dokázali konvergenci Jacobiho metody. Zároveň vidíme, že největší vlastní čísla odpovídají nejnižším a nejvyšším prostorovým frekvencím a tyto složky chybového vektoru jsou tedy redukovány nejpomaleji.\n", "\n", "**Př:** Kolik iterací Jacobiho metody je potřeba pro redukci chyby faktorem $10^{-p}$?\n", "Chyba po $n$ iteracích se redukuje faktorem $\\rho^n(\\mathbf M)$, takže řešíme rovnici\n", "$$\n", "10^{-p} = \\rho^n(\\mathbf M).\n", "$$\n", "Spektrální poloměr si můžeme aproximovat Taylorovým rozvojem \n", "$\\rho(\\mathbf M) = \\cos\\frac{k\\pi}{M+1}\\approx 1 - \\frac{1}{2}\\left(\\frac{\\pi}{M+1}\\right)^2\\approx1 - \\frac{1}{2}\\left(\\frac{\\pi}{M}\\right)^2$ a po zlogaritmování rovnice dostáváme\n", "$$\n", "-p\\ln10 = n\\ln\\left(1 - \\frac{1}{2}\\left(\\frac{\\pi}{M}\\right)^2\\right).\n", "$$\n", "Logaritmus můžeme opět aproximovat Taylorovým rozvojem jako\n", "$\\ln\\left(1 - \\frac{1}{2}\\left(\\frac{\\pi}{M}\\right)^2\\right)\\approx- \\frac{1}{2}\\left(\\frac{\\pi}{M}\\right)^2$\n", "a můžeme vyjádřit $n$ jako\n", "$$\n", "n = 2p\\ln10\\cdot\\frac{M^2}{\\pi^2}\\approx0.4666\\cdot pM^2\\approx\\frac{1}{2}pM^2.\n", "$$\n", "Pozn.: v $d$ dimenzích je počet iterací úměrný $N_g^{2/d}$, tedy druhé mocnině lineárního rozměru mříže.\n", "Výpočetní náročnost jedné iterace je $O(N_g)$ a celková výpočetní náročnost bude $O(N_g^{1+2/d})$. Konvergence je tedy pomalá. Určitou výhodou je, že počet iterací neroste s dimenzí mříže (ovšem počet operací ano, neboť jednotlivé iterace jsou výpočetně náročnější)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Konvergence Gauss Seidelovy metody\n", "Uvedeme pouze stručně bez odvození. Předpodniňovací matice Gauss Seidelovy metody je dolní trojúhelníková (včetně diagonály), $P=L$. Inverzi této matice lze snadno řešit dopřednou substitucí. Vlastní čísla iterační matice Gauss Seidelovy metody jsou druhou mocninou vl. čísel matice Jacobiho metody,\n", "$$\n", "\\rho_{\\rm GS}=\\rho_{\\rm Jac}^2\\approx 1-\\left(\\frac{\\pi}{M}\\right)^2.\n", "$$\n", "Potřebný počet iterací je tedy poloviční,\n", "$$\n", "n\\approx \\frac{1}{4}p M^2.\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Superrelaxace (Successive overrelaxation, SOR)\n", "Myšlenka superrelaxace obecně spočívá v tom, že korekci řešení předepsanou nějakým iteračním vztahem\n", "$$\n", "x^{n+1} = f(x^n) = x^n + (f(x^n)-x^n)\n", "$$\n", "vynásobíme faktorem $\\omega$ větším než 1:\n", "$$\n", "x^{n+1} = x^n + \\omega(f(x^n)-x^n) = (1-\\omega)x^n + \\omega f(x^n)\n", "$$\n", "Konkrétně pro Gauss Seidelovu metodu tak dostáváme\n", "$$\n", "u_i^{n+1} = (1-\\omega)u_i^n + \\omega\\left(\\frac{u_{i-1}^{n+1} + u_{i+1}^n}{2} - \\frac{h^2\\rho_i}{2}\\right).\n", "$$\n", "Optimální volba pro urychlení konvergence je \\[viz např reference v Numerical recipes\\]\n", "$$\n", "\\omega = \\frac{2}{1+\\sqrt{1-\\rho_{\\rm Jac}^2}}\n", "$$\n", "a odpovídající spektrální poloměr potom je \n", "$$\n", "\\rho_{\\rm SOR} = \\left(\\frac{\\rho_{\\rm Jac}}{1+\\sqrt{1-\\rho_{\\rm Jac}^2}}\\right)^2.\n", "$$\n", "Konkrétně pro Poissonovu rovnici potom máme\n", "$$\n", "\\omega\\approx\\frac{2}{1+\\pi/M},\\qquad\\rho_{\\rm SOR}\\approx1-\\frac{2\\pi}{M}.\n", "$$\n", "Potřebný počet iterací lze vypočítat analogicky k Jacobiho metodě a dostaneme\n", "$$\n", "n\\approx\\frac{1}{3}pM.\n", "$$\n", "Superrelaxací tedy dosahujeme řádového urychlení konvergence. Metoda je vhodná a například ve výpočtech časového vývoje, kde obvykle řešení v předchozím kroku je dobrým počátečním odhadem řešení v kroku následujícím." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Rychlé eliptické řešiče (rapid elliptic solvers, RES)\n", "Využívají specifickou strukturu matic pro eliptické PDR.\n", "\n", "### Cyklická redukce\n", "Prvním příkladem je metoda cyklické redukce. Napišme si tři řádky soustavy Poissonovy rovnice:\n", "$$\\begin{align}\n", "u_{i-2}-2u_{i-1}+&u_i &=&F_{i-1}, &\\rm(C1)\\\\\\\n", " u_{i-1}-2&u_i + u_{i+1} &=&F_i, &\\rm(C2)\\\\\\\n", " &u_i - 2u_{i+1} + u_{i+2} &=&F_{i+1} &\\rm(C3)\n", "\\end{align}$$\n", "a sečtěme $\\rm (C1) + 2(C2) + (C3)$, čímž dostaneme\n", "$$\n", "u_{i-2} - 2u_i + u_{i+2} = F_{i-1} + 2F_i + F_{i+1}.\n", "$$\n", "Obdobnou operaci můžeme provést pro všechna sudá $i$, čímž dostaneme soustavu s polovičním počtem rovnic, která bude mít ovšem stejný tvar jako původní Poissonova rovnice. Když z této redukované soustavy vypočteme hodnoty v sudých bodech, tak hodnoty v lichých bodech potom jednoduše vypočteme z odpovídajících rovnic neredukovanné soustavy.\n", "\n", "Tuto redukci můžeme rekurzivně opakovat. Pokud $M=2^m-1$, tak po $m$ redukcích nám zbyde pouze jedna jedna rovnice, kterou triviálně vyřešíme a potom postupně v $m$ iteracích doplníme chybějící hodnoty.\n", "\n", "Algoritmus lze implementovat i ve 2D, kde se pracuje s celými řadky 2D mříže. Celková náročnost výpočtu je $O(N_g\\log(N_g))$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Fourierovské metody\n", "Fourierovské metody jsou založeny na tom, že známe vlastní vektory matice $K$ a zároveň máme efektivní algoritmus pro rozklad libovolného vektoru do této báze (FFT). Konkrétně v případě Poissonovy rovnice s nulovými dirichletovými okrajovými jsme si ukázali, že vlastní vektory mají sinusový průběh,\n", "$$\n", "z_i^k = \\sin\\frac{ik\\pi}{M+1}.\n", "$$\n", "Příslušnou metodou je tedy diskrétní sinová transformace s předpisem\n", "$$\n", "u_i = \\frac{2}{M+1}\\sum_{k=1}^M\\sin\\left(\\frac{ik\\pi}{M+1}\\right)\\hat u_k\n", "= \\frac{2}{M+1}\\sum_{k=1}^Mz_i^k\\hat u_k\\qquad\\rm(eq:DST)\n", "$$\n", "a její inverze\n", "$$\n", "\\hat u_k = \\sum_{i=1}^M\\sin\\left(\\frac{ik\\pi}{M+1}\\right)\\hat u_i\n", "= \\sum_{i=1}^Mz_i^k\\hat u_i.\n", "$$\n", "Rovnici (eq:DST) můžeme také napsat vektorově jako\n", "$$\n", "U = \\frac{2}{M+1}\\sum_{k=1}^MZ^k\\hat u_k.\n", "$$\n", "Pokud si dosadíme tento rozklad do bázových vektorů za $U$ a $F$ v naší rovnici, dostaneme\n", "$$\n", "KU = \\frac{2}{M+1}\\sum_{k=1}^M\\lambda_kZ^k\\hat u_k = \\frac{2}{M+1}\\sum_{k=1}^MZ^k\\hat f_k.\n", "$$\n", "Z ortogonality bázových vektorů $Z^k$ dále vyplývá, že odpovídající koeficienty na levé a pravé straně se musí rovnat (můžeme například levou i pravou stranu vynásobit $Z^l$), tedy\n", "$$\n", "\\lambda_k\\hat u_k = \\hat f_k\n", "$$\n", "a hledané koeficienty $\\hat u_k$ můžeme snadno explicitně vyjádřit jako\n", "$$\n", "\\hat u_k = \\hat f_k\\left(2\\cos\\frac{k\\pi}{M+1} - 2\\right)^{-1}.\n", "$$\n", "Postup tedy spočívá v tom, že vhodnou variantou rychlé Fourierovy transformace vypočteme $\\hat F$, to vydělíme vlastními čísly a $\\lambda_k$ a poté zpětně transformujeme, čímž dostaneme $U$. Náročnost vypočtu je daná náročností FFT a je tedy $O(N_g\\log(N_g))$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### FACR (Fourier analysis and cyclic reduction)\n", "Asymptoticky nejefektivnější známý přímý algoritmus pro řešení Poissonovy rovnice ve 2D lze získat aplikací $l$ kroků cyklické redukce a následným vyřešením metodou fourierovy transformace v algoritmu zvaném FACR(l) [(detaily viz např. Schwarztrauber et al. 1977)](https://www.jstor.org/stable/2029616). Při vhodné volbě parametru $l(N)$ je asymptotická výpočetní náročnost řádu $O(N_g\\log\\log(N_g))$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To je konec sekce o RES, nyní se podíváme na obecnější metody pro řešení soustav rovnic. Přepneme na znační obvyklé v lineární algebře:\n", "$$\n", "KU = F \\qquad\\to\\qquad Ax=b\n", "$$\n", "\n", "## Metoda sdružených gradientů (Conjugate gradients, CG)\n", "Předpokládejme, že matice soustavy $A$ typu $n\\times n$ je positivně definitní ($\\forall x\\in \\mathbb R^n, |x|>0 \\Rightarrow x^TAx >0$) a symetrická. Potom řešení soustavy $Ax=b$ zároveň minimalizuje funkci\n", "$$\n", "f = \\frac{1}{2}x^TAx - x^Tb,\n", "$$\n", "neboť podmínka nulového gradientu funkce v minimu je ekvivalentní řešené rovnici:\n", "$$\n", "0 = \\nabla f = \\frac{1}{2}\\underbrace{x^TA}_{=A^Tx=Ax} + \\frac{1}{2}Ax-b = Ax-b.\n", "$$\n", "Minimum funkce $f$ můžeme obecně hledat iterativně gradientní metodou. Označíme li si reziduum v $k$-té iteraci\n", "$$r_k = b - Ax_k = -\\nabla f(x_k).$$\n", "Tak podle gradientní metody musíme v každé iteraci provést korekci řešení proti směru gradientu:\n", "$$x_{k+1} = x_k-\\gamma\\nabla f(x_k) = x_k + \\gamma r_k.$$\n", "ovšem v případě kvadratické funkce typu $f$ lze řešení hledat přímo v bázi tzv. sdružených vektorů. \n", "##### Vsuvka: řešení v bázi sdružených vektorů\n", "Dva vektory nazýváme sdružené, pokud jsou vzájemně ortogonální ve skalárním součinu definovaném jako\n", "$\\left_A = u^TAv$. Vzhledem k ortogonalitě tak $n$ sdružených vektorů $\\{p_1,\\ldots,p_n\\}$ tvoří bázi $\\mathbb R^n$.\n", "Vyjádříme-li si hledané řešení jako \n", "$$x = \\sum_{i=1}^n\\alpha_ip_i,$$\n", "vynásobíme tento vektor maticí soustavy\n", "$$Ax = \\sum_{i=1}^n\\alpha_iAp_i,$$\n", "a uvedenou rovnici zleva vynásobíme $p_k^T$, dostaneme\n", "$$p_k^TAx = \\sum_{i=1}^n\\alpha_ip_k^TAp_i.$$\n", "Dále můžeme na levé straně nahradit $Ax=b$ a s využitím ortogonality dostáváme výraz\n", "$$p_k^Tb = \\alpha_ip_k^TAp_k,$$\n", "ze kterého můžeme vyjádřit koeficient rozkladu\n", "$$\\alpha_i = \\frac{\\left}{\\left_A}$$\n", "\n", "Metoda sdružených gradientů je založena na tom, že rozklad řešení do báze sdružených vektorů můžeme konstruovat iterativně s velmi rychlou konvergencí.\n", "První sdružený vektor volíme ve směru gradientu a první korekce řešení tak odpovídá v gradientní metodě. V Dalších iteracích opět vypočítáváme gradient, ovšem ortogonalizujeme jej vůči všem předchozím korekcím. Při vhodné zvolené velikosti lze ukázat, že po $n$ iteracích dostaneme přesné řešení. Obvykle však postačuje mnohem menší počet iterací" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "execution": { "iopub.execute_input": "2024-06-14T19:38:43.801840Z", "iopub.status.busy": "2024-06-14T19:38:43.801704Z", "iopub.status.idle": "2024-06-14T19:38:44.096518Z", "shell.execute_reply": "2024-06-14T19:38:44.096110Z", "shell.execute_reply.started": "2024-06-14T19:38:43.801830Z" } }, "outputs": [], "source": [ "%matplotlib inline\n", "from pylab import *\n", "from numpy import *" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Definice problému: Konstantní nábojová hustota s Dirichletovou OP. Síť se 100 body, krok $h=1$, $\\epsilon_0=1$. Maximální chyba iteračních metod $10^{-3}$" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "execution": { "iopub.execute_input": "2024-06-14T19:38:44.097034Z", "iopub.status.busy": "2024-06-14T19:38:44.096916Z", "iopub.status.idle": "2024-06-14T19:38:44.098995Z", "shell.execute_reply": "2024-06-14T19:38:44.098807Z", "shell.execute_reply.started": "2024-06-14T19:38:44.097026Z" } }, "outputs": [], "source": [ "eps0 = 1\n", "dx = 1\n", "imax = 100\n", "M = imax\n", "rho = ones(imax)\n", "b = -rho/eps0\n", "epsilon = 1e-3" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Iterativní řešení Jacobiho methodou" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "execution": { "iopub.execute_input": "2024-06-14T19:38:44.099367Z", "iopub.status.busy": "2024-06-14T19:38:44.099284Z", "iopub.status.idle": "2024-06-14T19:38:44.552789Z", "shell.execute_reply": "2024-06-14T19:38:44.552397Z", "shell.execute_reply.started": "2024-06-14T19:38:44.099359Z" }, "slideshow": { "slide_type": "-" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "17398\n", "17398\n", "17398\n", "17398\n", "17398\n", "17398\n", "75.1 ms ± 149 µs per loop (mean ± std. dev. of 3 runs, 2 loops each)\n" ] } ], "source": [ "%%timeit -n 2 -r 3\n", "U = zeros(imax+2)\n", "iter = 0\n", "while True:\n", " res = 0.5*(U[:-2] + U[2:] - dx**2*b) - U[1:-1]\n", " U[1:-1] += res\n", " iter += 1\n", " if sum(res**2) < epsilon**2: break\n", "print(iter)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Iterativní řešení Gauss Seidelovou methodou. Nelze v pythonu jednoduše vektorizovat - používáme cyklus, což ji zpomaluje." ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "execution": { "iopub.execute_input": "2024-06-14T19:38:44.553614Z", "iopub.status.busy": "2024-06-14T19:38:44.553518Z", "iopub.status.idle": "2024-06-14T19:38:47.090579Z", "shell.execute_reply": "2024-06-14T19:38:47.090297Z", "shell.execute_reply.started": "2024-06-14T19:38:44.553606Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "9416\n", "9416\n", "9416\n", "9416\n", "9416\n", "9416\n", "422 ms ± 3.24 ms per loop (mean ± std. dev. of 3 runs, 2 loops each)\n" ] } ], "source": [ "%%timeit -n 2 -r 3\n", "U = zeros(imax+2)\n", "iter = 0\n", "while True:\n", " err=0\n", " for i in range(1,imax+1):\n", " res = 0.5*(U[i-1] + U[i+1] - dx**2*b[i-1]) - U[i]\n", " U[i] += res\n", " err += res**2\n", " iter += 1\n", " if err < epsilon**2: break\n", "print(iter)" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "execution": { "iopub.execute_input": "2024-06-14T19:38:47.091047Z", "iopub.status.busy": "2024-06-14T19:38:47.090942Z", "iopub.status.idle": "2024-06-14T19:38:47.166778Z", "shell.execute_reply": "2024-06-14T19:38:47.166477Z", "shell.execute_reply.started": "2024-06-14T19:38:47.091039Z" } }, "outputs": [], "source": [ "from numba import njit" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "execution": { "iopub.execute_input": "2024-06-14T19:38:47.167321Z", "iopub.status.busy": "2024-06-14T19:38:47.167165Z", "iopub.status.idle": "2024-06-14T19:38:47.172189Z", "shell.execute_reply": "2024-06-14T19:38:47.172010Z", "shell.execute_reply.started": "2024-06-14T19:38:47.167313Z" } }, "outputs": [], "source": [ "@njit\n", "def solve_GS(b, dx, imax):\n", " U = zeros(imax+2)\n", " iter = 0\n", " while True:\n", " err=0\n", " for i in range(1,imax+1):\n", " res = 0.5*(U[i-1] + U[i+1] - dx**2*b[i-1]) - U[i]\n", " U[i] += res\n", " err += res**2\n", " iter += 1\n", " if err < epsilon**2: break\n" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "execution": { "iopub.execute_input": "2024-06-14T19:38:47.172623Z", "iopub.status.busy": "2024-06-14T19:38:47.172538Z", "iopub.status.idle": "2024-06-14T19:38:47.513310Z", "shell.execute_reply": "2024-06-14T19:38:47.513022Z", "shell.execute_reply.started": "2024-06-14T19:38:47.172616Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "2.29 ms ± 3.97 µs per loop (mean ± std. dev. of 7 runs, 1 loop each)\n" ] } ], "source": [ "%%timeit\n", "solve_GS(b, dx, imax)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Superrelaxace (Successive overrelaxation (SOR)) s faktorem $\\omega$, bez vektorizace" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "execution": { "iopub.execute_input": "2024-06-14T19:38:47.513755Z", "iopub.status.busy": "2024-06-14T19:38:47.513664Z", "iopub.status.idle": "2024-06-14T19:38:47.666875Z", "shell.execute_reply": "2024-06-14T19:38:47.666606Z", "shell.execute_reply.started": "2024-06-14T19:38:47.513747Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "265\n", "265\n", "265\n", "265\n", "265\n", "265\n", "265\n", "265\n", "265\n", "265\n", "265\n", "265\n", "12.6 ms ± 187 µs per loop (mean ± std. dev. of 3 runs, 4 loops each)\n" ] } ], "source": [ "%%timeit -n 4 -r 3\n", "U = zeros(imax+2)\n", "omega = 2/(1+sqrt(1-cos(pi/imax)**2))\n", "\n", "iter = 0\n", "while True:\n", " err=0\n", " for i in range(1,imax+1):\n", " res = 0.5*(U[i-1] + U[i+1] - dx**2*b[i-1]) - U[i]\n", " U[i] += omega*res\n", " err += res**2\n", " iter += 1\n", " if err < epsilon**2: break\n", "print(iter)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Výše uvedené metody lze vektorizovat iterací zvlášť přes liché a sudé elementy v šachovnicovitém vzoru (red-black ordering)" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "execution": { "iopub.execute_input": "2024-06-14T19:38:47.667236Z", "iopub.status.busy": "2024-06-14T19:38:47.667154Z", "iopub.status.idle": "2024-06-14T19:38:49.766440Z", "shell.execute_reply": "2024-06-14T19:38:49.766023Z", "shell.execute_reply.started": "2024-06-14T19:38:47.667229Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "2.58 ms ± 12.7 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n" ] } ], "source": [ "%%timeit\n", "U = zeros(imax+2)\n", "omega = 2/(1+sqrt(1-cos(pi/imax)**2))\n", "iter = 0\n", "while True:\n", " err=0\n", " res = 0.5*(U[:-2:2] + U[2::2] - dx**2*b[::2]) - U[1:-1:2]\n", " U[1:-1:2] += omega*res\n", " err = sum(res**2)\n", " \n", " res = 0.5*(U[1:-2:2] + U[3::2] - dx**2*b[1::2]) - U[2:-1:2]\n", " U[2:-1:2] += omega*res\n", " err += sum(res**2)\n", " \n", " iter += 1\n", " if err < epsilon**2: break" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Fourier method" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "execution": { "iopub.execute_input": "2024-06-14T19:38:49.767471Z", "iopub.status.busy": "2024-06-14T19:38:49.767357Z", "iopub.status.idle": "2024-06-14T19:38:49.775410Z", "shell.execute_reply": "2024-06-14T19:38:49.775222Z", "shell.execute_reply.started": "2024-06-14T19:38:49.767463Z" } }, "outputs": [], "source": [ "from scipy.fftpack import dst, idst" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "execution": { "iopub.execute_input": "2024-06-14T19:38:49.775779Z", "iopub.status.busy": "2024-06-14T19:38:49.775692Z", "iopub.status.idle": "2024-06-14T19:39:03.663173Z", "shell.execute_reply": "2024-06-14T19:39:03.662685Z", "shell.execute_reply.started": "2024-06-14T19:38:49.775771Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "17.1 µs ± 27.4 ns per loop (mean ± std. dev. of 7 runs, 100,000 loops each)\n" ] } ], "source": [ "%%timeit\n", "rho_k = dst(b, type=1) # type 1 assumes odd func around n=-1 and n=imax\n", "k = arange(1, M+1)\n", "lambda_k = 2*cos(k*pi/(M+1))-2\n", "U_k = rho_k/lambda_k\n", "U = idst(U_k, type=1)/(2*(M+1))" ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "execution": { "iopub.execute_input": "2024-06-14T19:39:03.663722Z", "iopub.status.busy": "2024-06-14T19:39:03.663628Z", "iopub.status.idle": "2024-06-14T19:39:03.665652Z", "shell.execute_reply": "2024-06-14T19:39:03.665469Z", "shell.execute_reply.started": "2024-06-14T19:39:03.663714Z" } }, "outputs": [], "source": [ "k = arange(1, M+1)\n", "lambda_k = 2*cos(k*pi/(M+1))-2" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "execution": { "iopub.execute_input": "2024-06-14T19:39:03.666017Z", "iopub.status.busy": "2024-06-14T19:39:03.665934Z" } }, "outputs": [], "source": [ "%%timeit\n", "rho_k = dst(b, type=1) # type 1 assumes odd func around n=-1 and n=imax\n", "U_k = rho_k/lambda_k\n", "U = idst(U_k, type=1)/(2*(M+1))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In 1D, the solution can be obtained using the fast and simple Thomas algorithm. Instead of it we may use a general sparse solver, which is also easily generalized to 2D." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Direct solver" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from scipy.sparse import dia_matrix, eye\n", "from scipy.sparse.linalg import spsolve\n", "imax = 100\n", "def d2matrix(nelem):\n", " elements = ones((3,nelem))\n", " elements[1,:] *= -2\n", " return dia_matrix((elements, [-1,0,1]), shape=(nelem,nelem)).tocsc()\n", "d2x = d2matrix(imax)/dx**2\n", "rho = ones(imax)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The d2matrix is a matrix of second differences in sparse format:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### zero Dirichlet BC" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "b = -rho" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%%timeit\n", "U = spsolve(d2x, b)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "U = spsolve(d2x, b)\n", "plot(U)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Dirichlet BC" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "b = -rho\n", "BC = array([-1000, 1000])\n", "b[[0, -1]] -= BC/dx**2" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "U = spsolve(d2x, b)\n", "plot(U)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Try direct solver with separate factorization" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from scipy.sparse.linalg import factorized\n", "LU = factorized(d2x)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%%timeit\n", "U = LU(b)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We have obtained a solver, which is more than 1000 times faster than our initial attempts" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Thomas algorithm" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#%%timeit\n", "U = zeros(imax+2)\n", "# first and last element are zero in accord with our BC\n", "U[imax] = 1/(imax+1)*(U[0] - np.sum(arange(1,imax+1)*b) + imax/(imax+1)*U[-1])\n", "for i in range(imax-1, 0, -1):\n", " U[i] = 2*U[i+1] - U[i+2] + b[i]" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "U0 = U" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#%%timeit\n", "from scipy.signal import lfilter, lfiltic\n", "U = zeros(imax+2)\n", "# first and last element are zero in accord with our BC\n", "U[imax] = 1/(imax+1)*(U[0] - np.sum(arange(1,imax+1)*b) + imax/(imax+1)*U[-1])\n", "U[:imax+1] = lfilter([0, 1], [1, -2, 1], r_[b[imax:0:-1], 0, 0], zi = U[-2:])[0][::-1]" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#plot(U)\n", "plot(U0)" ] } ], "metadata": { "kernelspec": { "display_name": "sysenv", "language": "python", "name": "sysenv" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.12.3" } }, "nbformat": 4, "nbformat_minor": 4 }