Programování pro fyziky 2021/22 – Numerické recepty č. 1

Náš průvodní text je především v poznámkách k přednášce. Pomáhat si budeme i odkazy do knih Numerical Recipes in C a Numerical Recipes in Fortran (krátce: NR), které jsou (v Second Edition) na webu volně k dispozici. Vidět je ovšem lze jen z Adobe Readeru po instalaci pluginu FileOpen. Interní PDF prohlížeče Firefoxu a Chromu nám nepomohou.

Hádanky

  1. Jak velkou soustavu lineárních algebraických rovnic s plnou maticí vyřeší počítač za sekundu?

  2. Totéž pro soustavu s trojúhelníkovou maticí?

  3. Totéž pro soustavu s třídiagonální maticí?

  4. Totéž pomocí iterační metody?

Soustavy lineárních algebraických rovnic

(NR kap. 2) Úlohou je řešit soustavu

\[\mathbf{A}\cdot\mathbf{x}=\mathbf{b} \quad\mbox{neboli po řádcích}\quad \sum_{j=1}^n a_{ij} x_j = b_i,\ i=1,\dots,n,\]

kde \(\mathbf{A}\) je čtvercová matice soustavy, \(\mathbf{x}\) vektor neznámých a \(\mathbf{b}\) pravá strana. Gaussova eliminační metoda sice bývá algebraickou metodou první volby, v praktickém počítání jsou však upřednostněny (algebraicky ekvivalentní) faktorizační metody. Vedle těchto tzv. přímých metod stojí metody iterační, určené pro velké soustavy s řídkými maticemi.

Octave nabízí pro řešení soustav algebraických rovnic operátor \ („dělení maticí zleva“). Čtvercovou matici A vytvoříme pomocí konstruktoru pole, v kterém se řádky oddělují pomocí středníků, a podobně lze zapsat i sloupcový vektor b:

A=[0 1;1 0]; b=[1;0]; x=A\b
% x = 0
%     1

K větším maticím nám mohou pomoci mnohé předdefinované funkce, jako třeba rand(n) pro čtvercovou matici naplněnou pseudonáhodnými čísly nebo eye(n) pro identickou matici \(I_n\). Po vynásobení skalárem n*eye(n) a sečtení s rand(n) máme v ruce diagonálně dominantní matici a funkce ones(n,1) vytvoří sloupcový vektor naplněný jedničkami. Užitím dvojice funkcí tic a toc si můžeme změřit, kolik času Octave k řešení soustavy s maticí řekněme o několika tisících řádků spotřebuje:

n=5000; A=rand(n)+n*eye(n); b=ones(n,1); tic; x=A\b; toc; tic; b=A*x; toc

Zkouška A*x má vrátit vektor b, jehož několik prvních prvků si pro ověření můžeme vypsat pomocí např. b(1:3)', kde apostrof transponuje sloupcový vektor na řádkový.

Alternativou k operátoru \ je funkce linsolve(A,b,opts), umožňující svým třetím argumentem upřesnit typ matice a tím volbu vhodné metody. Obě cesty volají přímé řešiče z osvědčených numerických knihoven pro plné i řídké matice (LAPACK a UMFPACK).

Snadné případy

Snadno řešitelnými případy, co do složitosti algoritmu i časové složitosti řešení, jsou soustavy s diagonální a třídiagonální maticí (časová složitost \(O(n)\)) a soustavy s trojúhelníkovými maticemi (časová složitost \(O(n^2)\)).

Soustava s diagonální maticí, \(\mathbf{D}\cdot\mathbf{x}=\mathbf{b},\ d_{ij}=0\) pro \(i\ne j\) (a samozřejmě \(d_{ii}\ne0\)) je řešitelná lineárním cyklem \(x_i=b_i/d_{ii},\ i=1,\dots,n\). V Octavu funkce rand(1,n) vygeneruje (řádkový) vektor náhodných čísel, funkce diag vytvoří diagonální matici s vektorovým argumentem na diagonále. Pomocí funkcí tic; toc pak změříme, že lineární časová složitost, \(O(n)\), znamená pro n v řádu tisíců prakticky neměřitelnou dobu běhu.

n=5000; D=diag(rand(1,n)); b=ones(n,1); tic; x=D\b; toc          % 0.00 s

Soustava s dolní (lower) trojúhelníkovou maticí, \(\mathbf{L}\cdot\mathbf{x}=\mathbf{b},\ l_{ij}=0\) pro \(i<j\), se řeší algoritmem dopředná substituce – neznámé \(x_i\) se postupně získávají řešením rovnic směrem od první k poslední:

\[x_1=b_1/l_{11},\ x_i=\left(b_i-\sum_{j=1}^{i-1}l_{ij}x_j\right)/l_{ii},\ i=2,\dots,n\]

V Octavu si obohatíme repertoár o funkci tril pro lower triangular část matice A a posoudíme projev kvadratické časové složitosti pro n opět kolem tisíců:

n=5000; L=tril(rand(n))+n*eye(n); b=ones(n,1); tic; x=L\b; toc;  % 0.31 s

Soustava s horní (upper) trojúhelníkovou maticí, \(\mathbf{U}\cdot\mathbf{x}=\mathbf{b},\ u_{ij}=0\) pro \(i>j\), má svůj algoritmus zpětné substituce, kdy se neznámé \(x_i\) postupně získávají řešením rovnic směrem od poslední k první (NR 2.2.4):

\[x_n=b_n/u_{nn},\ x_i=\left(b_i-\sum_{j=i+1}^{n}u_{ij}x_j\right)/u_{ii},\ i=n-1,\dots,1\]

Octave o upper triangular matici požádáme voláním funkce triu a zkusíme zde i volání funkce linsolve s uvedením argumentu sdělujícího, že matice soustavy je horní trojúhelníková. Toto úsilí Octave odmění kratší dobou běhu:

n=5000; U=triu(rand(n))+n*eye(n); b=ones(n,1); tic; x=U\b; toc;  % 0.31 s
opts.UT=true; tic; x=linsolve(U,b,opts); toc;                    % 0.16 s

Soustava s třídiagonální maticí s nulami všude vyjma prvků \(a_{i-1,i},\ a_{ii}\) a \(a_{i+1,i}\) se řeší sekvencí dvou cyklů – eliminací subdiagonálních prvků a zpětnou substitucí pro matici se superdiagonálou. Třídiagonální matici získáme společnou akcí funkcí triu a tril, zde s druhými argumenty umožňujícími ponechat ze zdrojové matice několik subdiagonál, resp. superdiagonál v protějších trojúhelnících:

n=5000; A=triu(tril(rand(n),1),-1)+eye(n); b=ones(n,1); tic; x=A\b; toc; % 3.79 s
As=sparse(A); tic; xs=As\b; toc; max(abs(x-xs))                  % 0.001 s

V ukázkách vidíme, že operátor \ třídiagonální matici nerozeznal a strávil s ní \(O(n^3)\) času jako s plnou maticí. V takové situaci má smysl aktivovat interní řešič pro řídké matice, zde převodem matice soustavy voláním sparse. Prověřením odchylky absolutní hodnoty rozdílu řešení získaných oběma cestami si potvrdíme, že na výstupu máme totéž při podstatně jiných dobách běhu.

Gaussova eliminační metoda

(NR kap. 2.1–2.2) Nulováním nediagonálních, resp. poddiagonálních prvků matice soustavy pomocí odčítání vhodných násobků jiných řádků (eliminačními transformacemi) se kráčí ke tvaru

\[\mathbf{D}\cdot\mathbf{x}=\mathbf{b'}, \quad\mbox{resp.}\quad \mathbf{U}\cdot\mathbf{x}=\mathbf{b''},\]

kde \(\mathbf{D}\) je diagonální matice a \(\mathbf{U}\) horní trojúhelníková matice. Časová složitost eliminace je \(O(n^3)\), neboť \(n^2-n\), resp. \((n^2-n)/2\) prvků se eliminuje odčítáním \(n\)-prvkových řádků, časová složitost řešení pro \(\mathbf{D}\), resp. \(\mathbf{U}\) je uvedena výše. Pro numerickou stabilitu je nezbytná sloupcová (záměna jen sloupců), řádková (záměna jen řádků) nebo úplná pivotace tak, aby v odčítané rovnici prvek (pivot) s největší absolutní hodnotou (normalizovaných rovnic) ležel na diagonále.

LU rozklad

(NR kap. 2.3) Krok „A“: Nalezne se rozklad matice soustavy na součin dvou faktorů, \(\mathbf{A}=\mathbf{L}\cdot\mathbf{U}\), kde \(\mathbf{L}\) je dolní (lower) trojúhelníková (zde s jedničkami na diagonále, \(l_{ii}=1\)) a \(\mathbf{U}\) horní (upper) trojúhelníková matice. Úloha je vlastně soustavou \(n^2\) nelineárních rovnic pro \(n^2\) neznámých \(l_{ij},\ u_{ij},\)

\[\begin{split}\left(\begin{array}{ccc} 1 & 0 & 0 \\ & \ddots & 0 \\ l_{ik} & & 1 \end{array}\right) \cdot \left(\begin{array}{ccc} u_{11} & & u_{kj} \\ 0 & \ddots & \\ 0 & 0 & u_{nn} \end{array}\right) = \left(\begin{array}{ccc} \\ & a_{ij} & \\ \\ \end{array}\right),\end{split}\]

a je přímočará při vhodné volbě řazení rovnic. Croutův algoritmus „po sloupcích zleva“ (NR obr. 2.3.1) s časovou složitostí \(O(n^3)\) zajišťuje, že v každé rovnici

\[\left( \sum_{k=1}^{\min(i,j)}l_{ik}u_{kj}=a_{ij},\quad i=1,\dots,n \right), \quad j=1,\dots,n\]

se objeví právě jeden dosud nespočtený prvek matic \(\mathbf{L}\) a \(\mathbf{U}\). Hledané prvky tak lze snadno získat v pořadí \(u_{11},\ l_{21},\dots,\ l_{n1},\ u_{12},\ u_{22},\dots,\ l_{n2}\) atd.

Krok „b“: Soustavu \(\mathbf{A}\cdot\mathbf{x}=\left(\mathbf{L}\cdot\mathbf{U}\right)\cdot\mathbf{x}=\mathbf{L}\cdot\left(\mathbf{U}\cdot\mathbf{x}\right)=\mathbf{b}\) pak lze řešit jako dvě soustavy s trojúhelníkovými maticemi,

\[\mathbf{L}\cdot\mathbf{y}=\mathbf{b},\quad\mathbf{U}\cdot\mathbf{x}=\mathbf{y},\]

každou s časovou složitostí \(O(n^2)\). Vektor \(\mathbf{b}\) je potřeba až v kroku „b“.

V Octavu lze LU rozklad získat voláním funkce lu a soustavy s faktory L a U pak ponechat k řešení operátoru \. Doba běhu je blízká přímému volání A\b, výsledné vektory řešení jsou totožné:

n=5000; A=rand(n)+n*eye(n); b=ones(n,1);
tic; [L,U]=lu(A); y=L\b; x=U\y; toc
max(abs(x-A\b))                          % ans = 0

Pro determinant platí \(\det\mathbf{A}=\det\mathbf{L}\cdot\det\mathbf{U}=\prod_{i=1}^n u_{ii}\), neboť determinanty trojúhelníkových matic jsou snadné. V Octavu by bylo třeba zohlednit provádění pivotace, tedy uložit na výstupu funkce lu i permutační matici a zpracovat i ji, anebo volat přímo funkci det:

n=100; A=rand(n)+n*eye(n); b=ones(n,1);
[L,U,P]=lu(A);                           % P.A=L.U
det(A)-prod(diag(U))*det(inv(P))         % ans = 0

Iterační metody

Jsou užívány pro řešení velkých soustav (\(n\gg10^3\)), které ovšem mají řídkou matici (s mnoha nulovými prvky). Soustava \(\mathbf{A}\cdot\mathbf{x}=\mathbf{b}\) se upraví na tvar

\[\mathbf{x}^{(k+1)}=\mathbf{H}\cdot\mathbf{x}^{(k)}+\mathbf{g},\]

kde \(\mathbf{H}\) se nazývá iterační maticí. Efektivní iterační řešiče a předpodmiňující metody jsou součástí řady knihoven, často vyvíjených pro klastry s paralelizací pomocí knihovny MPI. Za nejrychlejší iterační metodu se považuje multigridová (MG) metoda, populární jsou varianty metody sdružených gradientů (CG). Dvě následující metody jsou jako samostatné metody spíše učebnicovými ukázkami, vyskytující se však frekventovaně jako komponenty složitějších metod.

(Jacobiova metoda) Pro rozklad \(\mathbf{A}=\mathbf{L'}+\mathbf{D}+\mathbf{U'}\), kde \(\mathbf{L'}\) a \(\mathbf{U'}\) jsou trojúhelníkové matice s nulami na diagonále, lze v řešené soustavě vše kromě \(\mathbf{D}\cdot\mathbf{x}\) převést na pravou stranu a doplněním iteračních indexů získat iterační vzorec

\[\mathbf{D}\cdot\mathbf{x}^{(k+1)}=-(\mathbf{L'}+\mathbf{U}')\cdot\mathbf{x}^{(k)}+\mathbf{b}.\]

Pokud posloupnost konverguje, limitou je řešení soustavy \(\mathbf{A}\cdot\mathbf{x}=\mathbf{b}\). Poznamenejme, že na pořadí vyčíslování prvků vektoru \(\mathbf{x}^{(k+1)}\) nezáleží a metoda je tak vhodná pro paralelizaci.

(Gaussova-Seidelova metoda) Převedením \(\mathbf{U'}\cdot\mathbf{x}\) na pravou stranu a doplněním iteračních indexů se získá soustava s dolní trojúhelníkovou maticí řešitelná dopřednou substitucí, tedy sériově, od prvního prvku k poslednímu,

\[(\mathbf{L'}+\mathbf{D})\cdot\mathbf{x}^{(k+1)}=-\mathbf{U'}\cdot\mathbf{x}^{(k)}+\mathbf{b}.\]

Obě metody konvergují pro ostře diagonálně dominantní matice, Gaussova-Seidelova metoda i pro symetrické pozitivně definitní matice, zkusit je lze i v jiných případech. Rychlost konvergence obou metod je však malá.

n=5;
A=rand(n)+n*eye(n); b=A*(1:n)';      % volíme vektor řešení ve tvaru (1:n)'
L=tril(A,-1); Dvec=diag(A); U=triu(A,1); D=diag(Dvec);

% Jacobi...
x=zeros(n,1);
% opakovat do konvergence...
x=D\(-(L+U)*x+b); x'                 % 5 platných míst po 11 iteracích

% Gauss-Seidel...
x=zeros(n,1);
% opakovat do konvergence...
x=(L+D)\(-U*x+b); x'                 % 5 platných míst po 5 iteracích