Numerické metody lineární algebry
Contents
3. Numerické metody lineární algebry#
Naimportujeme si knihovny potřebné pro následující příklady:
import numpy as np
3.1. Násobení matic#
A = np.array([
[0, -3.7, 0, 0],
[1, 2, 3, 4.5],
[5, 4, 3, 2],
[1, 2, 3, 2],
[1, 3, 1, 5]
])
B = np.array([
[1, 1, 2, 2, 6],
[1, 2, 3, 4, 7],
[5, 4, 3, 2, 0],
[1, 2, 3, 2, -4]
])
pocet_radku = A.shape[0]
pocet_sloupcu = B.shape[1]
C = np.zeros((pocet_radku,pocet_sloupcu))
for i in range(pocet_radku): #zvolim i-ty radek matice A
for j in range(pocet_sloupcu): # a j-ty sloupec matice B
for k in range(A.shape[1]): # projdu zaroven cely i-ty radek A a j-ty sloupec B
C[i,j] = C[i,j] + A[i,k]*B[k,j] # poscitam vynasobene hodnoty a ulozim je do prvku C[i,j]
print(C)
C_kontrola = np.dot(A,B)
print(C_kontrola)
[[ -3.7 -7.4 -11.1 -14.8 -25.9]
[ 22.5 26. 30.5 25. 2. ]
[ 26. 29. 37. 36. 50. ]
[ 20. 21. 23. 20. 12. ]
[ 14. 21. 29. 26. 7. ]]
[[ -3.7 -7.4 -11.1 -14.8 -25.9]
[ 22.5 26. 30.5 25. 2. ]
[ 26. 29. 37. 36. 50. ]
[ 20. 21. 23. 20. 12. ]
[ 14. 21. 29. 26. 7. ]]
3.2. Přímé metody pro řešení lineárních rovnic \(\mathbb{A}\mathbf{x}=\mathbf{b}\)#
Gaussova metoda
Gauss-Jordanova metoda
LU dekompozice
3.2.1. Gaussova metoda#
Přímý běh: matici \(\mathbb{A}\) převedeme na trojúhleníkový tvar.
Zpětný běh: vypočítáme složky vektoru \(\mathbf{x}\).
3.2.1.1. Řešení soustav s trojúhelníkovou maticí#
Při převodou matice \(\mathbb{A}\) na trojúhelníkový tvar hraje roli výběr hlavního prvku \(a_{11}\) (pivotu) matice \(\mathbb{A}\) (v důsledku omezené přesnosti čísel v počítači)
Příklad na pivoting
Následuje zpětný běh, ve kterém vypočítáme složky vektoru \(\mathbf{x}\) ve směru klesajícího indexu \(i\): \(x_{i}=\dfrac{b_{i}-\sum_{j=i+1}^{n-1}a_{ij}x_{j}}{a_{ii}}\) kde \(i= n-1, n-2,..., 0\)
#
# zpetny beh
A = np.array([[6,9,21],[0,21,-57],[0,0,78]]) # matice v hornim trojuhelnikovem tvaru
b = np.array([[141,-123,312]]).T # vektor b
x = np.zeros((3,1)) # neznamy vektor
n = x.size # pocet prvku (3)
for i in range(n-1,-1,-1): # hodnoty i: 2, 1 , 0
soucet = 0
for j in range(i+1,n,1):
soucet = soucet+x[j]*A[i,j]
x[i] = (b[i]-soucet)/A[i,i]
print('Nami vypocitany vektor x = ' , x)
x_kontrola = np.linalg.solve(A,b)
print('Kontrola: x = ',x_kontrola)
Nami vypocitany vektor x = [[2.]
[5.]
[4.]]
Kontrola: x = [[2.]
[5.]
[4.]]
3.2.2. Gauss-Jordanova metoda#
Matici \(\mathbb{A}\) převedeme na tvar, kdy jsou na hlavní diagonále samé jedničky
Touto metodou lze získat i inverzní matici \(\mathbb{A}^{-1}\)
3.2.3. LU dekompozice#
Matici \(\mathbb{A}\) lze rozepsat jako součin \(\mathbb{A}=\mathbb{L}\mathbb{U}\)
\(\mathbb{L}\) - matice v dolním trojúhelníkovém tvaru, na hlavní diagonále má jedničky
\(\mathbb{U}\) - matice v horním trojúhelníkovém tvaru, na hlavní diagonále má nenulové prvky
Poznámky k LU rozkladu
LU dekompozice pomocí Croutova algoritmu:
\(u_{ij}=a_{ij}-\sum_{k=1}^{i-1}l_{ik}u_{kj}\) pro \(i\le j\)
\(l_{ij}=\left(a_{ij}-\sum_{k=1}^{j-1}l_{ik}u_{kj}\right)/u_{jj}\) pro \(i>j\)
#
# matice A
A = np.array([
[4,3],
[6,3],
])
n = A.shape[0] # pocet radku matice A
U = np.zeros((n, n)) # nulova matice, ze ktere postupne vytvorime horni trojuhlenikovou matici
L = np.eye(n) # jednotkova matice, ze ktere postupne vytvorime dolni trojuhlenikovou matici
for i in range(n):
U[i, i:] = A[i, i:] - np.dot(L[i,:(i-2)],U[:(i-2),i:])
L[(i+1):,i] = (A[(i+1):,i] - np.dot(L[(i+1):,:(i-2)],U[:(i-2),i])) / U[i, i]
print('Matice L = ',L)
print('Matice U = ',U)
Matice L = [[1. 0. ]
[1.5 1. ]]
Matice U = [[ 4. 3. ]
[ 0. -1.5]]
3.3. Speciální typy matic#
3.3.1. Tridiagonální matice#
Matice s kladnými prvky na hlavní a první pod- a naddiagonále.
Ukážeme si implementaci algoritmu řešení soustavy lineárních rovnic s tridiagonální maticí.
#
A = np.array([
[0.0, 3.0, 8.0, 15.0],
[1.0, 14.0, 9.0, 0.3],
[0.8, 1.5, 7.0, 2.0],
[3.0, 6.0, 0.0, 11.0]
])
n = A.shape[0] # pocet radku vstupni matice
c = A[:,0] # pod diagonalou
a = A[:,1] # diagonala
b = A[:,2] # nad diagonalou
f = A[:,3] # prava strana
x = np.zeros((n+1,1)) # vektor reseni
# pomocne vektory
rho = np.zeros((n+1,1))
mu = np.zeros((n+1,1))
for i in range(n):# primy beh
mu[i+1] = -b[i] / ( c[i]*mu[i] + a[i] )
rho[i+1] = ( f[i] - c[i]*rho[i] ) / ( c[i]*mu[i] + a[i] )
for i in range(n-1,-1,-1): # zpetny beh
x[i] = mu[i+1]*x[i+1] + rho[i+1]
print(x)
[[14.54464286]
[-3.57924107]
[ 3.98497024]
[-0.15915179]
[ 0. ]]