4.4. Linearna algebra¶
Kao vektore i matrice koristit ćemo jednodimenzionalna i dvodimenzionalna NumPy polja te rutine NumPy modula za linearnu algebru.
>>> import numpy as np
>>> import numpy.linalg as la
>>> v1 = np.array([1,1,2])
>>> v2 = np.array((2,2,4))
Množenje skalarom je prirodno:
>>> 3*v1
array([3, 3, 6])
Skalarni i vektorski produkt vektora:
>>> np.dot(v1, v2)
12
>>> np.cross(v1, v2)
array([0, 0, 0])
Norma (“duljina”) vektora:
>>> la.norm(v2)
4.8989794855663558
Množenje matrica te množenje matrice i vektora ide na prirodan način:
>>> A = np.array([[1, 2, 1], [4, 3, 3], [9, 1, 7]]); A
array([[1, 2, 1],
[4, 3, 3],
[9, 1, 7]])
>>> np.dot(A, A)
array([[18, 9, 14],
[43, 20, 34],
[76, 28, 61]])
>>> np.dot(A, v1)
array([ 5, 13, 24])
Determinanta i inverz matrice:
>>> la.det(A) # = 7
-6.9999999999999947
>>> la.inv(A)
array([[-2.57142857, 1.85714286, -0.42857143],
[ 0.14285714, 0.28571429, -0.14285714],
[ 3.28571429, -2.42857143, 0.71428571]])
Numerika po defaultu ide s double precision točnošću, pa kad provjeravamo da li je \(A^{-1} A = 1\) ne dobijemo egzaktno jediničnu matricu:
>>> should_be_unit = np.dot(la.inv(A), A); should_be_unit
array([[ 1.00000000e+00, 1.60982339e-15, 4.44089210e-16],
[ -2.22044605e-16, 1.00000000e+00, -2.22044605e-16],
[ 0.00000000e+00, -1.33226763e-15, 1.00000000e+00]])
Moguće je zatražiti od NumPyja da ispisuje ovakve male brojeve kao nule:
>>> np.set_printoptions(suppress=True)
>>> should_be_unit
array([[ 1., 0., 0.],
[-0., 1., -0.],
[ 0., -0., 1.]])
Zadatak 1
Za datu matricu A definiramo svojstvene vektore (eigenvectors) \({\bf v}\) i njima pripadajuće svojstvene vrijednosti \(\lambda\) (eigenvalues) kao rješenja matrične jednadžbe
Pomoću funkcije la.eig() odredite svojstvene vrijednosti i svojstvene vektore matrice
i provjerite da dobivena rješenja zaista zadovoljavaju gornju jednadžbu.
Zadatak 2
Kreirajte 3x3 matricu sa slučajnim realnim brojevima između 0 i 10. Invertirajte je i pomnožite s originalnom matricom te se uvjerite da dobijete jediničnu matricu.
Dijagonalizacija matrice A je pronalaženje njenog rastava oblika
gdje je \(D\) dijagonalna matrica. Takav rastav se također izvodi fukcijom la.eig(), koja daje svojstvene vrijednosti i svojstvene vektore matrice. Naime, stupci matrice \(P\) su upravo svojstveni vektori od \(A\), a dijagonalna matrica \(D\) na dijagonali ima upravo odgovarajuće svojstvene vrijednosti:
>>> A = np.array([[3, 1], [1, 3]]); A
array([[3, 1],
[1, 3]])
>>> evs, P = la.eig(A)
>>> D = np.diag(evs); D
array([[ 4., 0.],
[ 0., 2.]])
>>> np.dot(np.dot(la.inv(P), A), P) # provjeravamo P^-1 A P = D
array([[ 4., 0.],
[ 0., 2.]])
Neke matrice nije moguće dijagonalizirati, u slučaju čega će matrice \(P\) i
\(D\) koje vraća la.eig()
i dalje zadovoljavati
ali \(P\) neće biti invertibilna:
>>> A = np.array([[1, 1], [0, 1]]); A
array([[1, 1],
[0, 1]])
>>> evs, P = la.eig(A)
>>> np.dot(A, P) - np.dot(P, np.diag(evs))
array([[ 0., 0.],
[ 0., 0.]])
Neinvertibilnost od \(P\) se ogleda u ogromnim brojevima koje dobijemo kad “invertiramo” tu matricu:
>>> print(la.inv(P))
[[ 1.00000000e+00 4.50359963e+15]
[ 0.00000000e+00 4.50359963e+15]]
Za rad sa simboličkim matricama možemo koristiti odgovarajuće objekte iz sympyja