Zápočtový úkol 3

Zápočtový úkol 3#

LU dekompozice a její využití.

Úkol - zápočet 3

Implementujte LU dekompozici a aplikujte ji na řešení soustavy linearních rovnic. Můžete využít zpětné substituce (zpětný chod Gaussovy eliminace) ze cvičení. Dopřednou substituci (\(\mathbb{L} y = b\)) je třeba doimplementovat zvlášť. Správnost ověřte pomocí knihovní funkce scipy.linalg.solve() aplikovanou na náhodnou matici \(\mathbb{A}\) a náhodné vektory \(\vec{b}_1\), \(\vec{b}_2\) a \(\vec{b}_3\).

Pomocí LU dekompozice spočítejte determinant matice \(\mathbb{A}\). Využite vlastnosti determinantu součinu matic a tvaru matic \(\mathbb{L}\) a \(\mathbb{U}\) pro odvození zjednodušeného vzorce. Výsledek opět zkontrolujte pomocí knihovní funkce numpy.linalg.det(). Jaká je složitost (v rámci big-O notace) výpočtu determinantu pomocí LU dekompozice oproti přímému použití vzorečku z definice determinantu?

Použijte řešení soustavy přes LU dekompozici v metodě iterativního zpřesnění. Nejprve spočítejte první odhad řešení a poté ho iterativně zpřesněte podle tohoto algoritmu. Otestujte na maticích o velikosti \(N = 50 - 500\). Měly by stačit pouze 1-2 iterace. Srovnejte chybu před a po zlepšení. Chybu určete jako \(||\mathbb{A}\vec{x} - \vec{b}||\), kde normu můžete spočítat např. pomocí scipy.linalg.norm().

Metoda iterativního zpřesnění funguje i pro velmi nepřesný nebo náhodný počáteční odhad řešení. Otestujte, kolik je potřeba iterací v takovém případě, aby chyba dosáhla podobné úrovně jako v předchozím případě.

import numpy as np
import matplotlib.pyplot as plt
import scipy.linalg as la

LU dekompozice

A = np.random.rand(4,4)
b1 = np.random.rand(4)
b2 = np.random.rand(4)
b3 = np.random.rand(4)
## DOPLŇTE ##
# kontrola reseni
try:
    np.testing.assert_array_almost_equal(x1, la.solve(A, b1), decimal=7)
except AssertionError as E:
    print(E)
else:
    print("The implementation is correct.")

Výpočet determinantu

## DOPLŇTE ##

Iterativní zpřesnění

N = 500

A = np.random.rand(N, N)
b = np.random.rand(N)
## DOPLŇTE ##

Náhodný odhad řešení

x = np.random.rand(N)

## DOPLŇTE ##