Ukázka Gaussovy eliminace

Tomáš Kalvoda, KAM FIT ČVUT, 2019

V souboru gauss.sage jsme implementovali základní verzi Gaussovy eliminace, tak jak je popsaná v textu BI-LIN (vlastní štěstí s psaním kódu můžete zkusit se souborem gauss_fixme.sage).

Sage umí samozřejmě také "Gaussovsky eliminovat". Slouží k tomu metoda echelon_form na objektu typu matice. Tato metoda na matici provede Gaussovu eliminaci malinko důkladněji:

dojde k normalizaci pivotů na $1$ (tj. vydělením řádku vhodným číslem tak, aby na místě pivota byla $1$ )
provede se i zpětný běh Gaussovy eliminace, tj. vyeliminují se i prvky nad pivoty zprava doleva.

# načteme kód
load('gauss.sage')

Příklad 1.8 z prvního cvičení

Gaussova eliminace pro soustavu z Příkladu 1.8 s následující rozšířenou maticí s prvky z $GF(2)$ :

m = matrix(GF(2),[
    [1, 0, 1, 0, 0,  1],
    [1, 1, 0, 1, 1,  1],
    [0, 0, 1, 1, 0,  1],
    [0, 1, 0, 1, 0,  1],
    [0, 1, 0, 0, 1,  1]
])
show(m)

\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrrrrr} 1 & 0 & 1 & 0 & 0 & 1 \\ 1 & 1 & 0 & 1 & 1 & 1 \\ 0 & 0 & 1 & 1 & 0 & 1 \\ 0 & 1 & 0 & 1 & 0 & 1 \\ 0 & 1 & 0 & 0 & 1 & 1 \end{array}\right)

Necháme si vypsat jednotlivé kroky (druhý parametr True). Všimněte si, že výpočty se provádějí v $GF(2)$ .

gauss(m, True)

[1 0 1 0 0 1]
[1 1 0 1 1 1]
[0 0 1 1 0 1]
[0 1 0 1 0 1]
[0 1 0 0 1 1]
====
[1 0 1 0 0 1]
[0 1 1 1 1 0]
[0 0 1 1 0 1]
[0 1 0 1 0 1]
[0 1 0 0 1 1]
====
[1 0 1 0 0 1]
[0 1 1 1 1 0]
[0 0 1 1 0 1]
[0 0 1 0 1 1]
[0 0 1 1 0 1]
====
[1 0 1 0 0 1]
[0 1 1 1 1 0]
[0 0 1 1 0 1]
[0 0 0 1 1 0]
[0 0 0 0 0 0]
====
[1 0 1 0 0 1]
[0 1 1 1 1 0]
[0 0 1 1 0 1]
[0 0 0 1 1 0]
[0 0 0 0 0 0]
====
[1 0 1 0 0 1]
[0 1 1 1 1 0]
[0 0 1 1 0 1]
[0 0 0 1 1 0]
[0 0 0 0 0 0]
====

[1 0 1 0 0 1]
[0 1 1 1 1 0]
[0 0 1 1 0 1]
[0 0 0 1 1 0]
[0 0 0 0 0 0]

Porovnejme tento výsledek s výsledkem metody echelon_form.

m.echelon_form()

[1 0 0 0 1 0]
[0 1 0 0 1 1]
[0 0 1 0 1 1]
[0 0 0 1 1 0]
[0 0 0 0 0 0]

Cvičení 2.3 c) z druhého cvičení

Matice v zadání má prvky z $\mathbb{R}$ , ale i v $\mathbb{Q}$ , tj. při výpočtu si vystačíme s tělesem racionálních čísel.

m = matrix(QQ, [
    [-1,  1,  1,  1,  1],
    [ 1, -1,  1,  1,  1],
    [ 1,  1, -1,  1,  1],
    [ 1,  1,  1, -1,  1],
    [ 1,  1,  1,  1, -1],
])
show(m)

\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrrrr} -1 & 1 & 1 & 1 & 1 \\ 1 & -1 & 1 & 1 & 1 \\ 1 & 1 & -1 & 1 & 1 \\ 1 & 1 & 1 & -1 & 1 \\ 1 & 1 & 1 & 1 & -1 \end{array}\right)

gauss(m, True)

[-1  1  1  1  1]
[ 1 -1  1  1  1]
[ 1  1 -1  1  1]
[ 1  1  1 -1  1]
[ 1  1  1  1 -1]
====
[-1  1  1  1  1]
[ 0  0  2  2  2]
[ 0  2  0  2  2]
[ 0  2  2  0  2]
[ 0  2  2  2  0]
====
[-1  1  1  1  1]
[ 0  2  0  2  2]
[ 0  0  2  2  2]
[ 0  0  2 -2  0]
[ 0  0  2  0 -2]
====
[-1  1  1  1  1]
[ 0  2  0  2  2]
[ 0  0  2  2  2]
[ 0  0  0 -4 -2]
[ 0  0  0 -2 -4]
====
[-1  1  1  1  1]
[ 0  2  0  2  2]
[ 0  0  2  2  2]
[ 0  0  0 -4 -2]
[ 0  0  0  0 -3]
====

[-1  1  1  1  1]
[ 0  2  0  2  2]
[ 0  0  2  2  2]
[ 0  0  0 -4 -2]
[ 0  0  0  0 -3]

m.echelon_form()

[1 0 0 0 0]
[0 1 0 0 0]
[0 0 1 0 0]
[0 0 0 1 0]
[0 0 0 0 1]

Cvičení 3.9 f) z třetího cvičení

Eliminace s komplexními čísly.

m = matrix(SR, [
    [I,     2,     0],
    [1,   2+I, 3-3*I],
    [1+I,   0,     7]
])
show(m)

\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrr} i & 2 & 0 \\ 1 & i + 2 & -3 i + 3 \\ i + 1 & 0 & 7 \end{array}\right)

gauss(m, True)

[       I        2        0]
[       1    I + 2 -3*I + 3]
[   I + 1        0        7]
====
[       I        2        0]
[       0  3*I + 2 -3*I + 3]
[       0  2*I - 2        7]
====
[               I                2                0]
[               0          3*I + 2         -3*I + 3]
[               0                0 -24/13*I + 55/13]
====

[               I                2                0]
[               0          3*I + 2         -3*I + 3]
[               0                0 -24/13*I + 55/13]

m.echelon_form()

[1 0 0]
[0 1 0]
[0 0 1]

Poznámky ke Gaussově eliminaci

Algoritmus popsaný v BI-LIN je opravdu jen základní verze. V praxi při počítání s čísly s pohyblivou desetinnou čárkou je nutné snažit se potlačit chyby vznikající při zaokrouhlování. Například výběr pivota ve sloupci (prvku, vůči kterému se bude eliminovat) popsaný výše není nejvhodnější. Typicky se volí co největší možný prvek.
Gaussova eliminace je numericky nestabilní. Pro některé matice může vést při výpočtu s čísly s pohyblivou desetinnou čárkou k velmi špatným výsledkům. Viz příklad níže.

eps = 1/10000000
A = matrix(QQ, [
    [1+eps, 1,     1],
    [1,     1-eps, 1]
])
B = matrix(RR, [
    [1+eps, 1,     1],
    [1,     1-eps, 1]
])

Nejprve počítejme s první maticí.

Aout = gauss(A)
show(Aout)

\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrr} \frac{10000001}{10000000} & 1 & 1 \\ 0 & -\frac{1}{100000010000000} & \frac{1}{10000001} \end{array}\right)

Čili řešení v tomto případě existuje právě jedno a jeho složky jsou:

y  = Aout[1,2] / Aout[1,1]
x  = (Aout[0,2] - Aout[0,1] * y) / Aout[0,0]
sA = (x,y)
show(sA)

\newcommand{\Bold}[1]{\mathbf{#1}}\left(10000000, -10000000\right)

Nyní stejný proces s druhou maticí.

Bout = gauss(B)
show(Bout)

\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrr} 1.00000010000000 & 1.00000000000000 & 1.00000000000000 \\ 0.000000000000000 & -9.88098491916389 \times 10^{-15} & 9.99999900663795 \times 10^{-8} \end{array}\right)

y  = Bout[1,2] / Bout[1,1]
x  = (Bout[0,2] - Bout[0,1] * y) / Bout[0,0]
sB = (x,y)
show(sB)

\newcommand{\Bold}[1]{\mathbf{#1}}\left(1.01204475834609 \times 10^{7}, -1.01204475955056 \times 10^{7}\right)

Absolutní hodnota rozdílů složek skutečného řešení (první postup) a "přibližného řešení" (druhý postup) je obrovská:

print(abs(sA[0] - sB[0]))
print(abs(sA[1] - sB[1]))

120447.583460858
120447.595505618