{ "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": 39, "metadata": {}, "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": 40, "metadata": {}, "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": 41, "metadata": { "slideshow": { "slide_type": "-" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "17398\n", "17398\n", "17398\n", "17398\n", "17398\n", "17398\n", "166 ms ± 6.95 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", " 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": 42, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "9416\n", "9416\n", "9416\n", "9416\n", "9416\n", "9416\n", "1.27 s ± 11.2 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": "markdown", "metadata": {}, "source": [ "Superrelaxace (Successive overrelaxation (SOR)) s faktorem $\\omega$, bez vektorizace" ] }, { "cell_type": "code", "execution_count": 51, "metadata": {}, "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", "37.4 ms ± 560 µ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 iterace zvlášť přes liché a sudé elementy v šachovnicovitém vzoru (red-black ordering)" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "6.95 ms ± 85.6 µ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": 12, "metadata": {}, "outputs": [], "source": [ "from scipy.fftpack import dst, idst" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "38.4 µs ± 476 ns per loop (mean ± std. dev. of 7 runs, 10000 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": null, "metadata": {}, "outputs": [], "source": [ "k = arange(1, M+1)\n", "lambda_k = 2*cos(k*pi/(M+1))-2" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "32.5 µs ± 1.69 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)\n" ] } ], "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": 35, "metadata": {}, "outputs": [], "source": [ "from scipy.sparse import dia_matrix, eye\n", "from scipy.sparse.linalg import spsolve\n", "imax = 100000\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", "b = -ones(imax)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The d2matrix is a matrix of second differences in sparse format:" ] }, { "cell_type": "code", "execution_count": 36, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[-2. 1. 0. 0. 0. 0. 0. 0. 0. 0.]\n", " [ 1. -2. 1. 0. 0. 0. 0. 0. 0. 0.]\n", " [ 0. 1. -2. 1. 0. 0. 0. 0. 0. 0.]\n", " [ 0. 0. 1. -2. 1. 0. 0. 0. 0. 0.]\n", " [ 0. 0. 0. 1. -2. 1. 0. 0. 0. 0.]\n", " [ 0. 0. 0. 0. 1. -2. 1. 0. 0. 0.]\n", " [ 0. 0. 0. 0. 0. 1. -2. 1. 0. 0.]\n", " [ 0. 0. 0. 0. 0. 0. 1. -2. 1. 0.]\n", " [ 0. 0. 0. 0. 0. 0. 0. 1. -2. 1.]\n", " [ 0. 0. 0. 0. 0. 0. 0. 0. 1. -2.]]\n" ] } ], "source": [ "print(d2matrix(10).todense())" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": 37, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " " ] } ], "source": [ "#%%timeit\n", "%prun U = spsolve(d2x, b)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Try direct solver with separate factorization" ] }, { "cell_type": "code", "execution_count": 29, "metadata": {}, "outputs": [], "source": [ "from scipy.sparse.linalg import factorized\n", "LU = factorized(d2x)" ] }, { "cell_type": "code", "execution_count": 30, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " " ] } ], "source": [ "# %%timeit\n", "%prun U = LU(b)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We have obtained a solver, which is more than 1000 faster than our initial attempts" ] }, { "cell_type": "code", "execution_count": 72, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[]" ] }, "execution_count": 72, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYEAAAEACAYAAABVtcpZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl4VOXZx/HvLUIVpSpvWwRXRKiguKKglRJxw2IFqRWo\nL4uopS4Fq68W1BbauuGuKLUVQRRBRVGDRSRSouKG4gZCBFSo0RJcikutsuR+/3gmkmIIMJnMMzPn\n97muXJw5c+bMnaM593l2c3dERCSZtoodgIiIxKMkICKSYEoCIiIJpiQgIpJgSgIiIgmmJCAikmC1\nJgEzG2dmFWY2v4b3LjSzSjNrWm3fcDNbYmZlZnZctf2HmNn81Hs3Z/ZXEBGRdG2qJDAe6LbhTjPb\nDTgWWF5tXzugN9Au9ZkxZmapt/8MnOHurYHWZvatc4qISPbVmgTc/RngXzW8dQNw8Qb7egCT3X2N\nuy8DlgIdzaw50MTd56aOuxvoWaeoRUQkI7a4TcDMegDl7v7GBm+1AMqrvS4Hdqlh//up/SIiEtnW\nW3KwmTUGLiFUBX2zO6MRiYhI1mxREgBaAXsCr6eq+3cF5plZR8IT/m7Vjt2VUAJ4P7Vdff/7NZ3c\nzDSRkYhIGtw9rQfyLaoOcvf57t7M3Vu6e0vCTf5gd68AioE+ZtbIzFoCrYG57r4C+MzMOqYaivsB\nj9TyHfpxZ8SIEdFjyJUfXQtdC12L2n/qYlNdRCcDzwFtzOw9Mzt9w3t2tZv3QuABYCHwOHCOr4/u\nHGAssARY6u4z6hS1iIhkRK3VQe7edxPv77XB6yuBK2s4bh7QPp0ARUSk/mjEcI4qKiqKHULO0LVY\nT9diPV2LzLC61idlkpl5LsUjIpIPzAzPRsOwiIgUFiUBEZEEUxIQEUkwJQERkQRTEhARSTAlARGR\nBFMSEBFJMCUBEZEEUxIQEUkwJQERkQRTEhARSTAlARGRBFMSEBFJsC1dXlIkr1VWQnk5LFkCS5fC\nsmXwz3/CihVQUQGffQZffBF+1q5d/7kGDWD77cNPkybwgx9A8+aw886wxx7QujXsvTfsvjtsrb8q\nySOaSloKljuUlcGzz8Krr4af+fNhhx3CDXvvvaFly/U382bNwntVN/uGDdefa+3a9cnhs89g5cqQ\nOP75z5BIli4NieXDD2HffeGgg8LPEUfAfvuFJCJSX+oylbSSgBSUt9+Gxx+HWbNgzpxwM//Rj6BD\nh3BTPvDAcKOvL198AW+8ERLOvHkhAa1cGZJB165wwgnQti1YWn+uIjVTEpDEqqyEF1+EKVPgscfg\n88+hWzc47jjo3Bl23TV2hKGaac4cePJJmD497OveHU45Bbp0USlB6k5JQBLntdfg7rvDzb9JE/j5\nz6FnTzjgANgqh7s7uMOiRVBcHGIvL4devaB/f+jUSSUESY+SgCTCv/4FEyfC+PHw8cfhxtmnT6iD\nz1dvvw333w8TJoTkNXAgDBgQ2ihENpeSgBS0BQtg9Gh44IFQpz5oUKhfz+Un/i3lDs8/D+PGwUMP\nwYknwq9/DYcdFjsyyQdaY1gKjnto3D3mmFC/v8suoRpl0qSwr5ASAIRqoCOOgLFj4Z13QgN2795w\n+OHw6KOh7UOkPtRaEjCzcUB3YKW7t0/tuxY4EVgNvA2c7u6fpt4bDgwC1gFD3H1mav8hwF3ANsB0\ndx+6ke9TSSDh3EN9+ZVXwqefwrBh8ItfQKNGsSPLvnXr4OGH4aqr4OuvYfjwkBg0DkE2VG/VQWbW\nGfgCuLtaEjgWmOXulWZ2NYC7DzOzdsAk4FBgF+BJoLW7u5nNBc5z97lmNh24xd1n1PB9SgIJ5Q4z\nZ8Kll4ab32WXhYZe9ZxZf20uvxw++gj+9KfQmFxopSFJX71VB7n7M8C/NthX4u5VhdMXgapOeD2A\nye6+xt2XAUuBjmbWHGji7nNTx90N9EwnWClML7wARUUwdGh48p83D372MyWAKmZw/PHw9NNw002h\nZHDooVBSEjsyKQR1fZYYBKR6PtMCKK/2XjmhRLDh/vdT+yXh3nsPTjst9JcfODA0AJ9yip5wN6Yq\nGbz8cqgaOvts+OlP4a23Ykcm+Szt2kUzuxRY7e6TMhgPI0eO/Ga7qKiIoqKiTJ5ecsBXX8GoUXDL\nLXDOOfCXv4SRvbJ5zEKy/OlPQ6+pH/0odJcdORK++93Y0Uk2lJaWUlpampFzbbKLqJntCUyrahNI\n7RsInAUc7e5fpfYNA3D3q1OvZwAjgOXAbHdvm9rfF+ji7r+q4bvUJlDgnnwyPMG2bw833hgmX5O6\nWbkyVKOVlITqol69NOgsabLaRdTMugEXAT2qEkBKMdDHzBqZWUugNTDX3VcAn5lZRzMzoB/wSDrB\nSv5auTJU/Zx5Zrj5T52qBJApP/hBGF9w773w+9+HMQbLl8eOSvJFrUnAzCYDzwE/NLP3zGwQMBrY\nHigxs1fNbAyAuy8EHgAWAo8D51R7rD8HGAssAZbW1DNICteUKbD//tCiBbz5ZrhJSeb9+Mdh4rqq\nCfP++tfQs0ikNhoxLPXmww9Dnf/8+XDXXWFuHMmON98Mje077hgGoKnUVdg0YlhyzowZYTK3PfcM\nT6dKANm1775hGoquXUN30smTY0ckuUolAcmor74KjZRTp4ZJ0Y46KnZEMm9eGHXdsSPceqt6EBUi\nlQQkJyxeHG405eVhqmclgNxwyCHwyiuwzTZhYZ1582JHJLlESUAyYsqU0CB59tlhu2nT2BFJddtt\nFxqKr7oqLLpz++1qNJZA1UFSJ6tXw0UXwbRp4eZ/yCGxI5JNWbw4DDZr314D9QqFqoMkiooKOPro\nMPXxvHlKAPmiTZswX1PDhmH66nffjR2RxKQkIGl5+eXQ6+Soo8J89zvtFDsi2RKNG4cV2s48M6xZ\nMGtW7IgkFlUHyRabPBmGDAlVCb16xY5G6urvfw+9hy65JPx3lfyj5SUlK9zDXPbjxoWFX/bfP3ZE\nkinvvhtGcnftGqb10MI1+UVJQOrd11/DWWeFJR6Li6F589gRSaatWgU//3lYxe2++6BJk9gRyeZS\nw7DUq1Wrwjz2n38OpaVKAIVqxx1h+vQwx1PnzvDBB7EjkmxQEpBaffBBmJhs//3hwQdDf3MpXA0b\nhvEEp54axn0sXhw7IqlvSgKyUW+9FboQ9u0LN9+s5R6Twiw0El92GXTpAi+9FDsiqU9qE5Aavfxy\naCi88koYNCh2NBJLcTGccQZMmgTHHhs7GtkYNQxLRs2ZE7p+3nEH9OgROxqJ7Zln4Gc/C1NSn3RS\n7GikJnVJAuoIJv+lpCSsAHbvvXryk6Bz59BgfOKJ8OWX0KdP7Igkk5QE5BvFxWEE6dSpcOSRsaOR\nXNKhQ1gf+vjjQyJQFWHhUBIQIEz98Mtfhie+Dh1iRyO5aL/9YPbsMF9UZWV4YJD8pyQg3ySAv/1N\nCUBq16ZNmGaia9fwWokg/ykJJJwSgGyp1q2VCAqJkkCCTZ+uBCDpqZ4IGjaEAQNiRyTpUhJIqNmz\nwx/utGlKAJKe1q1Db7KjjgojyU85JXZEkg4lgQR6/nno3TusBNapU+xoJJ/tsw/MmAHHHQfbbgvd\nu8eOSLZUrdNGmNk4M6sws/nV9jU1sxIzW2xmM81sx2rvDTezJWZWZmbHVdt/iJnNT713c/38KrI5\nXnstDACbMAGKimJHI4XggANC9+LTTw8lTMkvm5o7aDzQbYN9w4ASd28DzEq9xszaAb2BdqnPjDGz\nqhFsfwbOcPfWQGsz2/CckgVvvx2e1MaMgRNOiB2NFJKOHeGBB0IJ85VXYkcjW6LWJODuzwD/2mD3\nScCE1PYEoGdquwcw2d3XuPsyYCnQ0cyaA03cfW7quLurfUayZMWKUGT/3e9Udyv1o6gIbr89jCxe\nujR2NLK50mkTaObuFantCqBZarsF8EK148qBXYA1qe0q76f2S5Z8+ml48h84EH71q9jRSCHr1Qs+\n/jg8cDz7rNaeyAd1ahh2dzezjM74NnLkyG+2i4qKKFLFdZ18/TWcfHKYG/6yy2JHI0lw1lmwcmV4\n8Hj6afjud2NHVHhKS0spLS3NyLk2OYuome0JTHP39qnXZUCRu69IVfXMdvd9zGwYgLtfnTpuBjAC\nWJ46pm1qf1+gi7t/65lUs4hmljv06xfmepkyResBSPa4wznnwDvvwGOPhbEEUn+yvbxkMVA1NGQA\n8Ei1/X3MrJGZtQRaA3PdfQXwmZl1TDUU96v2GalHv/tdaAyeOFEJQLLLDEaPDjf/X/0qJAXJTZvq\nIjoZeA74oZm9Z2anA1cDx5rZYqBr6jXuvhB4AFgIPA6cU+2x/hxgLLAEWOruM+rjl5H17rgD7r8/\ndN1r3Dh2NJJEW28dFqx//XW4/PLY0cjGaFGZAlRSEqqBnnkmjOoUiWnFCjj88JAITjstdjSFSSuL\nyTfKysLC8A8+GP4VyQULFoR5hh59NCQEyaxstwlIjvr449BHe9QoJQDJLfvtB3fdFZapXL48djRS\nnUoCBWL16tA3+7DD4JprYkcjUrObboI774TnnoMmTWJHUzhUHSQMHhzqXqdOVU8gyV3uobfQihXw\n8MOwleoiMkLVQQn3l7/AnDnqCiq5r6rr6CefwB/+EDsaAZUE8t6cOaGedc4c9QSS/FFRAYceGqqH\nevWKHU3+U3VQQr33Xpi9cdw46KZ5WSXPvPxymFpi9uzQcCzpU3VQAn31VSgBDB2qBCD5qUMHuPHG\nMLfVqlWxo0kulQTy1ODBoUvolCmhnlUkX/361/CPf6ihuC5UEkiYu+6Cp54K1UBKAJLvrr8+zDo6\nalTsSJJJJYE889prcOyxUFoK++4bOxqRzCgvDw3FEyfC0UfHjib/qCSQEKtWhXaA0aOVAKSw7Lor\n3HtvmFuovHzTx0vmqCSQJ9xDAmjRAm69NXY0IvXjyivh8cdDj6Gt67TkVbKoJJAAt94a5ly5/vrY\nkYjUn2HDwtTnv/997EiSQyWBPPDyy/CTn8Dzz0OrVrGjEalfK1fCwQfD2LHq/ry5VBIoYKtWwamn\nwpgxSgCSDD/4AUyaBAMHqn0gG1QSyGHu0Lt3+KNQO4AkzRVXhAWSZs3SnFibopJAgRo3LiwSc911\nsSMRyb5hw8I4mKuuih1JYVNJIEeVlUHnzmFQWLt2saMRieP99+GQQ8IU6UccETua3KWSQIH5+mvo\n2zesyaoEIEm2yy5hqvTTTtP8QvVFJYEc9JvfhLlUHnxQ00KIAJx3Hnz0EUyerL+JmqgkUEBKSsLN\n/4479D+7SJVrr4X588OoYskslQRyyCefwAEHwPjxcMwxsaMRyS1V82a9/DLssUfsaHJLlJKAmQ03\nszfNbL6ZTTKz75hZUzMrMbPFZjbTzHbc4PglZlZmZsel+72Fqmrt1VNOUQIQqcmBB8JFF0H//rBu\nXexoCkdaScDM9gTOAg529/ZAA6APMAwocfc2wKzUa8ysHdAbaAd0A8aYmaqiqpk4ERYuVHc4kdpc\neGH4V9OnZE66N+LPgDVAYzPbGmgMfACcBExIHTMB6Jna7gFMdvc17r4MWAoclm7QhWb5crjggpAI\nttkmdjQiuatBA7j77tBG8PrrsaMpDGklAXf/BLge+Afh5r/K3UuAZu5ekTqsAmiW2m4BVB8AXg7s\nklbEBaayEgYNCkngwANjRyOS+/bYIySB/v1h9erY0eS/tCZrNbNWwPnAnsCnwBQz+9/qx7i7m1lt\nrbw1vjdy5MhvtouKiigqKkonxLzx5z/Dv/8d6jpFZPMMGBAGkP3xj2E8TdKUlpZSWlqakXOl1TvI\nzHoDx7r7manX/YBOQFfgKHdfYWbNgdnuvo+ZDQNw96tTx88ARrj7ixucN1G9g5YuhU6d4Nln4Yc/\njB2NSH5ZsSL0pnvssbAqWZLF6B1UBnQys23NzIBjgIXANGBA6pgBwCOp7WKgj5k1MrOWQGtgbprf\nXRDWrYPTT4dLL1UCEEnHzjvDzTeHUsF//hM7mvyV9jgBM7uYcKOvBF4BzgSaAA8AuwPLgFPdfVXq\n+EuAQcBaYKi7P1HDORNTErjxRnj44bBW8FbqJyWSlqqZdqvaCZKqLiUBDRaLoKoa6IUXYO+9Y0cj\nkt8+/BDat4fiYjgsoX0ONW1EHqmshDPOCNVASgAidff978NNN4Xq1a+/jh1N/lESyLLbbw/d2oYM\niR2JSOHo3Rtat05mT6G6UnVQFi1bBh06wDPPQNu2saMRKSwffBDG2jzxBBx0UOxoskvVQXnAHQYP\nDoPClABEMq9FC7jmmjD4cs2a2NHkDyWBLJk4ESoqNChMpD4NGBDaCG68MXYk+UPVQVnw0Uew335h\nUEuHDrGjESls77wTegm9+CK0ahU7muxQF9Ec178/fO97cMMNsSMRSYZrr4WZM8NPEhZnUptADps5\nMzQE//GPsSMRSY7f/AY+/hjuuSd2JLlPJYF69OWXoRrottvghBNiRyOSLPPmwU9+AgsWhHaCQqbq\noBw1fHjoFjp5cuxIRJLpggvCsq133RU7kvqlJJCD3nwTiorC4tg77xw7GpFk+vxz2HffsBBNIc9K\nrzaBHFNZGdYL/sMflABEYmrSJMw0evbZWoBmY5QE6sGECWEOk8GDY0ciIj17hnm6rrsudiS5SdVB\nGfbRR6H4+fjjcPDBsaMREVg/ZcvcubDXXrGjyTy1CeSQM8+E7bYLRVARyR2jRsHTT4dBm4U2dkBJ\nIEe88AL06gWLFsEOO8SORkSqW706LEc5ahScdFLsaDJLDcM5YN06OPfcMIGVEoBI7mnUCEaPhqFD\nwxgeCZQEMuSOO0I10GmnxY5ERDbmmGPCvEJXXx07ktyh6qAM+PDD0Bj85JOw//6xoxGR2pSXh3UH\nCml5V7UJRHbmmbD99mGJOxHJfddcA089BX/7W+xIMkNJIKKXXgqNTGVlagsQyRerV4dS+/XXQ/fu\nsaOpOzUMR1JZGdYKvuIKJQCRfNKoUSi5n3++FqdXEqiDe++FtWth4MDYkYjIlurWDfbZR2N60q4O\nMrMdgbHAvoADpwNLgPuBPYBlwKnuvip1/HBgELAOGOLuM2s4Z95UB33+efgf6KGHoFOn2NGISDqW\nLIHDDw8TPTZvHjua9MWqDroZmO7ubYH9gTJgGFDi7m2AWanXmFk7oDfQDugGjDGzvC6FXHEFHH20\nEoBIPmvdOnTsGDYsdiTxpFUSMLMdgFfdfa8N9pcBXdy9wsx2BkrdfZ9UKaDS3UeljpsBjHT3Fzb4\nfF6UBJYuDTf/N96AFi1iRyMidVFVqp86FTp2jB1NemKUBFoCH5rZeDN7xczuMLPtgGbuXpE6pgJo\nltpuAZRX+3w5sEua3x3dxRfDhRcqAYgUgiZNQsn+/PMhD55BM27rOnzuYOA8d3/JzG4iVfVTxd3d\nzGq7pDW+N3LkyG+2i4qKKMqxlSBKS+GVV2DSpNiRiEim9O8fppS47z7o2zd2NJtWWlpKaWlpRs6V\nbnXQzsDz7t4y9fpIYDiwF3CUu68ws+bA7FR10DAAd786dfwMYIS7v7jBeXO6OmjdOjj0UPjtb6F3\n79jRiEgmPf009OsXxvxsu23saLZM1quD3H0F8J6ZtUntOgZ4E5gGDEjtGwA8ktouBvqYWSMzawm0\nBuam890xTZgQ/uc49dTYkYhIpv34x+Eh74YbYkeSXXXpInoAoYtoI+BtQhfRBsADwO58u4voJYQu\nomuBoe7+RA3nzNmSwOefww9/CI8+Gv5HEZHC88474e97wYL86jKqaSOy4He/C6sT3XNP7EhEpD5d\nfHFYIXDcuNiRbD4lgXr2/vthnpFXX4Xdd48djYjUp08/hTZtoKQkf2YFVhKoZ4MGQbNmcNVVsSMR\nkWwYPTrMMDpjRuxINo+SQD16/XU4/nh46y1NEieSFKtXhzVCbrsNjjsudjSbpllE69HFF8NllykB\niCRJo0ZhLeKLLgpdwwuZkkAtnngC3n0XBg+OHYmIZNvJJ4fRxIXeGUTVQRuxbh0cfDCMGAG9esWO\nRkRieOEF+PnPQ3Vw48axo9k4VQfVg3vvDQvHn3xy7EhEJJZOncKkcqNHx46k/qgkUIOvvgoDwyZO\nhM6dY0cjIjG99RYceWT4t2nT2NHUTCWBDPvzn+GAA5QARCQ8EP7sZ3D11bEjqR8qCWzg00/DQhN/\n/zvst1/UUEQkR/zzn+F+8NprsNtusaP5NpUEMuiaa+DEE5UARGS95s3h7LNDR5FCo5JANR98AO3b\n5262F5F4qqaTmDUr9x4SNWI4Q845J3QDu+66aCGISA674QZ45hl4+OHYkfw3JYEMePttOOyw0APg\ne9+LEoKI5Lj//CeUBh58MLfWI1abQAaMHAlDhigBiMjGbbttmFb+0ktjR5I5KgkQFpA4+mhYujQM\nExcR2Zg1a6BdO7j99nDfyAUqCdTRZZeFdYOVAERkUxo2hD/+ES65BHLoGTptiU8CL74I8+aFRmER\nkc3Ru3eYWaC4OHYkdZf4JHDZZaGOb5ttYkciIvliq63giivC/aOyMnY0dZPoJPDUU2Fh6dNPjx2J\niOSb7t3DJJNTpsSOpG4S2zDsDl26wBlnwIABWflKESkwM2eGXoULFsDWW8eLQw3DaXjySaiogNNO\nix2JiOSrY4+F738fJk2KHUn6ElkScIfDD4ehQ6Fv33r/OhEpYKWloUahrCz0HIohWknAzBqY2atm\nNi31uqmZlZjZYjObaWY7Vjt2uJktMbMyM4u6dPP06fDFF6GFX0SkLoqKYM89YcKE2JGkp67VQUOB\nhUDV4/swoMTd2wCzUq8xs3ZAb6Ad0A0YY2ZRqqLc4fe/hz/8IbTwi4jU1Z/+FH5Wr44dyZZL+zZo\nZrsCPwHGAlXFkJOAqnw4AeiZ2u4BTHb3Ne6+DFgKHJbud9fFtGlh/WAtGykimXLEEdC2LYwfHzuS\nLVeXZ+EbgYuA6r1km7l7RWq7AmiW2m4BlFc7rhzYpQ7fnRb3MEfQiBEqBYhIZo0cCVdemX+lgbQ6\nNZnZicBKd3/VzIpqOsbd3cxqa+Wt8b2RI0d+s11UVERRUY2nT8u0aWFgR48eGTuliAgQFqVv1y6U\nBgYPrt/vKi0tpbS0NCPnSqt3kJldCfQD1gLbAN8FpgKHAkXuvsLMmgOz3X0fMxsG4O5Xpz4/Axjh\n7i9ucN566x3kDoccEkYHqypIROrDCy+EDidLlkCjRtn73qz3DnL3S9x9N3dvCfQB/u7u/YBioGro\n1QDgkdR2MdDHzBqZWUugNTA3ne9Ol0oBIlLfqpcG8kWmasarHt+vBo41s8VA19Rr3H0h8AChJ9Hj\nwDnZnDNabQEiki0jRuRX20AiBosVF4dqoFdfVRIQkfp3wgnQs2f9tw1U0fKStXAPy0b+9rdwyikZ\nPbWISI2eey5MSbN4cXZGEWvuoFo88QR8+SX06hU7EhFJiiOOgL32gokTY0eyaQVdEnCHI4+Ec8+F\nX/wiY6cVEdmk0lI46yxYtKj+ZxhVSWAjZs+GDz/UHEEikn1dusDOO8P998eOpHYFXRI46igYOFDr\nBYhIHDNnwvnnh/UG6rNTikoCNZgzB5YvVzWQiMRz7LHQpAk89FDsSDauYJPA5ZfDsGHx5vcWETEL\n3dMvvzy0UeaigkwC8+aF4peqgUQktu7dw7/Tp8eNY2MKMglcdRX83//Bd74TOxIRSTozuOQSuOKK\n3CwNFFzD8KJFoVX+3Xdhu+0yFJiISB2sWwf77ANjx4b7U6apYbiaUaNgyBAlABHJHQ0ahDbKK6+M\nHcm3FVRJYNmyMF300qWw006Zi0tEpK5Wr4ZWreDhh6FDh8yeWyWBlOuug1/+UglARHJPo0Zw0UWh\nzTKXFExJoKIirPG5aBE0a7bp40VEsu3LL6FlyzClRNu2mTuvSgLALbdAnz5KACKSuxo3hvPOg2uv\njR3JegVREvjsszBj39y54V8RkVz1ySew997wxhuw666ZOWfiSwJ//WsYnq0EICK5rmnTMKfZjTfG\njiTI+5LA11+Hm//f/gYHHlhPgYmIZFB5Oey/f+jJ2LRp3c+X6JLAxInQvr0SgIjkj113hR49YMyY\n2JHkeUlg3Tpo1w5uvz1MGy0iki8WLYKiojC7QePGdTtXYksCjz4KO+wQLqSISD5p2xYOPxzGj48b\nR96WBNzDOp4XXqgF5EUkPz37LPTvHxakb9Ag/fMksiTw7LNh6ciTT44diYhIen70ozC2aerUeDGk\nlQTMbDczm21mb5rZAjMbktrf1MxKzGyxmc00sx2rfWa4mS0xszIzO66ugV93HVxwQd2yp4hIbBdd\nFAaPxaqUSas6yMx2BnZ299fMbHtgHtATOB34yN2vMbPfAju5+zAzawdMAg4FdgGeBNq4e+UG592s\n6qC33oIf/zgzDSoiIjGtWxfaB8aODfe1dGS9OsjdV7j7a6ntL4BFhJv7ScCE1GETCIkBoAcw2d3X\nuPsyYClwWDrfDXD99XD22UoAIpL/GjQIbZuxppKoc5uAme0JHAS8CDRz94rUWxVA1Uw+LYDyah8r\nJySNLVZRAVOmwLnnphWuiEjO6d8/THuzcGH2v3vrunw4VRX0EDDU3T83W18acXc3s9rqdmp8b+TI\nkd9sFxUVUbRB/89bbw0TxX3/++nHLSKSS7bdNjzYXn893Hnnpo8vLS2ltLQ0I9+ddhdRM2sIPAY8\n7u43pfaVAUXuvsLMmgOz3X0fMxsG4O5Xp46bAYxw9xc3OGetbQJffgl77AHPPQetW6cVtohITvro\no3BfKyvb8tmQs94mYOGR/05gYVUCSCkGBqS2BwCPVNvfx8wamVlLoDUwd0u/d8KE0KVKCUBECs33\nvhdqOW67Lbvfm27voCOBp4E3WF+tM5xwY38A2B1YBpzq7qtSn7kEGASsJVQfPVHDeTdaEqisXL9Q\nc7ot6CIiueytt6BzZ1i+PFQRba66lATyZsRwcTH86U+h8cTS+lVFRHLfSSdB9+4wePDmfyYRI4Zv\nuCEMDlMCEJFCdsEFYa2ByspNH5sJeZEE5s2Dd97RHEEiUvi6dIHttoPp07PzfXmRBG64AYYMgYYN\nY0ciIlLXp8zDAAAFwklEQVS/zEJp4IYbsvR9ud4mULUCz7vvhmmjRUQK3Zo1YcXE4mI46KBNH1/Q\nbQK33Qb9+ikBiEhyNGwYBo/dfHP9f1dOlwSqBoc9/zzsvXfEwEREsuzjj8N9b3MGjxVsSWDixLBw\njBKAiCTN//wP9O4dls+tTzlbEnCHffcN1UFaP1hEkmjhQujaNQwe+853Nn5cQZYESkpCvZjWDxaR\npGrXDg48EO67r/6+I2eTwE03wfnna3CYiCTb+eeH+2F9VdrkZBIoKwsDxPr2jR2JiEhcxx0H//kP\nPPVU/Zw/J5PA6NHwy1/CNtvEjkREJK6ttoKhQ+GWW+rn/DnXMLxqldOyJSxYAC1axI5IRCS+L74I\n3eVfeSX8u6GCahgePx6OP14JQESkyvbbw4ABMGZM5s+dcyWBVq2ce+6Bww+PHY2ISO54+23o2BH+\n8Q9o3Pi/3yuoksBOO0GnTrGjEBHJLa1ahcGz996b2fPmXBIYMkTdQkVEajJkSOg4k8kKnJxLAqee\nGjsCEZHcdPTRsHZtZruL5lwSqG1otIhIkpnBr3+d2e6iW2fuVCIiUt/69cvspJo51zsol+IREckH\nBdU7SEREsierScDMuplZmZktMbPfZvO7RUTk27KWBMysAXAr0A1oB/Q1s7bZ+v58U1paGjuEnKFr\nsZ6uxXq6FpmRzZLAYcBSd1/m7muA+4AeWfz+vKL/wdfTtVhP12I9XYvMyGYS2AV4r9rr8tQ+ERGJ\nJJtJQN1+RERyTNa6iJpZJ2Cku3dLvR4OVLr7qGrHKFGIiKQh3S6i2UwCWwNvAUcDHwBzgb7uvigr\nAYiIyLdkbcSwu681s/OAJ4AGwJ1KACIiceXUiGEREcmunBgxnORBZGa2m5nNNrM3zWyBmQ1J7W9q\nZiVmttjMZprZjrFjzRYza2Bmr5rZtNTrRF4LM9vRzB40s0VmttDMOib4WgxP/Y3MN7NJZvadpFwL\nMxtnZhVmNr/avo3+7qlrtSR1Tz1uU+ePngQ0iIw1wG/cfV+gE3Bu6vcfBpS4extgVup1UgwFFrK+\nR1lSr8XNwHR3bwvsD5SRwGthZnsCZwEHu3t7QnVyH5JzLcYT7o/V1fi7m1k7oDfhXtoNGGNmtd7n\noycBEj6IzN1XuPtrqe0vgEWE8RMnARNSh00AesaJMLvMbFfgJ8BYoKq3Q+KuhZntAHR293EQ2tTc\n/VMSeC2AzwgPS41THUwaEzqXJOJauPszwL822L2x370HMNnd17j7MmAp4R67UbmQBDSILCX1xHMQ\n8CLQzN0rUm9VAM0ihZVtNwIXAZXV9iXxWrQEPjSz8Wb2ipndYWbbkcBr4e6fANcD/yDc/Fe5ewkJ\nvBbVbOx3b0G4h1bZ5P00F5KAWqYBM9seeAgY6u6fV38vNb92wV8nMzsRWOnur7K+FPBfknItCD33\nDgbGuPvBwL/ZoLojKdfCzFoB5wN7Em5y25vZ/1Y/JinXoiab8bvXel1yIQm8D+xW7fVu/HcmK3hm\n1pCQAO5x90dSuyvMbOfU+82BlbHiy6IjgJPM7F1gMtDVzO4hmdeiHCh395dSrx8kJIUVCbwWHYDn\n3P1jd18LTAUOJ5nXosrG/iY2vJ/umtq3UbmQBF4GWpvZnmbWiNCoURw5pqwxMwPuBBa6+03V3ioG\nBqS2BwCPbPjZQuPul7j7bu7ektDw93d370cyr8UK4D0za5PadQzwJjCNhF0LQoN4JzPbNvX3cgyh\n40ASr0WVjf1NFAN9zKyRmbUEWhMG5m6cu0f/AU4gjCZeCgyPHU+Wf/cjCfXfrwGvpn66AU2BJ4HF\nwExgx9ixZvm6dAGKU9uJvBbAAcBLwOuEp98dEnwtLiYkwfmEhtCGSbkWhFLxB8BqQvvp6bX97sAl\nqXtpGXD8ps6vwWIiIgmWC9VBIiISiZKAiEiCKQmIiCSYkoCISIIpCYiIJJiSgIhIgikJiIgkmJKA\niEiC/T9DYzjOtYjFfQAAAABJRU5ErkJggg==\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "U = LU(b)\n", "plot(U)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "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.6.7" } }, "nbformat": 4, "nbformat_minor": 1 }