Octave¶
Az octave a matlab ingyenes utánzata. Egyaránt telepíthető linuxra, windowsra, és használható a jupyter rendszerben is.
A python-hoz hasonlóan nagy előnye, hogy egy utasításban tud műveleteket hatékonyan végezni adatsorok, vektorok illetve mátrixok összes elemével.
Néhány hasznos link:
- hivatalos octave dokumentáció
- ELTE-s octave dokumentáció
Alapok¶
A octave-ot a python-hoz hasonlóan használhatjuk számológépként:
1+2
3*4.0
octave-ban a "^"-al hatványozunk:
2^3 # Figyelem itt így hatványozunk!
Komplex számokat az i-betűvel jelöljük:
1i*1i
Változókat itt is az =-jellel deklarálunk:
x=7/4
x+1/x
A változók itt is lehetnek számok vagy karakterláncok:
x=1/55 # ez egy szám
s='Hy' # ez egy string
Fontos különbség a python-hoz képest hogy bool típusú változók külön nincsenek definiálva. A 0 és 1 egész számok reprezentálják a True és False értékeket:
2*2==4
Ahogy azt a fenti példák mutatjk az octave minden sor után produkál kimenetet. Egy sor kimenetét akkor nem írja ki az interpreter ha ";"-vel zárjuk. Egyébként pedig egy sorbeli utasítások vesszővel és pontosvesszővel is elválaszthatók:
x=4, y=3^2, 5*5
p=1/2; q=2*3;
Az értékadás nélküli kifejezés az 'ans' speciális változónak való értékadást generál, amit következő utasításokban felhsználhatunk. A 'who' paranccsal pedig lekérdezhetjük a már létrehozott változókat:
who
Range és for¶
A python-beli rang() függvényt az octave-ban az eleje:vege konstrukció implementálja:
(1:5)
A beadott értékek persze nem csak egészek lehetnek. Lépésközt 3 elemű ":" kifejezés középső elemében adhatjuk meg. ilyenkor a sorozat utolsó eleme az lesz, ami még nem lépett túl a limitként megadott értéken:
x = 2.48:-0.111:2
For ciklust írhatjuk egy vagy több sorba. A tabulálás itt nem számít, de hasznos. Itt is a ';'-t használjuk, ha csak a végeredményt akarjuk látni. Példaként számoljuk ki első néhány egész hatványát, összegüket és szorzatukat:
for i=1:5 i^3 endfor
s=0;
p=1;
for i=1:10
s=s+i;
p=p*i;
endfor
s
p
Számtani és mértani sorozatok pedig a linspace és logspace függvényekkel is állíthatók elő, a pythonnal azonos módon:
ari=linspace(0.5,3,5)
geo=logspace(0.5,3,5)
if¶
Az octave-ban sok más programnyelvhez hasonlóan az if konstrukció segítségével állítjuk döntések elé a számítógépet. Az alábbi példa illusztrálja az if- elseif-else konstrukció legálltalánosabb alakját, ha nincs szükség az else illetve az elseif ágakra akkor azok természetesen elhagyhatóak.
x=6
if (rem (x, 2) == 0)
printf ("x is even\n");
elseif (rem (x, 3) == 0)
printf ("x is odd and divisible by 3\n");
else
printf ("x is odd\n");
endif
Tömbök¶
Számokból vagy stringekből itt is alkothatunk tömböket, de fontos eltérés, hogy nincs külön list, tuple, array, matrix típus, amikre bizonyos műveletek a pythonban másképp működnek. Itt a műveleti jel elé tett ponttal lehet kérni az adott művelet elemenkénti változatát (ld. később). Másik eltérés, hogy itt egy tömbben egyszerre vagy csak számok, vagy csak stringek lehetnek. De ehhez nézzük előbb, hogyan lehet létrehozni elemeik felsorolásával:
Sorvektor, vagy adatsor az adatok közé tett ','-vel vagy spaceszel is képezhető.
t=[ 1, 3, 5, 7, 9]
u=[ 2 4 6 8]
Mátrix sorait ';'-vel vagy sortördeléssel kell elválasztani (extra [ ] nem szükséges).
M=[ 1, 2, 3; 2,3,2; 3 2 1]
R=[ 1 2 3
4 5 -6]
Adjungálást (transzponálás konjugálással) az aposztróf jellel végezhetünk, így oszlopvektort 3-féleképp tudunk létrehozni legegyszerűbben:
v=[1;2;-2]
w=[3,-4,5]'
b=[-5
-12
13]
[1i 2 3]' #figyelem itt van konjugalas
transpose([1i 2 3]) # ez csak sima transzponalas!
A tömbök méretét a size és length függvényekkel kérdezhetjük le.
size(R), size(v), length(v)
Azonos méretű tömbök egymással való összeadása, kivonása itt természetesen az elemenkénti művelet. Tömbnek számmal való szorzása, osztása is elemenkénti (érdekesség, hogy szám hozzáadása is).
A szorzások és hatványozások viszont a vektorszámításban tanult értelmezés szerintiek, láthatóan különböznek az elemenkénti művelettől, amit a *,/,^ jelek elé tett ponntal tudunk kérni:
M
#Ez matrix szorzás!!
M*M
# Ez elemenkénti szorzás!!
M.*M
# ez matrix hatvanyozas
M^3
# ez elemenkenti hatvanyozas
M.^3
2^M, 2.^M
Vektorok skaláris és vektoriális szorzatára pedig van külön függvény, a dot és a cross.
Nagyon fontos 3 eltérés:
- a tömbök indexe nem 0-val, hanem 1-gyel kezdődik
- tömbök eleminek kiolvasásánál kerek zárójel kell
z=v(3), c=M(1,3)
- Itt is lehet szeletelni, azonban a : utáni szám az index tényleges utolsó értékét jelöli:
M, bc=M(1,2:3), row12=M(1:2,:)
Itt is kereshetünk bennük minimumot, maximumot. Ha két eredményt kérünk, akkor megkapjuk a szélsőértékhez tartozó indexet is. Octaveban ez úgy lehetséges, hogy listát alkotunk azokból a változókból, amelyekbe az eredményeket akarjuk tenni. A szögletes zárójel tehát itt, az egyenlőségjeltől balra más jelentéssel bír, mint jobbra, amikor tömböt jelölt.
minbc=min(bc)
[maxbc,imaxbc]=max(bc)
Mátrix generáló függvények¶
Az eye, zeros, ones, rand, diag függvények hasonlóan működnek a pythonhoz:
E=eye(3), L=ones(1,5), rand(1,3), D=diag([1,10,100]), diagpart=diag(M)
Így például előállíthatjuk a kvantummechanikai $\hat a^+$ keltő operátor mátrixalakját, amiben a pozitív egészek négyzetgyökei vannak elhelyezve átlós vonalban (természetesen a végtelen méretű mátrix véges részletét hozzuk létre):
n=7; aplus = diag( sqrt(1:n-1), 1)
Még hasznos ez a függvény, ami Töplitz típusú mátrixot hoz létre, melynek az átlós vonalak mentén azonosak az elemei. A lineáris algebra anyagrészben láttunk ilyen szerkezetű mátrixot, ennek előálítását megkönnyíti. (E függvénnyel nemszimmetrikus mátrix is létrehozható, ha 2 vektort adunk be).
toeplitz( [ 2.5, 2, 1.5, 1, 0.5] )
Saját függvények¶
Írjunk például olyan függvényt, ami egy kisbolygó földközelponti és földtávolponti távolsága függvényében megadja a földközelponti pályára állításhoz szükséges sebességet. Függvények definálása a pythontól annyiban különbözik, hogy más kulcsszóval kell kezdeni, a visszatérési értéknek megfelelő változót meg kell adni az első sorban az = jel előtt, és a kezdő sor végére nem szabad kettőspontot tenni. Ezen felül azonban figyelni kell, hogy az értékadásokat ';'-vel zárjuk le, hacsak nem akarjuk, hogy a belső értékadás eredménye már a függvény végrehajtása közben kiíródjon (ilyen kiíratás persze tesztelés közben hasznos, így példaként meghagytam a félnagytengely $a$ hossszának kiírtását):
function z=v1(r1,r2)
gamma=6.67408e-11;
M=5.97237e24;
a=(r1+r2)/2
c=(r2-r1)/2;
b=sqrt(a^2-c^2);
z=sqrt(gamma*M/a)*b/r1;
end
Így kapjuk pl. a körpályához szükséges sebességet:
v1(6.371e6,6.37e6)
Így pedig a szökési sebességet:
v1(6.371e6,1e11)
Egyenletrendszer megoldása¶
Mint láttuk, ${\bf M} x = v$ alakú, azaz lineáris egyenletrendszer megoldását pythonban a solve függvény végzi el, itt ennek megfelelője az alábbi művelet:
M\v
Sajátérték-probléma¶
A mátrix sajátértékeit pedig így kapjuk meg:
Lambda=eig(M)
Ha a sajátvektorokat is meg akarjuk kapni, akkor két mennyiséget kell a függvénynek kiadnia. Octaveban ez úgy lehetséges, hogy listát alkotunk azokból a változókból, amelyekbe a sajátvektorokat és sajátértékeket akarjuk tenni. A szögletes zárójel tehát itt, az egyenlőségjeltől balra más jelentéssel bír, mint jobbra, amikor tömböt jelölt.
[V,D]=eig(M)
A pythontól eltérően az első mátrixban vannak a sajátvektorok, a másodikban a sajátértékek. A sajátvektorok itt is az oszlopokat alkotják, a sajátértékek pedig a vektoroknak megfelelő sorrendben következnek a második mátrix átlójában.
Szinguláris-érték dekompozíció¶
Ez egy hasonló probléma a sajátérték-kereséshez. Pl. arra vagyunk kíváncsiak, hogy egy mátrix legfeljebb milyen hosszúra tud nyújtani egy $v$ egységvektort, ha ${\bf M} v$ szorzás szerint hat rá. Ennek a szemléletes jelentését látjuk, ha felrajzoljuk a síkban az összes egységvektor végpontját, az az egységkör (3D-ben gömb). Minden pontjának vesszük a képét az előbbi szorzással véve, ezek általánosságban ferde tengelyű ellipszist alkotnak (3D-ben ellipszoidot). Tehát a legjobban megnyújtott vektor az, ami az ellipszis nagytengelyének végébe képződött, a nyújtás aránya pedig az ellipszis fél nagytengelye. A különbség a sajátvektor-kereséshez képest, hogy most nem kötöttük ki, hogy a vektor az irányát megtartsa. Könnyen belátható, hogy ez a nyújtási arány éppen az ${\bf M}^+{\bf M}$ mátrix legnagyobb sajátértékének négyzetgyöke, a keresett vektor pedig a hozzá tartozó sajátvektor. Általánosabban pedig a legkisebb, legnagyobb nyújtási arányokat (és 2-nél több dimenzióban az aránynak az irány függvényében található nyeregpontjait) e mátrix sajátértékeinek négyzetgyökei adják, a hozzájuk tartozó szinguláris vektorok pedig a megfelelő sajátvektorok. Ezeket adja meg nekünk az svd függvény, a szinguláris értékek az S mátrixba kerülnek, a szinguláris vektorok pedig a V oszlopaiba (egyúttal a mátrix ${\bf M}={\bf U}{\bf S}{\bf V}^+$ szorzatra bontódik):
[U,S,V]=svd(M)
Közel szinguláris mátrixok¶
Próbáljunk ki valamit octaveban, amit persze pythonban is kiszámolhattunk volna. Ha kimértük egy potenciál értékét $x_1,x_2,\ldots x_n$ pontokban, arra pl. interpolálás, extrapolálás érdekében illeszthetünk n-1 fokú polinomot. Ez persze megy a polyfit v. curvefit függvénnyel, de hogy a probléma lényegét lássuk, írjuk fel, hogy ez a következő egyenletrendszer megoldását jelenti: $$ \sum{j=1}^n a_j x_i^j = v_i, \;\; i=1,2,\ldots n$$ Ez ${\bf M}\cdot a=v$ egyenletrendszer alakjába írható, ahol az ismeretlenek most az $aj$ polinomegyütthatók, az egyenletrendszer együtthatómátrixának elemei pedig: $M{ij}=x_i^j$. Ha pl. $x_i=1,2,\ldots 5$, és a mért értékeke az alábbiak, akkor így állíthatjuk elő az egyenletrendszer mátrixát:
n=5; X=(1:n)'*ones(1,n); P=ones(n,1)*(0:n-1); M=X.^P
# gyártunk egy példa adatsort ilyen módon:
v=M*[1,1/2,1/3,1/4,1/5]';
Az ${\bf M}\cdot a=v$ egyenletrendszer megoldását pedig így kapjuk az itt szereplő mért $v_i$ értékekhez:
v=[2.2833;8.5333;28.4500;75.5333;168.0833];
a = M\v; a'
Ha azonban az egyik v értéket csupán 0.03 hibával mérjük, jelentős eltérést kapunk az eredményben:
a = M \ ( v + [ 0; 0; 0.03; 0; 0] ); a'
Ez a mátrix közel szinguláris jellegének köszönhető, amit könyebb megérteni 2x2 mátrix esetére. Készítsünk belőlük olyat, aminek a 2. oszlopa kis eltéréssel az 1. oszlop számszorosa:
M2=[ 2, 4
3, 6.01]
Az $\left(\begin{array}{r}1\\0\end{array}\right)$, $\left(\begin{array}{r}0\\1\end{array}\right)$ bázisvektorok képei pedig épp a mátrix két oszlopa. Így a korábbi ábra jobb felső ellipszisében a nyilak közel egyirányúak. Emiatt az ellipszis nagyon megnyúlt, azaz a két szinguláris érték hányadosa nagy. Az $\bf{M}\cdot x=v$ egyenletrendszer $x$ megoldását keressük, és $v$, pl. mérési hiba miatt megváltozik $v+dv$-re. és a megváltozott megoldást $x+dx$-szel jelöljük, akkor nyilván $\underline M\cdot \underline dx=\underline dv$ is igaz. Ha $v$ közel van az ábra egyik nyilának irányához, akkor láthatóan kis megváltozással viszonylag közel kerülhet a másik nyílhoz, így a hozzá tartozó $x$ iránya viszonylag sokat változik. Tehát a hiba lényegesen megnő.
Általánosságban is próbáljuk megérteni! Ha a megoldást megváltozását $dx$-szel jelöljük, akkor nyilván $\bf{M}\cdot dx=dv$ is igaz. Az $x$ben tapasztal hibát, pontosabban annak $|dx|/|x|$ relatív mértékét akarjuk viszonyítani a mért értékek $|dv|/|v|$ relatív hibájához. Az arány nyilván akkor a legnagyobb, ha $|dv|/|dx|$ a lehető legkisebb, és $|v/|x|$ a lehető legnagyobb, azaz $x$ és $dx$ a megfelelő szinguláris vektorok irányához közel állnak, s az utóbbi arányokat ekkor a szinguláris értékek adják. Így a relatív hiba növekedésének maximális aránya a maximális és minimális szinguláris érték hányadosa. Ezt hívjuk kondíciós számnak, $${\rm cond}(\bf{M})=s_{max}/s_{min}$$ Így ez az ${\bf M}^+{\bf M}$ sajátértékeivel is kifejezhető. Nagy számérték itt is azt jelenti, hogy a mátrix oszlopai között közelítőleg teljesül egy lineáris kapcsolat. Az octaveban egyenletrendszer megoldásakor érdemes ellenőrizni nagyságát, hogy az eredmény pontosságára következtessünk:
cond(M2)
A polinom-illesztést leíró mátrixra szintén nagy szám adódik, sőt ez a mátrix méretével gyorsan növekszik:
cond(M)