from __future__ import unicode_literals #ékezetes betűk
Egyszerű lineáris egyenletrendszerek¶
Írjunk egy egyszerű rutint ami egy előre megadott \(A\) 2x2-es mátrix-ból és egy 2 elemű \(b\) vektorból megoldja az \(Ax=b\) lineáris egyenletet és vissza adja az \(x\) vektort! Használjuk a szimbolikus mintapélda utolsó feladatának eredményét!
def megold22(A,b):
x=array([0.0,0.0]);
determinans=A[0,0]*A[1,1]-A[0,1]*A[1,0];
x[0]=(A[1,1]*b[0]-A[0,1]*b[1])/determinans;
x[1]=(A[0,0]*b[1]-A[1,0]*b[0])/determinans;
return x
Vizsgáljuk meg a rutin viselkedését! Legyen \[A=\begin{pmatrix} \cos(\varphi) & -\sin(\varphi)\\ \sin(\varphi) & \cos(\varphi) \end{pmatrix}\] valamint \[b=\begin{pmatrix} 0 \\ 1 \end{pmatrix}\] illetve a \(\phi=\pi/6\). Azaz keressük azt a vektort amit ha elforgatunk \(30^\circ\)-al akkor egy \(y\) irányú egységvektort kapunk!
#A b iletve phi értékadás
phi=pi/6;
A=array([[cos(phi), -sin(phi)],
[sin(phi), cos(phi)]])
b=array([0,1])
Ezután már készek vagyunk meghívni a rutint!
megold22(A,b) # A megoldás
Szemfülesebb halgatók már előre is megsejthették a helyes megoldást:
array([cos(pi/3),sin(pi/3)]) #Vajon ez hogy jött ki.. ?
Vizsgáljunk meg még egy példát! Legyen \[A=\begin{pmatrix} 3 & 4\\ 9 & 12 \end{pmatrix}\] valamint \[b=\begin{pmatrix} 1 \\ 3 \end{pmatrix}\]
# Az A és b tömbök inicializálása
A=array([[3, 4],
[9, 12]])
b=array([1,3])
# a megoldás...
megold22(A,b) # Vajon van megoldás ?
Amint a fenti hibaüzenet is mutatja az algoritmusunkban bizony valahol nullával osztottunk! Ha jobban megnézzük a mátrixot akkor észre vehetjük hogy a determinánsa zérus mivel a második sor az első háromszorosa! Ilyen esetben tehát nem is várunk megoldás, hiszen ilyen esetben a vizsgált két egyenlet nem független egymástól!
Lineáris algebra rutinok¶
Oldjuk meg az előző forgatós feladatot a numpy lineáris egyenlet rendszer megoldó rutinjával is! A rutin neve most is solve(), figyelnünk kell azonban hogy ez nem ugyan az a függvény amit a szimbolikus egyenletrendszerek megoldására használtunk! Ez egy jó példája annak hogy attól függően hogy éppen milyen csomagokat töltöttünk be az import parancsal bizonyos függvények bizony másképp viselkednek! A jelen rendszerben (sagemathcloud) a numpy mindíg betöltődik automatikusan minden egyes újy notebookhoz. Hogy elkerüljük a sympy hasonló nevű rutinjával a keveredést célszerű új notebook-ot nyitni, vagy legalább újra betölteni a numpy csomagot " from numpy import * " paranccsal!
#A változók inicializálása
phi=pi/6;
A=array([[cos(phi), -sin(phi)],
[sin(phi), cos(phi)]])
b=array([0,1])
#A függvény meghívása, az eredmény a várt vektor...
solve(A,b)
A numpy a lineáris algebrában használt függvények illetve konvenciók használatának könnyítésére "matrix" típusú objektumokat használ. Vizsgáljuk meg mi a különbség sima "array" illetve "matrix" típusú tömbök között!
#Tesztelő tömbök
A=array([[1,2],[3,4]])
B=array([[11,22],[33,44]])
A*B # A "*" itt elementkénti szorzást jelent
#Tesztelő tömbök
A=matrix(array([[1,2],[3,4]]))
B=matrix(array([[11,22],[33,44]]))
A*B # A "*" itt valódi mátrix szorzás!!!!
Egy pár hasznos mátrix művelet:
A.T # Egy mátrix traszponáltja
trace(A) # A mátrix trace-e, vagy spur-ja vagyis a diagonális elemek összege
det(A) # A mátrix determinánsa
inv(A) # A mátrix inverze.. ha létezik..
A*inv(A) # ennek az egységmátrixnak kell lennie!!!
(2.422+1.342j) # numpy-ban a képzetes számokat így lehet leírni
C=A+B*1.0j # Ez egy komplex mátrix, Így lehet két mátrixot össze adni.. a + jel itt is elemenként működik!
C
C.H # Ez pedig az adjungáltja azaz a transzponált mátrix komplex konjugáltja
b=matrix(array([[1],[2]])) # ez egy 1 oszopból álló mátrix.. vagyis oszlopvektor
A*b # Amit meg tudunk szorozni balról egy mátrix-al
b=matrix(array([[1,2]])) # ez egy 1 sorból álló mátrix.. vagyis sorvektor
b*B # Amit meg tudunk szorozni jobbról egy mátrix-al
Oldjunk meg "matrix"-okal egy lineáris egyenletrenszert! Válasszunk véletlenszerű mátrixokat!
A=matrix(rand(10,10)) # Ez egy 10x10-es véletlen mátrix
b=matrix(rand(10,1)) # Ez egy 10x1-es véletlen mátrix.. egy oszlopvektor
x=solve(A,b) # A solve változatlanul működik..
A puding próbája az evés! Vizsgáljuk meg a kapott x menyire teljesíti a megoldandó problémát:
A*x-b
Ezt vártuk! Az eltérés csak numerikus szemét !!
Vizsgáljunk meg most egy olyan lineáris egyenletrendszert ahol a mátrix egy 10x10-es mátrix és mindenkét első mellékátlója -1 a többi elem pedig nulla! A \(b\) vektor elemeit pedig válaszuk a következő módon:
b=sin(2*array(range(1,11))*pi/11)
A=matrix(zeros((10,10))) # Inicializálás
for i in range(9):#A mellékátlók 1-el rövidebbek mint a mátrix mérete!
A[i,i+1]=-1.0 #A mellékátlók feltöltése -1 -el
A[i+1,i]=-1.0
A
x=solve(A,b) # A megoldás
Hasonlítsuk össze a bemenő \(b\) vektort és az \(x\)-et! A számértékeiket is!
#A vektorok "alakja"
plot(b,'o-',label='b') # A bemenő vektor
plot(x,'o-',label='x') # A kijövő vektor
legend();# ábra magyarázat
Irassuk ki az értékeket illetve az értékek hányadosát is!
for i in range(10):
print x[i],b[i],b[i]/x[i]
Úgy tűnik ez egy saját vektora a mátrixnak!!! Azaz \(b=\lambda x\) illetve az eredeti egyenlet \[Ax=\lambda x\] alakú volt!
Egy mátrix sajátértékeit és ajaátvektorait is meg tudjuk határozni! Ehez az eig() illetve eigvals() rutinokat használhatjuk!
eigvals(A) # Ez az A mátrix sajátértékeit adja, keressük meg a fent megtalált sajátértéket!
vals,vecs=eig(A) #Ez a sajátvektorokat is megadja
vals #itt vannak a sajátértékek
vecs #itt pedig a saját vektorok
A saját vektorok mátrixának i-edik oszlopa tartalmaza az i-edik sajátértékhez \(\lambda_i\) tartozó sajátvektort normált sajátvektort $x_i $. Tehát igaz az hogy \[x^\dagger_i A x_i=\lambda_i\]
vecs[:,3].H * A * vecs[:,3]-vals[3] # a 4. sajátvektor leellenőrzése
Jelenítsük meg az első 5 sajátvektort!
#itt egy ciklus ahelyett hogy ötször kiírnánk a plot parancsot
for i in range(5):
plot(array(vecs[:,i]), # a plot parancs csak array-t eszik
'o-', # pont és bogyó minden görbére
label='az '+str(i)+'. sajátvektor' # a görbék címkéje legyen "az i. sajátvektor" alakú
) #
legend();
