Pokazano je nekoliko standardnih primjera korištenja numpy i sympy modula u rješavanju zadataka iz vektorskih prostora.
import platform
platform.platform()
platform.python_version()
import numpy as np
np.__version__
import numpy.linalg as LA
import sympy as sp
sp.init_printing()
sp.__version__
U egzaktnoj aritmetici je $\frac{7}{3}-\frac{4}{3}-1=0$, dok u konačnoj aritmetici dvostruke preciznosti zbog grešaka zaokruživanja rezultat ne ispada nula, nego je broj jako blizu nule.
a=7/3-4/3-1
a
a==0
a==0.0
U egzaktnoj aritmetici je $\frac{10}{3}-\frac{7}{3}-1=0$, ali i u konačnoj aritmetici dvostruke preciznosti dobivamo točni rezultat. Sve ovisi o pripadnoj binarnoj reprezentaciji realnih brojeva i broju bitova koji se čuvaju za mantisu i eksponent.
b=10/3-7/3-1
b
b==0
Gornja dva jednostavna primjera pokazuju da moramo biti iznimno oprezni kada radimo u konačnoj aritmetici računala. Testiranje je li rezultat nekog numeričkog izračuna jednak nula može naš program odvesti u potpuno krivom smjeru ako nismo dobro upoznati s problemom i ako samo testiranje ne obavimo na inteligentniji način.
if a==0:
print('Rezultat je nula')
else:
print('Rezultat nije nula')
if np.abs(a) < 1e-12:
print('Rezultat je nula')
else:
print('Rezultat nije nula')
Sympy modul omogućuje nam da u pythonu provodimo izračune u egzaktnoj aritmetici. Takvi izračuni su bitno sporiji jer ne koriste direktno hardversku aritmetiku računala, već dodatno implementirane algoritme koji troše znatno više resursa.
sp.Rational(7,3)-sp.Rational(4,3)-1
sp.Rational(10,3)-sp.Rational(7,3)-1
Numpy koristi hardversku aritmetiku računala u dvostrukoj preciznosti na 64-bitnom računalu, Sympy koristi egzaktni simbolički račun koji je sporiji od hardverske aritmetike.
np.sqrt(2)+np.sqrt(5)
br=sp.sqrt(2)+sp.sqrt(5)
br
sp.N(br,50)
Matrica $A$ je singularna matrica pa nema inverznu matricu. Pomoću modula sympy to se lagano provjeri.
A=sp.Matrix([[4,-2,2],[-4,-3,3],[-9,-3,3]])
A
sp.det(A)
A.inv()
Međutim, numpy modul zbog grešaka zaokruživanja daje potpuno drukčije rezultate.
A=np.array([[4,-2,2],[-4,-3,3],[-9,-3,3]])
LA.det(A)
LA.det(A)==0
Čak dobivamo i inverznu matricu
LA.inv(A)
To nije bug u numpy modulu, to je jednostavno posljedica konačne aritmetike. Trebamo li onda izbjegavati korištenje numpy modula? Naravno da ne trebamo, jer bismo onda trebali izbjegavati hardversku aritmetiku računala i raznorazne kalkulatore, excel i slične programe. Ovi jednostavni primjeri samo pokazuju da moramo biti oprezni s rezultatima koje daje računalo, a ne ih uzimati samo zdravo za gotovo. Konkretno na ovom primjeru matrice, vidimo da je sama determinanta reda veličine $10^{-15}$, što je zaista broj jako blizu nule i već iz tog podatka možemo reći da je zadana matrica najvjerojatnije singularna. Čak ako i u stvarnosti nije možda singularna, bolje ju je proglasiti numerički singularnom jer ćemo na taj način napraviti manje štete. Također, brojevi u navodnoj inverznoj matrici su reda veličine $10^{15}$, što su iznimno veliki brojevi s obzirom na jako male brojeve u početnoj matrici pa nam već to mora sugerirati da sumnjamo u točnost rezultata koje je dalo računalo numeričkim putem. Numerički alati poput Octave će tom prilikom obavijestiti korisnika da je početna matrica najvjerojatnije singularna i da dobivenu inverznu matricu treba zanemariti.
Dobar programer može napraviti veliku štetu ukoliko se upusti u implementaciju numeričkih algoritama ako pritom ne poznaje iznimno detaljno područje numeričke matematike. Implementacije numeričkih algoritama treba prepustiti ljudima koji su veliki stručnjaci u području numeričke matematike: dakle, prije svega moraju dobro znati matematiku, a nakon toga još iznimno jako dobro znati numeričku matematiku. Stoga je uvijek najbolje koristiti gotove numeričke biblioteke i alate ukoliko pojedinac nije detaljno upoznat s numeričkom matematikom, nego se upuštati u izradu vlastitih implementacija jer smo vidjeli da i na sasvim banalnom primjeru sve može krenuti u pogrešnom smjeru. Štoviše, ne mora dobiveni rezultat uopće izgledati čudno, a može biti u potpunosti pogrešan. Kako bismo izbjegli takve situacije, najbolje je koristiti gotove numeričke biblioteke i alate. Jedne od takvih numeričkih biblioteka su numpy i scipy za programski jezik python.
Pokazano je na nekoliko zadataka kako ih brzo riješiti pomoću dostupnih naredbi iz numpy i sympy modula. Kao što je bilo rečeno ranije, numpy koristi hardversku aritmetiku računala, a sympy koristi implementirani egzaktni simbolički račun.
U $\mathcal{P}_3(t)$ zadani su polinomi $p_1(t)=t^2+t,$ $p_2(t)=t^2-2t+3.$ Prikažite, ako je moguće, polinome $p_3(t)=-t^2+8t-9$ i $p_4(t)=t+2$ kao linearne kombinacije polinoma $p_1$ i $p_2$.
A=sp.Matrix([[1,1,0],[1,-2,3],[-1,8,-9],[0,1,2]])
A=A.T
A
Iz reducirane ešalonske forme vidimo da se polinom $p_3$ može prikazati kao linearna kombinacija polinoma $p_1$ i $p_2$, dok se polinom $p_4$ ne može prikazati kao linearna kombinacija polinoma $p_1$ i $p_2$.
A.rref()
$p_3(t)=2\cdot p_1(t)-3\cdot p_2(t)$
sp.var('x y')
sp.linsolve((A[:,:2],A[:,2]),x,y)
Polinom $p_4$ nije linearna kombinacija polinoma $p_1$ i $p_2$.
sp.linsolve((A[:,:2],A[:,3]),x,y)
Neka je $W$ potprostor od $\mathbb{R}^5$ razapet s vektorima
$u_1=(1,2,-1,3,4),\ $ $u_2=(2,4,-2,6,8),\ $ $u_3=(1,3,2,2,6),\ $ $u_4=(1,4,5,1,8),\ $ $u_5=(2,7,3,3,9).$
Odredite jednu bazu i dimenziju vektorskog prostora $W$.
M=sp.Matrix([[1,2,-1,3,4],[2,4,-2,6,8],[1,3,2,2,6],[1,4,5,1,8],[2,7,3,3,9]])
M=M.T
M
Jedna baza za vektorski prostor $W$
M.columnspace()
$\dim{W}=3$
M.rank()
Možemo sve dobiti i iz reducirane retčane ešalonske forme
M.rref()
Pomoću numpy modula možemo odrediti dimenziju prostora $W$. Za određivanje neke baze pomoću numpy modula mogli bismo koristiti određene numeričke metode, ali ovdje nećemo dublje ulaziti u numeričku linearnu algebru.
M=np.array([[1,2,-1,3,4],[2,4,-2,6,8],[1,3,2,2,6],[1,4,5,1,8],[2,7,3,3,9]])
M
LA.matrix_rank(M)
Odredite dimenziju i jednu bazu vektorskog prostora $R$ svih realnih rješenja homogenog sustava linearnih jednadžbi
\begin{align*} x+2y+z-3t&=0\\ 2x+4y+4z-t&=0\\ 3x+6y+7z+t&=0 \end{align*}
i nadopunite dobivenu bazu do baze za $\mathbb{R}^4$.
A=sp.Matrix([[1,2,1,-3],[2,4,4,-1],[3,6,7,1]])
A
Jedna baza za $R$
baza=A.nullspace()
baza
$\dim{R}=2$
len(baza)
Jedna nadopuna do baze za $\mathbb{R}^4$: $\big\{(-2,1,0,0),\big(\frac{11}{2},0,-\frac{5}{2},1\big),(1,0,0,0),(0,0,1,0)\big\}$
M=sp.Matrix([list(baza[0]),list(baza[1]),[1,0,0,0],[0,1,0,0],[0,0,1,0],[0,0,0,1]])
M=M.T
M
M.rref()
Pomoću numpy modula možemo odrediti $\dim{R}$ jer je ta dimenzija jednaka broju parametara u općem rješenju homogenog sustava.
A=np.array([[1,2,1,-3],[2,4,4,-1],[3,6,7,1]])
A
$\text{broj parametara}=\text{broj nepoznanica}-r(A)$
4-LA.matrix_rank(A)
U $M_2(\mathbb{R})$ zadan je skup $\mathcal{B}=\bigg\{\begin{bmatrix}2&3\\ 0&1\end{bmatrix}, \begin{bmatrix}5&4\\ -1&-1\end{bmatrix}, \begin{bmatrix}2&2\\ 1&1\end{bmatrix}, \begin{bmatrix}8&1\\ -3&4\end{bmatrix}\bigg\}.$
Dokažite da je $\mathcal{B}$ baza za $M_2(\mathbb{R})$ i odredite koordinate matrice $A=\left[\begin{smallmatrix}2&0\\ 0&8\end{smallmatrix}\right]$ u bazi $\mathcal{B}$.
B=sp.Matrix([[2,3,0,1],[5,4,-1,-1],[2,2,1,1],[8,1,-3,4]])
B=B.T
B
$\mathcal{B}$ je baza za $M_2(\mathbb{R})$ jer je matrica $B$ maksimalnog ranga
B.rank()
Ili pomoću determinante, $\mathcal{B}$ je baza za $M_2(\mathbb{R})$ jer je $\det{B}\neq0.$
B.det()
Inverzna matrica matrice $B$
B.inv()
Koordinate matrice $A$ u bazi $\mathcal{B}$
v=B.inv()*sp.Matrix([2,0,0,8])
v
sp.N(v)
B=np.array([[2,3,0,1],[5,4,-1,-1],[2,2,1,1],[8,1,-3,4]])
B
$\mathcal{B}$ je baza za $M_2(\mathbb{R})$ jer je matrica $B$ maksimalnog ranga
LA.matrix_rank(B)
Ili pomoću determinante, $\mathcal{B}$ je baza za $M_2(\mathbb{R})$ jer je $\det{B}\neq0.$
LA.det(B)
Inverzna matrica matrice $B$
T=LA.inv(B)
T
Koordinate matrice $A$ u bazi $\mathcal{B}$
u=np.array([2,0,0,8])
np.matmul(u,T)
Ako transponiramo formulu $Y=BX$, tada dobivamo $Y^T=X^TB^T$. Kako su $X$ i $Y$ jednostupčane matrice, tada su $X^T$ i $Y^T$ jednoredne matrice. Isto tako, stupci matrice $B$ postaju retci matrice $B^T$. Drugim riječima, to znači da kod formiranja matrice prijelaza $T$ možemo koordinate pripadnih vektora umjesto u stupce pisati u retke matrice prijelaza $T$. U tom slučaju formula za računanje koordinata vektora umjesto oblika $Y=TX$ poprima oblik $Y=XT$ pri čemu su $X$ i $Y$ jednoredne matrice. Upravo smo tu činjenicu koristili kod računanja koordinata vektora pomoću numpy modula pa nismo trošili vrijeme na transponiranje matrica.
Koji ćemo pristup odabrati je stvar ukusa, ovisno o tome želimo li radije koordinate vektora pisati u jednostupčane ili jednoredne matrice. Oba pristupa daju na kraju iste koordinate, razlika je samo u zapisu vektora i formi formule.
Ako nam se više sviđa stupčani zapis vektora, možemo to isto realizirati u numpy modulu. Međutim, s obzirom na unošenje vektora u numpy modulu, jednostavnije nam je preko formule koja koristi retčani zapis vektora.
B
B=B.transpose()
B
T=LA.inv(B)
T
u=np.array([[2],[0],[0],[8]])
np.matmul(T,u)
Kod množenja matrica, jednodimenzionalnu listu numpy zapravo prilagođava kao retčanu ili stupčanu matricu, ovisno o tome nalazi li se ta lista lijevo ili desno od matrice s kojom se množi. Stoga možemo vektor pisati kao jednodimenzionalnu listu i koristiti stupčani zapis linearnog operatora. Na izlazu se opet vraća jednodimenzionalna lista.
u1=np.array([2,0,0,8])
np.matmul(T,u1)
U $\mathcal{P}_4(t)$ zadani su polinomi
$$\begin{aligned}p_1(t)&=9t^3-13t^2-t+3\\ p_2(t)&=17t^2-6\\ p_3(t)&=6t^3+12t^2+26t+17\\ p_4(t)&=7t^3-4t^2-21t-6\\ q_1(t)&=14t^3+5t^2+6t-1\\ q_2(t)&=9t^3-16t^2-5t+17\\ q_3(t)&=-21t^3+6t^2+2t-10\\ q_4(t)&=-23t^3+15t^2+15t-22\end{aligned}$$
Odredite matricu prijelaza iz baze $\mathcal{B}_2=\big\{q_1,q_2,q_3,q_4\big\}$ u bazu $\mathcal{B}_1=\big\{p_1,p_2,p_3,p_4\big\}.$
$\mathcal{B}_{\tiny\text{kan}}\stackrel{\ T_1\ }{\longrightarrow}\mathcal{B}_1,\quad \mathcal{B}_{\tiny\text{kan}}\stackrel{\ T_2\ }{\longrightarrow}\mathcal{B}_2,\quad\mathcal{B}_2\xrightarrow{\ T_2^{-1}T_1\ }\mathcal{B}_1$
T1=sp.Matrix([[9,-13,-1,3],[0,17,0,-6],[6,12,26,17],[7,-4,-21,-6]])
T1=T1.T
T2=sp.Matrix([[14,5,6,-1],[9,-16,-5,17],[-21,6,2,-10],[-23,15,15,-22]])
T2=T2.T
T1,T2
Matrica prijelaza iz baze $\mathcal{B}_2$ u bazu $\mathcal{B}_1$
T=T2.inv()*T1
T
sp.N(T)
T1=np.array([[9,-13,-1,3],[0,17,0,-6],[6,12,26,17],[7,-4,-21,-6]])
T1=T1.transpose()
T2=np.array([[14,5,6,-1],[9,-16,-5,17],[-21,6,2,-10],[-23,15,15,-22]])
T2=T2.transpose()
T1
T2
Matrica prijelaza iz baze $\mathcal{B}_2$ u bazu $\mathcal{B}_1$
np.matmul(LA.inv(T2),T1)