Soustavy lineárních algebraických rovnic

Mnoho algoritmů numerické matematiky, jako např. soustavy diferenciálních rovnic, míří k řešení soustav lineárních algebraických rovnic,

\[\mathbf{A}\cdot\mathbf{x}=\mathbf{b} \ ,\qquad \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. Někdy jsou to soustavy s plnou maticí, jindy s pásovou maticí, často s řídkou maticí, soustava rovnic ani nemusí být regulární. Řešení soustav algebraických rovnic se obvykle neprogramuje, spíše vybírá z nabídky numerických knihoven; příkladem budiž Matlab, přizpůsobující algoritmus zapozdřený v operátoru \ vlastnostem matice soustavy:

% MATLAB algebraic solver
n=5000;
A=rand(n); b=ones(n,1);
tic; x=A\b; toc;       % backslash operator
err=max(abs(b-A*x))    % error estimate

Ukážeme si ve Fortranu.

Studijní text: Numerické recepty kap. 2.1–2.5.

Soustavy s diagonální a trojúhelníkovou maticí

Snadno řešitelnými případy, co do složitosti algoritmu i časové složitosti řešení, jsou soustavy s diagonální a trojúhelníkovou maticí.

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\)) má zjevné řešení \(x_i=b_i/d_{ii},\ i=1,\dots,n\). Jediný cyklus implikuje lineární časovou složitost vzhledem k počtu rovnic, tedy \(O(n)\).

! System with diagonal matrix
! Matlab ... n=5000; D=diag(rand(1,n)); b=ones(n,1); tic; x=D\b; toc; err=max(abs(b-D*x))
real(8),allocatable :: D(:),b(:),x(:)
nmax=5000
allocate (D(nmax),b(nmax))
call random_number(D); b=1.
x=b/D
err=maxval(abs(b-D*x))
print *,'err =',err
end program

Pro soustavu s dolní (lower) trojúhelníkovou maticí, \(\mathbf{L}\cdot\mathbf{x}=\mathbf{b},\ l_{ij}=0\) pro \(i<j\), stačí použít tzv. dopřednou substituci, kdy neznámé prvky vektoru \(x_i\) získáváme jeden po druhém řešením rovnic směrem od první k poslední,

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

Vidíme cyklus s proměnnou \(i\) procházející jednotlivé rovnice a vnořený součtový cyklus přes \(j\), algoritmus bude mít kvadratickou časovou složitost, \(O(n^2)\).

! System with lower triangular matrix
! Matlab ... n=5000; L=tril(rand(n))+n*eye(n); b=ones(n,1); tic; x=L\b; toc; err=max(abs(b-L*x))
real(8),allocatable :: L(:,:),b(:),x(:)
nmax=5000
allocate (L(nmax,nmax),b(nmax),x(nmax))
call random_number(L); forall (i=1:nmax,j=1:nmax,i<j); L(i,j)=0.; endforall; b=1.
forall (i=1:nmax); L(i,i)=L(i,i)+nmax; endforall  ! diagonally dominant matrix
x(1)=b(1)/L(1,1)                                  ! algorithm of forward substitution
do i=2,nmax
  x(i)=(b(i)-dot_product(L(i,1:i-1),x(1:i-1)))/L(i,i)
enddo
err=maxval(abs(b-matmul(L,x)))
print *,err
end program

Ničím už nepřekvapí analogický algoritmus zpětné substituce pro soustavu s horní (upper) trojúhelníkovou maticí, \(\mathbf{U}\cdot\mathbf{x}=\mathbf{b},\ u_{ij}=0\) pro \(i>j\), pro získání neznámých \(x_i\) průchodem rovnic směrem od poslední k první,

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

Soustavy s plnou maticí

Gaussova eliminace

V matematice je klíčovým algoritmem k řešení algebraických soustav Gaussova eliminace. Je i výchozím algoritmem Numerických receptů (kap. 2.1), kde lze pro zopakování látky nalézt jeho popis. Podstatnou zprávou je, že v případě plných matic se jedná o algoritmus s kubickou časovou složitostí, \(O(n^3)\), neboť pro eliminaci každého z kvadratického počtu prvků se používají řádkové operace o lineární složitosti. Je zde i ukázkový recept – použitelný zdrojový kód ve Fortranu (alternativně i v C). Je zde ovšem i varování – hlavní cíle autorů jsou pedagogické, nikoliv výkonostní. Pro vyzkoušení receptu gaussj lze užít následující ukázku:

! Testing Numerical Recipe gaussj
! gfortran -O3 testGaussj1.f90 DGAUSSJ.FOR
program testGaussj
integer,parameter :: wp=8  ! working precision
real(wp),allocatable :: a(:,:),a_orig(:,:),b(:),b_orig(:),x(:)
np=500                     ! physical size of arrays
n=np                       ! number of equations
mp=1                       ! max number of right-hand-side vectors
m=mp                       ! actual number of rhs vectors
allocate (a(np,np),b(np))  ! matrix, rhs vector
call random_number(a)      ! random elements
a_orig=a                   ! saved copy
b=1.                       ! rhs vector
b_orig=b                   ! saved copy
call gaussj(a,n,np,b,m,mp) ! Gauss-Jordan elimination
x=b                        ! solution vector
err=maxval(abs(b_orig-matmul(a_orig,x))) ! error estimate
print *,'err =',err
end program

Na vstupu je matice soustavy a, matice pravých stran b (v ukázce jí je předáván jen první sloupec, tedy jeden vektor pravé strany) a dále – jak je v mnohých starších knihovnách zvykem – velikost aktuálního problému (zde počet rovnic n a počet pravých stran m) a alokační velikost polí (zde velikost np první dimenze pole a a velikost mp druhé dimenze pole b). Na výstupu recept přepisuje pole a inverzní maticí a pole b příslušnými vektory řešení. Jak je zřejmé, voláme double precision neboli real(8) variantu receptu – to je při řešení algebraických soustav běžné, šířící se zaokrouhlovací chyby by mohly z platných číslic real(4) dat ukousnout příliš. Recept jsme ponechali v samostatném souboru, neboť naše ukázky budeme zapisovat ve volném formátu Fortranu 90, zatímco recepty budeme přebírat v pevném formátu starších norem Fortranu, což vyjadřují i přípony souborů. Nebudeme však nechávat recepty jako samostatné procedury. To by nás ochuzovalo o řadu možností a kontrol, které překladače nabízejí, když recepty vnoříme do fortranských modulů. Soubor s receptem tedy upravíme takto:

MODULE mNR                         ! module header
INTEGER,PARAMETER :: NMAX=10000    ! module constant (max expected value of np)
CONTAINS

SUBROUTINE gaussj(a,n,np,b,m,mp)
INTEGER m,mp,n,np                  ! deleted declaration of local constant NMAX
DOUBLE PRECISION a(np,np),b(np,mp)
! PARAMETER (NMAX=50)              ! deleted initialization of local NMAX
INTEGER i,icol,irow,j,k,l,ll,indxc(NMAX),indxr(NMAX),ipiv(NMAX)
! (rest of subroutine)
END SUBROUTINE

END MODULE

Zde jsme provedli některé z bodů následujícího seznamu (NR-in-modules checklist), kterého se budeme držet i v budoucnu:

  • vnoření receptů do modulu,

  • doplnění závěrečného END o typ programové jednotky (END SUBROUTINE, END FUNCTION),

  • nahrazení zastaralého příkazu pause téměř ekvivalentem print*,,

  • přesun vhodných lokálních konstant do záhlaví modulu,

  • přesun obsahu COMMON bloků do záhlaví modulu a odstranění popisů COMMON,

  • odstranění popisů EXTERNAL a popisů typu funkcí pro recepty zahrnuté do modulu,

  • udržení limitu 72 znaků na rozeznávanou délku řádku.

Ukázka použití je téměř stejná jako výše, navíc je jen řádek s připojením modulu, use mNR. Při překladu je vhodné přeložit jako první soubor s modulem a teprve pak soubor modul využívající.

! Testing Numerical Recipe gaussj
! gfortran -O3 mNR.for testGaussj2.f90
program testGaussj
use mNR
integer,parameter :: wp=8  ! working precision
! (as above)

(NEDOKONČENO)

Zdrojové soubory:

  • soustavy s plnou maticí (přednáška č. 2)

  • soustavy s plnou maticí (přednáška č. 3)

  • soustavy s pásovou maticí (přednáška č. 3)