Programování pro fyziky 2023/24 – Numerické recepty č. 1

Při výkladu si budeme často pomáhat odkazy do knih Numerical Recipes in C++, Numerical Recipes in C a Numerical Recipes in Fortran (krátce: NR) a také Numerical Computing with MATLAB, které jsou na webu volně k dispozici. Výklad budou provázet ukázky skriptů pro Python s balíčky NumPy a SciPy a také pro Matlab a GNU Octave. V těchto prostředích je k demonstraci numerických metod krátká cesta.

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:

% Octave: operátor \ (backslash)
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:

% Octave: soustava s obecnou maticí
n=10000; A=rand(n)+n*eye(n); b=ones(n,1);
tic; x=A\b; toc; tic; b=A*x; toc    % wallclock time: 7.4 s, 0.04 s

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.

V pythonovském modulu NumPy jsou inicializační funkce pro matice podobné jako v Octavu. Algebraické řešiče solve jsou umístěny v submodulech numpy.linalg a scipy.linalg.

# Python: soustava s obecnou maticí
import numpy as np
from time import time
A=[[0,1],[1,0]]; b=[1,0]; x=np.linalg.solve(A,b); print(x)

n=10000; A=np.random.rand(n,n)+n*np.eye(n); b=np.ones(n)
tic=time(); x=np.linalg.solve(A,b); toc=time(); dt1=toc-tic
tic=time(); b=A@x; toc=time(); dt2=toc-tic
print(dt1,dt2)                      # wallclock time: 6.6 s, 0.0 s

Octave i Python volají algebraické řešiče pro obecné plné matice a další algebraické funkce z numerických knihoven BLAS a LAPACK v některé z osvědčených implementací, jako je OpenBLAS nebo Intel MKL. Jak vidíme, tyto implementace vracejí pro matice o 10000 řádcích řešení v několika sekundách. Algoritmus má kubickou časovou složitost, \(O(n^3)\), a tedy při desetinásobném počtu řádků by výpočetní čas vzrostl 1000krát.

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í s lineární časovou složitostí, \(O(n)\), a soustavy s trojúhelníkovými maticemi s kvadratickou časovou 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 desetitisíců prakticky neměřitelnou dobu běhu. Na rozdíl od Octavu, solver z NumPy diagonální matici nepozná; řešení však můžeme rychle a snadno získat vektorovým zápisem algoritmu, x=numpy.diag(D)/b.

% Octave: soustava s diagonální maticí
n=10000; D=diag(rand(1,n)); b=ones(n,1); tic; x=D\b; toc          % 0.00 s
# Python: soustava s diagonální maticí
import numpy as np
from time import time
n=10000; D=np.diag(np.random.rand(n)); b=np.ones(n)
tic=time(); x=np.linalg.solve(D,b); toc=time(); print(toc-tic)    # 5.8 s
tic=time(); x=b/np.diag(D); toc=time(); print(toc-tic)            # 0.0 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 nazývaným dopředná substituce – prvky vektoru řešení \(x_i\) se postupně získávají řešením jednotlivých 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\]

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\]

V Octavu si obohatíme repertoár o funkci tril pro lower triangular část matice A a triu pro upper triangular a posoudíme projev kvadratické časové složitosti pro n opět pro desítku tisíc. Zkusíme 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. Python s obecnými NumPy a SciPy solvery solve běží s trojúhelníkovými maticemi pomalu, na rozdíl od solveru scipy.linalg.solve_triangular. Můžeme zde oba algoritmy zapsat vlastnoručně, dokonce s obdobnou dobou běhu.

% Octave: soustavy s trojúhelníkovými maticemi
n=10000; L=tril(rand(n))+n*eye(n); b=ones(n,1); tic; x=L\b; toc;  % 1.0 s
n=10000; U=triu(rand(n))+n*eye(n); b=ones(n,1); tic; x=U\b; toc;  % 1.0 s
opts.UT=true; tic; x=linsolve(U,b,opts); toc;                     % 0.2 s
# Python: soustavy s trojúhelníkovými maticemi
import numpy as np
import scipy
from time import time

def MySolveTriL(L,b):
  x=np.empty(b.shape)
  x[0]=b[0]/L[0,0]
  for i in range(1,b.size):
    x[i]=(b[i]-np.dot(L[i,:i],x[:i]))/L[i,i]
  return x

def MySolveTriU(U,b):
  n=b.size
  x=np.empty(b.shape)
  x[n-1]=b[n-1]/U[n-1,n-1]
  for i in range(n-2,-1,-1):
    x[i]=(b[i]-np.dot(U[i,i+1:],x[i+1:]))/U[i,i]
  return x

n=10000; L=np.tril(np.random.rand(n,n))+n*np.eye(n); b=np.ones(n)
tic=time(); x=np.linalg.solve(L,b); toc=time(); print(toc-tic)     # 5.8 s
tic=time(); x=scipy.linalg.solve(L,b); toc=time(); print(toc-tic)  # 8.0 s
tic=time(); x=scipy.linalg.solve_triangular(L,b,lower=True); toc=time(); print(toc-tic) # 0.1 s
tic=time(); x=MySolveTriL(L,b); toc=time(); print(toc-tic)         # 0.05 s
n=10000; U=np.triu(np.random.rand(n,n))+n*np.eye(n); b=np.ones(n)
tic=time(); x=scipy.linalg.solve_triangular(U,b); toc=time(); print(toc-tic) # 0.1 s
tic=time(); x=MySolveTriU(U,b); toc=time(); print(toc-tic)         # 0.05 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 z protějších trojúhelníků. V ukázkách vidíme, že operátor \ ani funkce numpy.linalg.solve třídiagonální matici nerozeznaly a strávily s ní \(O(n^3)\) času jako s plnou maticí. V takové situaci má smysl aktivovat interní řešič pro řídké matice (sparse matrix), v Octavu převodem matice soustavy voláním sparse a užitím operátoru \, v Pythonu konverzní funkcí csr_array (pro jeden ze standardních formátů pro ukládání řídkých matic) a řešičem linalg.spsolve ze submodulu scipy.sparse.

% Octave: soustava s třídiagonální maticí
n=10000; A=triu(tril(rand(n),1),-1)+eye(n); b=ones(n,1); tic; x=A\b; toc;        % 7.0 s
As=sparse(A); tic; xs=As\b; toc; max(abs(x-xs))                                  % 0.0 s
# Python: soustava s třídiagonální maticí
import numpy as np
import scipy
from time import time
n=10000; A=np.triu(np.tril(np.random.rand(n,n),1),-1)+np.eye(n); b=np.ones(n)
tic=time(); x=np.linalg.solve(A,b); toc=time(); print(toc-tic)               # 5.7 s
As=scipy.sparse.csr_array(A)
tic=time(); xs=scipy.sparse.linalg.spsolve(As,b); toc=time(); print(toc-tic) # 0.0 s
print(max(abs(x-xs)))

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 soustavy s maticemi \(\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) Řešení soustav algebraických rovnic pomocí rozkladu na součin neboli faktorizace matice soustavy je konkurenční cestou k eliminačním metodám. Obě cesty jsou algebraicky ekvivalentní, v numerických knihovnách se však ujaly spíše faktorizace, jako např. LU rozklad pro obecné plné matice nebo Choleského rozklad pro symetrické pozitivně definitní matice. Algoritmus má dvě fáze – náročnější faktorizaci a jednodušší řešení soustav se získanými faktory.

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) zajišťuje s časovou složitostí \(O(n^3)\), že seřadíme-li rovnice podle prvků na pravé straně v pořadí \(a_{11},\ a_{21},\dots,\ a_{n1},\ a_{12},\ a_{22},\dots,\ a_{nn}\),

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

v každé rovnici 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é. Ve SciPy zavoláme funkci scipy.linalg.lu a dvakrát scipy.linalg.solve_triangular.

% Octave: LU rozklad
n=10000; A=rand(n)+n*eye(n); b=ones(n,1);
tic; [L,U]=lu(A); y=L\b; x=U\y; toc        % 11.6 s
tic; max(abs(x-A\b)); toc                  % 8.0 s
# Python: LU rozklad
import numpy as np
import scipy
from time import time
n=10000; A=np.random.rand(n,n)+n*np.eye(n); b=np.ones(n)
tic=time(); L,U=scipy.linalg.lu(A,permute_l=True)
y=scipy.linalg.solve_triangular(L,b,lower=True)
x=scipy.linalg.solve_triangular(U,y); toc=time(); print(toc-tic)                     # 7.6 s
tic=time(); print(max(abs(x-scipy.linalg.solve(A,b)))); toc=time(); print(toc-tic) # 9.5 s

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:

% Octave: determinant pomocí LU rozkladu a voláním funkce 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^4\)), 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. Populární a v knihovnách běžně obsažené jsou varianty metody sdružených gradientů (CG), za nejrychlejší iterační metodu se považuje multigridová metoda (MG). 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

Ukázka metody sdružených gradientů (Matlab: pcg, Octave: pcg, SciPy: cg) při řešení diskretizované 2D Poissonovy rovnice, zahrnující 2D Laplaceův operátor. Příslušná Poissonova matice je v Octavu součástí galerie matic, v Pythonu si pro ni lze sáhnout do balíčku PyAMG s implementací multigridové metody. Pro zrychlení konvergence naší ukázky matici upravíme přičtením jedniček na diagonálu, tedy zajištěním diagonální dominance. Na strukturu matice se lze podívat v Octavu funkcí spy a podobně i v Pythonu: matplotlib.pyplot.spy.

% Octave: metoda předpodmíněných sdružených gradientů
format compact; format long
tic
npoisson=2000;
n=npoisson^2;
A=gallery('poisson',npoisson);
A=A+diag(sparse(ones(n,1)));
x=(1:n)';
b=A*x;
tol=1e-12;
maxit=n;

[x,flag,relres,iter,resvec]=pcg(A,b,tol,maxit); % bez předpodmínění
% L=ichol(A);                                   % s předpodmíněním
% [x,flag,relres,iter,resvec]=pcg(A,b,tol,maxit,L,L');

normb=norm(b);
for it=0:10:iter
  disp(sprintf('%5i%10.3e',it,resvec(it+1)/normb))  % norm(b-A*x)
end
disp(sprintf('%5i%10.3e',iter,resvec(iter+1)/normb))  % norm(b-A*x)
x(1:5)'
disp(sprintf('%s %i %i %10.3e','flag,iter,relres =',flag,iter,relres))
toc
# Python: metoda předpodmíněných sdružených gradientů
import numpy as np
import scipy
import pyamg
from time import time
tic=time()
npoisson=2000
n=npoisson**2
A=pyamg.gallery.poisson((npoisson,npoisson),format='csr')
A.setdiag(A.diagonal()+np.ones(n))
x=np.arange(1,n+1)
b=A*x
tol=1e-12
maxit=n

(x,info)=scipy.sparse.linalg.cg(A,b,tol=tol)  # bez předpodmínění

normb=np.linalg.norm(b)
err=np.linalg.norm(b-A*x)/normb
print('err =',err)
print('x =',*[x[i] for i in range(5)])
toc=time()
print('dt =',toc-tic)