Matlab Hogyan. c GPL. BME Hidrodinamikai Rendszerek Tanszék február 17. belső használatra

Matlab Hogyan c

GPL

BME Hidrodinamikai Rendszerek Tanszék bels˝o használatra

2005. február 17.

Tartalomjegyzék 1. Alapok

5

1.1. Az alapbeállítás menürendszere . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

5

1.2. Parancssor, scriptek, függvények . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

6

1.3. Script fájlok . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

6

1.4. Függvények . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

6

1.5. Lokális és globális változók . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

7

1.6. Függvény mutatók . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

8

2. Muveletek ˝ mátrixokkal, vektorokkal

10

2.1. Alapm˝uveletek . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

10

2.2. Lineáris Algebra . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

11

3. Ciklusutasítások, elágazások

13

3.1. for . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

13

3.2. while . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

13

3.3. if . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

14

3.4. switch . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

15

4. Nemlineáris feladatok

16

4.1. Polinomok zérushelyei . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

16

4.2. Optimalizációs feladatok . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

16

4.3. Nemlineáris függvények zérushelyei . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

19

5. I/0

21

5.1. Formázott kivitel (output) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

21

5.1.1. Mátrix kiírása fájlba . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

22

5.2. Formázott beolvasás (intput) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

23

5.2.1. Adatfájlból beolvasás . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

23

6. Grafika

25

1

TARTALOMJEGYZÉK

7. Interpoláció

2

27

7.1. Egydimenziós interpoláció . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

27

7.2. Többdimenziós interpoláció . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

27

8. Görbeillesztés

28

9. Közönséges differenciálegyenletek (ODE) megoldása

29

9.1. Kezdeti érték feladatok (IVP) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

29

9.1.1. Szabadsugár sebességeloszlása . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

31

9.1.2. Egy merev feladat: a Van der Pol egyenlet . . . . . . . . . . . . . . . . . . . . . . . . . . . .

32

9.1.3. Eseménykezelés: vízcsepp alakja . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

33

9.2. Perem érték feladatok (BVP) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

33

9.2.1. Couette áramlás sebességgradienssel . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

34

10. Fourier analízis

37

11. Példaprogramok

39

11.1. Karakterisztikák módszere . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

39

11.2. Nyíltfelszín˝u csatornaáramlás . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

42

11.3. Függ˝olegesb˝ol vízszintesbe vezet˝o ívben mozgó anyagdugó . . . . . . . . . . . . . . . . . . . . . . .

44

Tárgymutató Fourier analízis fft_demo.m , 36 függvény kilincsek handle_pelda.m , 7 függvények masodfok_fv.m , 6 for demo_for.m , 12 formázott i/o demo_io_1.m , 21 demo_io_2.m , 22 görbeillesztés polyfit_pelda.m , 27 grafika demo_2d_grafika.m , 24 demo_3d_grafika.m , 25 kezdeti érték feladatok csepp.m , 32 szabadsugar.m , 30 van_der_Pol.m , 31 lokális és globális változók masodfok_fv2.m , 6 masodfok_fv3.m , 6 nemlineáris egyenleternedszerek ventillator.m , 18 optimalizálás fminbnd_demo.m , 16 fminsearch_demo.m , 17 peremérték feladatok couette.m , 34 nyiltfelszin.m , 38 scriptek masodfok_script.m , 5 while stirling.m, 13

3

TÁRGYMUTATÓ

4

El˝oszó Ez a rövid összefoglaló a Budapesti M˝uszaki és Gazdaságtudományi Egyetem Hidrodinamikai Rendszerek Tanszékén (volt Vízgépek Tanszék) készült a Matlab programcsomag oktatásának és használatának megkönnyítésére. Igyekeztem nagyon egyszer˝u példákkal kezdeni és onnan fokozatosan bemutatni a program által kínált lehet o˝ ségeket, ugyanakkor sokkal inkább ugródeszkát és ötleteket szeretnék adni, mintsem egy magyar nyelv˝u dokumentációt. Mivel nem vagyok profi felhasználó, esetleg egyes programozástechnikai megoldásaim nem elegánsak/szépek/gazdaságosak, ezért nagyon fontosnak tartom a beépített Help használatát, ami szerintem annyira alapos és érthet o˝ , hogy akár azon keresztül is gyorsan el lehet sajátítani a programnyelvet. A példaprogramokkal kapcsolatban szívesen fogadok minden észrevételt, javaslatot. Ezen írásos anyag és a hozzá tartozó kiegészíto˝ és magyarázó példaprogramok letöltheto˝ k a www.vizgep.bme.hu oldalról nyíló link segítségével. Meg szeretném jegyezni, hogy a Matlab scriptek nagy része kompatibilis a Unix/Linux rendszerek alatt hozzáférheto˝ és ingyenes Octave programcsomaggal, mely letölthet o˝ pl. a www.octave.org honlapról. A dokumentáció magyarosított LATEX-ben készült, ehhez további információk találhatók pl. a következ o˝ weblapon: http://ebizlab.hit.bme.hu/ tisanyi/latex/ Az anyaggal kapcsolatos minden észrevételt, hibát, stb. szívesen fogadok. H˝os Csaba [email protected] 2004. június 27.

Szeretném megköszönni Kiss Atillának és Pandula Zolinak a hasznos észrevételeket! Igyekeztem tanácsaik alapján javítani a magyarázatokon ill. kijavítottam jópár hibát (de azért még hagytam benne....). Újdonságak számít a Pneumatikus anyagszállítás példaprogram, Dr. Váradi Sándor tárgyához segítség. 2005. február 15.

1. fejezet

Alapok A program neve a Matrix Laboratory rövidítése, ami nagyon találó, ugyanis az egész program mátrix-struktúrákban gonolkodik. Ez alatt nem csupán a numerikus mátrixokra gonolunk hanem pl. egy grafikus ábra beállításai is ilyen struktúrában kerülnek tárolásra. Ez elso˝ hallásra furcsának t˝unhet, de (a) ha kello˝ gyakorlatra teszünk szert, nagyban megkönnyíti a munkát és (b) úgysem tudunk mást csinálni, ilyenre írták a programot és kész. A csomag maga egy matematikai eljárásgy˝ujteményként is felfogható, ahol a felhasználónak kell a konkrét kódot megírnia, de mikor pl. egy sajárértékfeladathoz ér, egyszer˝uen meghív egy eljárást. Ebb o˝ l fakad az a tulajdonság is, hogy a Matlab f˝o er˝ossége a numerika, analitikus számítások elvégzésére nem kényelmes (arra pl. ott a Mathematica vagy a Maple rendszer). A program másik hatalmas el˝onye a fejlett grafikai rendszer, mellyel könnyen és gyorsan elkészíthet o˝ k a grafikonok, ábrák és ezek gyakorlatilag minden fontosabb formátumban (.ps, .eps, .jpg, .gif, stb.) menthet o˝ k. Elöljáróban még annyit érdemes tudni a Matlab-ról, hogy a technikai numerikus számításoknak pillanatnyilag ez a nemzetközileg legjobban elfogadott és elterjedt programnyelve.

1.1. Az alapbeállítás menürendszere A Matlab indításakor megjeleno˝ menürendszer a következ˝o elemekb˝ol áll (alapbeállításban): • Command Window: itt adjuk ki a parancsokat, pl. egyszer˝ubb számítások, függvényhívások vagy programok futtatása. • Launch Pad: hozzáférhetünk az installált Toolbox-okhoz (amik általában külön megvásárolható ’extra’ eljárásgy˝ujtemények.) • Workspace: nyomon követhetjük az aktuális változókat, azok nagyságát, esetleg törölhetjük o˝ ket. Ez komolyabb programok futtatása esetén lehet fontos, ha tudni akarjuk a Matlab által lefogalt memóriaterületet. • Command History: az eddig kiadott parancsok. • Current Directory: aktuális könyvtár. Ha esetleg más könyvtárban elhelyezett eljárást, adatfile-t szeretnénk használni, a File ⇒ Set Path... menüben hozzáadhatjuk a saját könyvtárainkat. Az ablakkiosztást megváltoztathatjuk a View ⇒ Desktop Layout paranccsal (vagy kézzel). A paranccsorba máris begépelhetjük: » 2*2, mire a program boldogan válaszol: ans = 4 5

1. FEJEZET. ALAPOK

6

1.2. Parancssor, scriptek, függvények A Matlab-bal végzett számítások során - a nagyon egyszer˝u, zsebszámológép szint˝u számításoktól eltekintve - programozói szemlélet szükséges, azaz az ember gyakorlatilag egy kódot ír. A bonyolultságtól függ o˝ en (elméletileg) háromféle megközelítése lehet egy számítási feladatnak: • Parancssor (ld. el˝oz˝o példa) • Script: egymás után, lépésr˝ol lépésre végrehajtandó parancsok, külön fájlban tárolva. • Egymásba ágyazott függvények, eljárások, stb., hasonlóan más programnyelvekhez (pl. Pascal, C, C++, stb.).

1.3. Script fájlok Ez egy lineárisan, lépésr˝ol lépésre végrehajtandó utasítássorozatot tartalmazó fájl. Példaként lépjünk be a Matlab saját szerkeszt˝ojébe (File ⇒ New ⇒ M-file) és készítsünk programot a másodfokú egyenlet megoldására. masodfok_script.m a=1; b=2; c=1; D= b^2 - 4*a*c; if D>0 uzenet=’valos gyokok:’; elseif D<0 uzenet=’komplex gyokok’; else uzenet=’ketszeres multiplicitasu gyok, Paal tanar ur kedveert’; end gyok1=(-b+sqrt(D))/2/a; gyok2=(-b-sqrt(D))/2/a; disp(uzenet); disp(gyok1); disp(gyok2); Magyarázat: Az els˝o sor értékadás a, b és c változóknak; a Matlab-ban nem kell deklarálni a változók típusát. Ha kett o˝ spontot teszünk egy sor végére, a program futtatásakor nem jelenik meg az aktuális sor kimenet a parancsablakban. Az if ciklus eldönti, a gyökök valósak-e vagy komplexek és az uzenet változó ennek megfelel o˝ en string típusú értéket kapott. Következik a két gyök (gyok1, gyok2) kiszámítása, majd a disp(x) paranccsal megjelenítjük az x változó értékét (ami bármilyen típus lehet). Mentsük el a fájlt pl. masodfok_script.m-nek , majd lépjünk vissza a Command Window-ra és gépeljük be: >> masodfok_script komplex gyokok -1.0000e+000 +1.4142e+000i -1.0000e+000 -1.4142e+000i

1.4. Függvények A függvényeket paraméterlistával hívjuk meg és visszatérési értékeik vannak. Az el o˝ z˝o másodfokú példa függvényként:

1. FEJEZET. ALAPOK

7

masodfok_fv.m function gyokok = masodfok_fv(a,b,c) D=b^2-4*a*c; gyok1=(-b+sqrt(D))/2/a; gyok2=(-b-sqrt(D))/2/a; gyokok = [gyok1 gyok2];

Magyarázat Ezt a fájlt már nem lehet akármilyen néven elmenteni, a fájlnévnek masodfok_fv.m-nek kell lennie (azaz az eljárárs nevének). A függvény bemenete a, b és c, kimenete pedig a gyokok vektor. A függvény hívása parancssorból: >> mo=masodfok_fv(1,2,5); >> mo mo = -1.0000 + 2.0000i

-1.0000 - 2.0000i

A D diszkrimináns bels˝o változó, a parancssorból nem érheto˝ el: >> D ??? Undefined function or variable ’D’.

1.5. Lokális és globális változók Induljunk ki az el˝oz˝o esetb˝ol, ahol adott egy függvény és ezen belül definiáltunk változókat. Ha most azt szeretnénk, hogy valamilyen bels˝o változóhoz másik függvény is hozzáférjen, két leheto˝ ségünk van: • a függvények argumentumlistájában definiáljuk az értéket (ld. masodfok_fv2.m) • globális változóként definiáljuk a változókat (ld. masodfok_fv3.m) A következ˝o program az argumentumlistában megkapott a, b és c változókat átadja a diszkr(a,b,c) eljárásnak, amely szintén az argumentumlistában veszi át az értékeket. A b2 és negyac véltozók csak a diszkr eljáráson belül hozzáférhet˝ok. masodfok_fv2.m function gyokok=masodfok_fv2(a,b,c) D=diszkr(a,b,c); gyok1=(-b+sqrt(D))/2/a; gyok2=(-b-sqrt(D))/2/a; gyokok = [gyok1 gyok2];

function ertek=diszkr(a,b,c)

1. FEJEZET. ALAPOK

8

b2=b^2; negyac=4*a*c ertek=b2-negyac;

A masodfok_fv3 eljárásnál az a, b és c változók globálisnak vannak deklarálva, így az összes olyan függvény hozzáfér ill. módosíthat rajtuk, amelyben szintén globálisként definiáltuk o˝ ket. (Tehát nem elég a ’f˝oprogramban’ globálisként definiálni, az összes olyan eljárásban meg kell tenni ugyanezt, ahol hozzá szeretnénk férni.) masodfok_fv3.m function gyokok=masodfok_fv3 global a b c a=1; b=2; c=3; D=diszkr; gyok1=(-b+sqrt(D))/2/a; gyok2=(-b-sqrt(D))/2/a; gyokok = [gyok1 gyok2]; function ertek=diszkr global a b c b2=b^2; negyac=4*a*c; ertek=b2-negyac;

1.6. Függvény mutatók A függvény mutatók segítségével függvényekre hivatkozhatunk. Tipikusan egy eljárás argumentumlistájában szerepelnek a változók mellett, pl. a beépített fzero(fun,x0) gyökkeres o˝ eljárás a fun változóban várja annak a függvénynek a nevét, melynek az x0-hoz közeli gyökét visszaadja. A mutatók használatának rejtelmeir o˝ l és el˝onyeir˝ol a Help b˝o felvilágosítást ad. Tipikusan akkor érdemes használni, ha sokszor kell ugyanazt az eljárást lefuttatni különböz o˝ függvényeken és nem akarjuk mindig átírni a programot. Az alábbi egyszer˝u példa (handle_pelda.m) szemlélteti a fentieket. handle_pelda.m function handle_pelda a=2; b=3; fhandle=@osszead; c=feval(fhandle,a,b) fhandle=@szoroz; c=feval(fhandle,a,b) %-----------------------function y=osszead(a,b) y=a+b;

1. FEJEZET. ALAPOK

9

function y=szoroz(a,b) y=a*b;

A program kimenete a parancssorban: >> handle_pelda c = 5

c = 6

2. fejezet

Muveletek ˝ mátrixokkal, vektorokkal 2.1. Alapmuveletek ˝ A Matlab használata során talán a mátrixm˝uveletek megismerése és elsajátítása alapvet o˝ , hiszen az egész programnyelv erre épül (a név is innen ered: Matrix Laboratory). A mátrixok elemeit szögletes zárójel fogja közre, a sorok elemeit üres karakter választja el, az új sort pedig pontosvesszo˝ jelöli: >> A = [16 3 2 13; 5 10 11 8; 9 6 7 12; 4 15 14 1] A = 16 5 9 4

3 10 6 15

2 11 7 14

13 8 12 1

A transzponálás fels˝o vessz˝ovel történik: >> A’ ans = 16 3 2 13

5 10 11 8

9 6 7 12

4 15 14 1

Egy mátrix elemeire pedig az alábbi módon hivatkozhatunk: >> A(2,3) ans = 11 Gyakran szükséges a mátrix egy sorát vagy oszlopát ’kivonni’, ez a kett˝ospont operátorral történik, a sor-oszlop konvenció figyelembevételével. Tehát a A(:,2) parancs a mátrix második oszlopát, a A(1,:) az els o˝ sorát adja eredményül. A kett˝ospont operátor jelentése lehet <-tól:lépés:-ig> is, azaz: 10

˝ 2. FEJEZET. MUVELETEK MÁTRIXOKKAL, VEKTOROKKAL

11

>> 0:pi/4:pi ans = 0

0.7854

1.5708

2.3562

3.1416

Fontos szerepet játszik még a pont operátor is, amely az elemenkénti m˝uveleteket jelenti. Tehát ha A és B két mátrix, akkor A*B a hagyományos mátrixszorzatot adja, míg a A.*B az elemenkénti szorzatot. Mátrixokat beolvashatunk adatfájlból is a load parancs segítségével. Bármilyen szövegszerkeszt o˝ ben (pl. Notepad) készítsük el az alábbi sorokat: 1 4. -5.56 344 -.454 33 -0 1 1 majd mentsük el ize.txt néven. A load parancs beolvassa az adatfile-t és a kiterjesztés levágásával kapott változóba teszi az adatokat: >> load ize.txt >> ize ize = 1.0000 344.0000 0

4.0000 -0.4540 1.0000

-5.5600 33.0000 1.0000

2.2. Lineáris Algebra A fontosabb parancsok rövid áttekintése: Parancs zeros(n,m) A+B A*B A’ det(A) trace(A) inv(A) eig(A) poly(A) A.*B A./B [L,U] = lu(A) A\b

Például: >> [V,lambda]=eig(ize);

Leírás n × m-es nullmátrixot hoz létre Összeadás (elemenként) Mátrixszorzás Transzponált Determináns Nyom (f˝oátlóbeli elemek összege, spur) Inverz Sajátértékek, sajátvektorok A karakterisztikus polinom együtthatói csökken˝u kitev˝u szerint Elemenkénti mátrixszorzat Elemenkénti osztás A mátrix LU felbontása Az Ax=b lineáris egyenletrendszer megoldása. Ugyanezt az eredményt adja az inv(A)*b parancs is. b-nek oszlopvektornak kell lennie!

˝ 2. FEJEZET. MUVELETEK MÁTRIXOKKAL, VEKTOROKKAL

12

>> ize*V(:,1) ans = -4.0124 37.6942 -0.9685 >> lambda(1)*V(:,1) ans = -4.0124 37.6942 -0.9685 Példaként számítsuk ki egy adott általános, térbeli feszültségi állapothoz tartozó f o˝ feszültségeket és a f˝ofeszültségi irányokat! A szilárdságtan tanítása szerint a T feszültségi tenzorhoz tartozó σ i f˝ofeszültségek a feszültségi tenzor sajátértékei: T v i = σ i vi , ahol vi jelöli a σi -hez tartozó sajátvektort. A feszültségi mátrix diagonalizálható egy V transzformációs mátrixszal, melynek oszlopait a sajátvektorokból képezzük: V = [v1 , v2 , v3 ] (az eig() parancs éppen ezt adja vissza):  σ1  V −1 T V =  0  0

Valóban, >> inv(V)*ize*V ans = -37.9196 0.0000 -0.0000

-0.0000 37.1044 -0.0000

0.0000 -0.0000 2.3611

0 σ2 0

 0   0  .  σ3

3. fejezet

Ciklusutasítások, elágazások 3.1. for Egy utasítás ismétlése megadott alkalommal. demo_for.m n=4; a = zeros(n,n); for i = 1:n for j = 1:n a(i,j) = 1/(i+j -1); end; end; disp(a);

A script futásának eredménye : >> demo_for 1.0000 0.5000 0.3333 0.2500

0.5000 0.3333 0.2500 0.2000

0.3333 0.2500 0.2000 0.1667

0.2500 0.2000 0.1667 0.1429

3.2. while Egy utasítás ismétlése bizonytalan számú alkalommal, valamilyen feltétel teljesüléséig. Példaként számítsuk ki π értékét eps relatív pontossággal a Stirling-formula segítségével: 1 n! en n→∞ 2n nn

π = lim

!2

.

Fontos azonban, hogy n!, nn és en értékét ne közvetlenül számítsuk ki a túlcsordulás veszélye miatt, hanem az alábbi alakban: (n − 1)e e 2e 3e n! en × ×···× × e. = × nn n n n n 13

3. FEJEZET. CIKLUSUTASÍTÁSOK, ELÁGAZÁSOK

14

Az eljárást függvényként definiáljuk kimenet nélkül, max_rel_hiba bemenettel stirling.m : stirling.m function stirling(max_rel_hiba) fprintf(’\n pi szamitasa Stirling-fromulaval \n’); eps=1.e5; n=1; while eps>max_rel_hiba a=1; for i=1:n a=a*i/n*exp(1); end; pii=a^2/2/n; eps=abs(pi-pii)/pi; n=n+1; end; fprintf(’\n n = %d \n kozelito ertek: %10.8f \n ... hiba: %10.8e \n \n’,n,pii,eps);

A futtatás eredménye: >> stirling(1e-4) pi szamitasa Stirling-fromulaval n = 1668 kozelito ertek: 3.14190677 hiba: 9.99850013e-005 Az eredmények: eps n

10−1 3

10−2 18

10−3 168

10−4 1668

3.3. if A parancsok feltételes elvégzése. Szintaxis: if kifejezés1 parancs1; parancs2; elseif kifejezés2 parancsok else parancsok end Természetesen az elseif rész kihagyható, az else sem kötelez o˝ egyszintes if esetén. A feltételvizsgálatnál több kifejezés is összekapcsolható a logikai és &, vagy | ill. a nem ~ operátorokkal. Majd ha lesz valami jó, nemtriviális ötletem példaprogramra, írok egyet.

3. FEJEZET. CIKLUSUTASÍTÁSOK, ELÁGAZÁSOK

15

3.4. switch Kapcsolás egy kifejezés több, elo˝ re megadott értéke esetén (ez általában megoldható if-el is, de akkor az eseményvizsgálat lassítja a programot). Tegyük fel pl., hogy adva van egy nap nev˝u string változónk: switch nap case {’hetfo’,’szerda’} program=’mozi’ case ’kedd’ program=’szinhaz’ case ’pentek’, program=’koncert’ otherwise, program=’Meg nem tudom...’ end C/C++ programozók figyelmébe: a Matlab-ban nem kell break parancs minden case után, csak az aktuális parancsok kerülnek hívásra.

4. fejezet

Nemlineáris feladatok 4.1. Polinomok zérushelyei Legyen adva az alábbi m-edrend˝u polinom:

p(x) =

m X i=0

a i xi = a 0 + a 1 x + a 2 x2 + a 3 x3 + · · · + a m xm .

Ekkor a fenti polinomnak pontosan m darab zérushelye van (multiplicitásokkal számolva), valósak és esetleg komplexek. A gyököket a roots(a) paranccsal kereshetjük meg, ahol a egy vektor, mely az együtthatókat tartalmazza x hatványai szerint csökken˝o sorrendben, azaz a = [am , am−1 , . . . a1 , a0 ]. A függvény kimenete szintén egy vektor, mely a gyököket tartalmazza. Példaként számítsuk ki a p(x) = −4 + (3 + 5i)x − 2ix2 + x3 polinom zérushelyeit (ahol i a képzetes egységgyök): >> a=[ 1 -2i 3+5i -4]; roots(a) ans = -1.2403 + 3.4313i 1.1345 - 0.5762i 0.1058 - 0.8551i

4.2. Optimalizációs feladatok Logikailag ez a feladat megel˝ozi a nemlineáris egyenletrendszerek zérushelyeinek kérdését, ugyanis azokat is erre fogjuk visszavezetni. Nézzük el˝ore az egy ismeretlen szerinti optimalizációt! Keressük azt az x0 értéket, melyre f (x0 ) = Min! és

x 1 ≤ x0 ≤ x2 .

(Egy maximumfeladatot analóg módon, − f (x0 ) minimumfeladataként lehet megfogalmazni.) A alkalmazandó Matlab 16

4. FEJEZET. NEMLINEÁRIS FELADATOK

17

1 0.8 0.6

y

0.4 0.2 0 −0.2

PSfrag replacements

−0.4 −0.6 0

0.5

1

x

1.5

2

4.1. ábra. Az f (x) = (x − 1)2 sin(10x) függvény és a két megtalált minimumhely.

függvény a fminbnd(fun,x1,x2), mely a fun függvény x1 és x2 közötti egyik, lokális minimumhelyét próbálja megkeresni. Példaként keressük meg a f (x) = (x − 1)2 sin(10x) függvény minimumát az 0 ≤ x ≤ 2 intervallumon! A kapcsolódó példaprogram (fminbnd_demo.m) két iterációt indít, egyet a 0 ≤ x ≤ 2 intervallumon, másikat a 1 ≤ x ≤ 2-en. Mint az a 4.1 ábrán látszik, két különböz˝o minimumhelyet talált meg az eljárás, ami er˝osen óvatósságra int a könnyelm˝u használattal kapcsolatban. A függvényhívás bal oldalán elegend˝o egy érték is (tehát pl. x1 az [x1,y1] helyett), ekkor csak minimumhely koordinátáját adja meg az eljárás. Az optimset paranccsal különbözo˝ opciókat állíthatunk be, itt az összes iterációs lépésro˝ l lekérdeztük az információt. Hasonlóan beállítható a tolerancia (TolX), a maximális iterációk száma (MaxIter) és a függvény kiértékeléseinek maximális száma (MaxFunEvals). A grafikával kapcsolatban a részleteket a 6. fejezet tartalmazza. fmin1D_demo.m function fmin1D_demo [x1,y1] = fminbnd(@fv,0,2); [x2,y2] = fminbnd(@fv,1,2,optimset(’Display’,’iter’)); fplot(@fv,[0,2]), hold on plot(x1,y1,’rx’,’Markersize’,15), hold on plot(x2,y2,’r+’,’Markersize’,15), hold off xlabel(’x’), ylabel(’y’), grid on function y=fv(x) y=(x-1)^2*sin(10*x);

A program kimenete: >> fmin1D_demo Func-count 1

x 1.38197

f(x) 0.138606

Procedure initial


18

200

z

150

100

50

rag replacements

0 1.5 1 0.5

1 0 −0.5

y

0.5

−1

x

4.2. ábra. A banánfüggvény, az iteráció kiindulópontja (piros csillag) és a megtalált minimumhely (sárga csillag). 2 3 4 5 6 7 8 9 10

1.61803 1.76393 1.8541 1.73995 1.75356 1.75369 1.75381 1.75384 1.75377

-0.173796 -0.546067 -0.22152 -0.543538 -0.549226 -0.549227 -0.549228 -0.549228 -0.549228

golden golden golden parabolic parabolic parabolic parabolic parabolic parabolic

Ha egy többváltozós függvény minimum- vagy maximumhelyeit szeretnénk megkeresni, az fminsearch eljárást kell használnunk. Az fminserach_demo program az ún. Rosenbrock-féle banánfüggvény minimumát keresi: f (x, y) = 100(y − x2 )2 + (1 − x)2 , mely a minimumkeres˝o eljárások egyik tesztfeladata. A minimumhely koordinátái (x, y) = (1, 1), az eljárást az (x0 , y0 ) = (−1.2, 1) pontból indítjuk, így annak át kell kelnie a két pont közti „csúcson” (ld. 4.2 ábra). fmin2D_demo.m function fmin2D_demo x0=-1.2; y0=1; z0=banan([x0 y0]); [x,z] = fminsearch(@banan,[x0,y0]) ezsurf(’100*(y-x^2)^2+(1-x)^2’,[-1.3,1.3,0.5,1.5]), hold on plot3(x0,y0,z0,’r*’,’Markersize’,30), hold on plot3(x(1),x(2),z,’y*’,’Markersize’,30), hold off xlabel(’x’), ylabel(’y’), zlabel(’z’), grid on function y=banan(x) y=100*(x(2)-x(1)^2)^2+(1-x(1))^2;


19

4.3. Nemlineáris függvények zérushelyei Itt megint azzal az egyszer˝u esettel kezdjük, amikor egyváltozós nemlineáris függvény zérushelyét kell megtalálnunk. A rendelkezésünkre álló beépített Matlab függvény az fzero(fun,x0), mely a fun függvény zérushelyét keresi x0 közelében. Használata teljesen analóg az fminbnd használatával. Amennyiben x0 egy vektor, pl. [x1 x2], az eljárás feltételezi, hogy fun(x1) és fun(x2) el˝ojele különböz˝o (ellenkez o˝ esetben hibaüzenettel tér vissza) és az adott intervallumon azt a pontot keresi meg, ahol fun elo˝ jelet vált. Meglep˝o módon, nincsen1 az fminsearch-el analóg eljárárás, azaz olyan függvény, mely nemlineáris egyenletrendszerek zérushelyeit keresi meg. Egy egyszer˝u trükkel azonban könnyen felhasználhatjuk az fminsearch eljárást ilyen célra. Legyen adva ugyanis az alábbi nemlineáris egyenletrendszer: F1 (x1 , x2 , . . . , xn ) = 0, F2 (x1 , x2 , . . . , xn ) = 0, .. . Fn (x1 , x2 , . . . , xn )

= 0.

Ha most valamilyen (x01 , x02 , . . . , x0n ) értékeket behelyettesítünk a fenti egyenletbe, nyilván nem fog teljesülni, hanem valamilyen nemnulla hiba lesz az eredmény: F i (x01 , x02 , . . . , x0n ) = ri . Ezek után a fenti egyenlet zérushelyének koordinátáira igaz lesz, hogy

R(x1 , x2 , . . . , xn ) =

n X

ri2 = Min! = 0.

i=1

Példaként tekintsük az alábbi differenciálegyenlet rendszert (DER), mely egy sorbakapcsolt ventillátor, cs o˝ , tartály és fojtás dinamikus viselkedését írja le:

d Q = f (Q) − ∆p = F1 (Q, ∆p) dt p d ∆p = Q − a ∆p = F2 (Q, ∆p), dt ahol a ventillátor jelleggörbéjét a 1 3 f (Q) = 0.27 + 0.18 1 + (4Q − 1) − (4Q − 1)3 2 2

!

(4.1)

függvény írja le. Ki szeretnénk számolni a DER egyensúlyi helyzeteit, melyekben a rendszer nyugalomban van, azaz minden változó id˝o szerinti deriváltja zérus. Tehát meg kell keresnünk az F 1 (Q, ∆p) = 0, F 2 (Q, ∆p) = 0 nemlineáris egyenletrendszer megoldásait. Grafikailag a megoldást a ventillátor jelleggörbe ∆p v = f (Q) és a fojtás jelleggörbe ∆p f = a12 Q2 metszéspontja adja. A kapcsolódó Matlab program ventillator.m közvetlenül a nemlineáris egyenletrendszert oldja meg az el˝oz˝oekben leírt séma szerint. ventillator.m 1

function ventillator

2 3

a=0.6;

4 5 6

1

q0=0.5; dp0=0.6; [x,z] = fminsearch(@hiba,[q0,dp0],[],a);

azaz én nem találtam meg...


20

0.8

0.7

0.6

∆p

0.5

0.4

0.3

0.2

PSfrag replacements

0.1

0

0

0.1

0.2

0.3

0.4

Q

0.5

0.6

0.7

0.8

0.9

4.3. ábra. Ventillátor- és fojtásjelleggörbe ill. a kiszámolt metszéspont.

7 8 9 10 11 12 13

xx=0:0.01:0.9; yv=jelleggorbe(xx); yf=1/a^2*xx.^2; plot(xx,yv,xx,yf,x(1),x(2),’ro’) axis([0 0.9 0 0.8]), grid on xlabel(’Q’), ylabel(’\Delta p’)

14 15

function R=hiba(x,a)

16 17 18 19 20

Q =x(1); dp=x(2); r1 = jelleggorbe(Q)-x(2); r2 = Q-a*sqrt(dp); R=r1^2+r2^2;

21 22

function dp=jelleggorbe(q)

23 24

dp=0.27+0.18*(1+1.5*(4*q-1)-0.5*(4*q-1).^3);

Külön magyarázatot igényel a 6. sorban az üres szögletes zárójel. Az a változó szeretnénk a f o˝ program 3. sorában beállítani és aztán átadni az összes többi eljárásnak is. Ezt megtehetjük az fminsearch argumentumlistájában, de ott a kezdeti értékek után az options változónak kell állnia és csak azután kerülhetnek be az egyéb extra argumentumok. Ezért, bár semmilyen opciót nem kívánunk igénybe venni, ki kell tennünk a szögletes zárójelet és csak azután kerülhet be az a változó. Ártatlannak t˝unhet még a 24. sor végén kitett pont a ∧3 elo˝ tt, azonban ez nagyon is fontos. Ugyanis ha a jellegörbe eljárást csak skalár értékkel hívnánk meg, ez a pont nem kellene, de a 9. sorban az xx vektor szerepel az argumentumban és így az eljárásnak egy vektort kell a harmadik hatványra emelnie a 24. sorban. Természetesen ezt elemenként szertnénk elvégezni, ezért ki kell tenni a pontot. A grafikával kapcsolatban a részleteket a 6. fejezet tartalmazza.

5. fejezet

I/0 Legegyszer˝ubben a sorvégi pontosvesszo˝ elhagyásával vagy a disp() paranccsal jeleníthetjük meg a parancsablakban egy változó, string, stb. értékét. Formázott adatkivitelt a C nyelvhez hasonló operátorokkal valósíthatunk meg. A formázás szintaktikája a következo˝ :

% + 5.3 e

PSfrag replacements

formázás módja

formázás kezdete zászló

mez˝oszélesség

pontosság

Fontosabb értékek: zászló

mez˝oszélesség pontosság formázás módja

+ 0

%d %e %f %g %s

A szám el˝ojelét mindig kiírja. Balra igazítja a számot a mez˝oben. A lefoglalt mez˝oben a kihasználatlan karaktereket 0-val tölti fel (és nem az alapértelmezett szóközzel). Az eredménynek lefoglalt karakterek száma. A tizedesvessz˝ot˝ol jobbra megjelenítend˝o jegyek száma. Decimális megjelenítés (egész típus számára). Exponenciális (tudományos) formátum. Lebeg˝opontos megjelenítés. Legjobb formátum az adott szám számára. String megjelenítése.

Pozicionálás mez˝ok között: \n \t

új sor tabulátor

5.1. Formázott kivitel (output) A formázott adatkivitel legfontosabb parancsa az fprintf(), amellyel mind képerny o˝ re, mind fájlba írhatunk:

21

5. FEJEZET. I/0

22

képerny˝ore: fprintf(formátum,a,b,c) Itt a formátum egy sztring, pl. ’a=%g \t b=%+10.7e vegul egy egesz: c=%d es jon a sorleptetes:\n’ fájlba: azonosító = fopen(filenév,jogok) fprintf(azonosító,formátum,a,b,c) fclose(azonosító) Az fopen() parancs megnyitja a filenév állományt a megadott jogokkal (amit el lehet hagyni, ekkor csak olvasásra). Megnyitás után az fprintf() a már leírt szintaktikával m˝uködik, de az argumentumlista elején a kimenet azonosítóját meg kell adni. Végül pedig lezárjuk a fájlt az fclose() paranccsal. A jogok lehetnek: r w a r+ w+ a+

olvasás (read), ez az alapbeállítás. file megnyitása/létrehozása írásra (write), már létez˝u file esetén a tartalom elvész. létez˝o file megnyitása/létrehozása hozzáírásra (append). file megnyitása írásra és olvasásra. file megnyitása/létrehozása írásra és olvasásra, már létez˝u file esetén a tartalom elvész. file megnyitása/létrehozása írásra és olvasásra, már létez˝u file esetén a tartalom nem vész el, az új adatokat a végéhez f˝uzi.

5.1.1. Mátrix kiírása fájlba A kovetkezo program formazva, elemrol elemre menti el M mátrix tartalmat (demo_io_1.m). demo_io_1.m % Az adatok: M= [ 1.2345 2.56 3.009 -4.01 -5.45e4 6245 76543 -453458 9 8 1.23e-5 0]; % Fajlnyitas adatfajl = fopen(’demo_io_1.dat’,’w’); % Iras tabulatorral tagolva, legalkalmasabb formatumban for i=1:length(M(:,1)) for j=1:length(M(1,:)) fprintf(adatfajl,’%g \t’,M(i,j)); end; fprintf(adatfajl,’\n’); end; % Kihagyunk ket sort.... fprintf(adatfajl,’\n\n’); % Iras fix szelesen, tabulatorral, tudomanyos formatumban

5. FEJEZET. I/0

23

for i=1:length(M(:,1)) for j=1:length(M(1,:)) fprintf(adatfajl,’%+10.5e \t’,M(i,j)); end; fprintf(adatfajl,’\n’); end; % Fajl lezarasa fclose(adatfajl);

5.2. Formázott beolvasás (intput) A legegyszer˝ubb eset az, ha a beolvasandó fájlban struktúráltan, mátrixszer˝uen, csak numerikus értékek vannak, ekkor megúszhatjuk egy load() paranccsal. Bonyolultabb esetben (pl. mérési eredmények) azonban a fájlnak fejléce van, nem egy mátrixot akarunk kiolvasni hanem több számot, stb...

5.2.1. Adatfájlból beolvasás Tegyük fel pl., hogy az adatfájlunk így néz ki (adatok.dat): adatok.dat Meres helye: BME Vizgepek labor Meres ideje: 2002.12.25 Merest vegezte: Hos Csaba Mintaveteli frekvencia: 10200 Hz adatfile felepitese: t,p1,p2 0.1 0.2 0.3 0.4

1.1 1.3 1.6 1.8

2.1 2.01 1.9 1.4

Ekkor az alabbi script segitsegevel kiolvashatunk mindent (demo_io_2.m): demo_io_2.m fp=fopen(’adatok.dat’,’r’); % A 13. karaktertol a sorvegig olvassuk a nevet: temp=fseek(fp,13,’bof’); % bof >> beginning of file hely=fgetl(fp); % az aktualis poziciotol sorvegig disp(hely); temp=fseek(fp,13,’cof’); % cof >> current position in the file datum=fscanf(fp,’%4f.%2f.%2f\n’); % a datum 3 elemu vektorban van az ev, honap, nap disp(datum’); temp=fseek(fp,16,’cof’); % cof >> current position in the file nev=fgetl(fp); disp(nev);

5. FEJEZET. I/0

temp=fseek(fp,24,’cof’); % cof >> current position in the file frek=fscanf(fp,’ %g Hz’); disp(frek); % 2 sort ugrunk for i=1:2 temp=fgetl(fp); end; % 3 oszlopu matrixot olvasunk, ameddig csak lehet (>> inf) A=fscanf(fp,’%g %g %g’,[3 inf]); A=A’; % igy lesz jo... disp(A); % lezarjuk az adatfajlt, ahogy azt illik fclose(fp);

24

6. fejezet

Grafika A grafikák készítését is áthatja a vektor-alapú szemlélet: pl. a plot(x,y) parancs az x és y vektorok által megadott pontokat rajzolja ki. Parancssorból (scriptbo˝ l) az ábra formázása elkészíthet˝o és a már elkészült ábra is tovább alakítható, de a grafikus ablakon beállított formázások újrarajzoláskor elvesznek. A rajzok elmenthet o˝ ek xxx.fig formátumban, ami a Matlab-bal bármikor megnyitható és tovább alakítható. Exportálni is lehet számos grafikai típusban (.bmp, .gif, .jpg, .eps, .ps), de ezek utólag már nem alakíthatóak (Matlab-bal). Tehát az ábrákat érdemes el o˝ ször xxx.fig formátumban teljesen elkészíteni és csak utána exportálni. Íme egy kis példa 2D grafikára (demo_2d_grafika.m) : Az y(φ)=sin(φ)cos(φ) fuggveny grafikonja 0.5

0.4

0.3

0.2

sin(φ)

0.1

0

−0.1

−0.2

−0.3

−0.4

−0.5

−pi

−pi/2

0 −π ≤ φ ≤ π

pi/2

6.1. ábra. A demo_2d_grafika.m script kimenete.

demo_2d_grafika.m % x és y vektorok elkeszitese: x = -pi:.1:pi; y = sin(x).*cos(x); % Rajzolas:

25

pi

6. FEJEZET. GRAFIKA

26

plot(x,y,’r-.’); %Beallitasok: grid; xlabel(’-\pi \leq \phi \leq \pi’); ylabel(’sin(\phi) cos(\phi)’); title(’Az y(\phi)=sin(\phi)cos(\phi) fuggveny grafikonja’) % x tengelyen a jelolo atallitasa: set(gca,’XTick’,-pi:pi/2:pi); set(gca,’XTickLabel’,{’-pi’,’-pi/2’,’0’,’pi/2’,’pi’});

Egy másik példa 3D grafikára (demo_3d_grafika.m) : demo_3d_grafika.m [X,Y] = meshgrid([-2:0.1:2]); Z = sin(X.*Y); % Elso abra rajzolasa (vonalasan) figure(1); plot3(X,Y,Z); % Masodik abra rajzolasa (feluletkent) figure(2); ezmeshc(’sin(x*y)’,[-2,2,-2,2]); xlabel(’x’); ylabel(’y’); zlabel(’z’); title(’A sin(xy) fuggveny.’);

A sin(xy) fuggveny.

1

z

0.5

0

−0.5

−1 2 2

1 1

0 0 −1 y

−1 −2

−2

x

6.2. ábra. A demo_3d_grafika.m script második ablaka.

7. fejezet

Interpoláció 7.1. Egydimenziós interpoláció Egydimenziós interpoláció esetén az yi = interp1(X,Y,xi,módszer) eljárást használjuk, ahol X és Y tartalmazza az adatokat és az xi pontokban szeretnénk az interpolált értékeket megkapni. A beépített módszerek: ’nearest’ ’linear’ ’spline’ ’pchip’ ’v5cubic’

„Legközelebbi szomszéd” interpoláció. Lineáris interpoláció (alapbeállítás). Harmadfokú spline interpoláció. Harmadfokú szakaszos Hermite interpoláció. Harmadfokú interpoláció az Matlab 5-b˝ul.

A ’nearest’, ’linear’, és ’v5cubic’ módszerek esetén interp1(X,Y,xi,módszer) NaN (Not-a-Number) üzenettel tér vissza xi azon elemeire, melyek az x értékei által meghatározott szakaszon kívül esnek. A többi módszer ez esetben extrapolál.

7.2. Többdimenziós interpoláció Kétdimenziós interpoláció esetén az zi = interp2(x,Y,z,xi,yi,módszer) eljárást használhatjuk, teljesen analóg módon az elo˝ z˝oekkel. Mindezek egyszer˝u általánosítása a vi = interp3(x,Y,z,v,xi,yi,zi,módszer) háromdimenziós és a yi = interp3(x1,x2,x3,y,x1i,x2i,x3i,módszer) sokdimenziós interpoláció.

27

8. fejezet

Görbeillesztés Az alap Matlab csomagban csak a polinomiális görbeillesztés található meg el o˝ re kész eljárás formájában. A p = polyfit(x,y,n) parancs eredménye p vektor, mely az n-edfokú polinom együtthatóit tartalmazza csökken˝u kitev o˝ szerint. Ha egyszer ezeket az együtthatókat kiszámítottuk, tetszo˝ leges xi pontokban yi = polyval(p,xi) paranccsal számíthatjuk ki a függvény értékét. A polyfit_pelda.m program az f (x) = x2 sin(5x) függvényhez illeszt poinomot a (0, π/2) intervallumon, majd kiszámítja a kapott polinom eltérését az osztáspontokban a norm parancs segítségével. polyfit_pelda.m function polyfit_pelda X=0:0.01:pi/2; Y=sin(5*X).*X.^2; for i=1:1:10 p=polyfit(X,Y,i); yy=polyval(p,X); hiba=norm(Y-yy,2); fprintf(’\nfokszam: %2d, hiba=%5.3e’,i,hiba); end

>> polyfit_pelda fokszam: 1, fokszam: 2, fokszam: 3, fokszam: 4, fokszam: 5, fokszam: 6, fokszam: 7, fokszam: 8, fokszam: 9, fokszam: 10,

hiba=9.571e+00 hiba=6.386e+00 hiba=2.374e+00 hiba=2.341e+00 hiba=8.731e-01 hiba=2.577e-01 hiba=1.375e-01 hiba=1.321e-02 hiba=1.030e-02 hiba=4.593e-04

28

9. fejezet

Közönséges differenciálegyenletek (ODE) megoldása Amennyiben közönséges differenciálegyenleteket szeretnénk numerikusan megoldani, az els o˝ lépés mindenképpen az, hogy els˝orend˝u alakra kell átírni az ún. Cauchy-átírással. Ez általánosan úgy történik, hogy ha adva van a d(n) z (n−1) (n−2) = G x; z , z , . . . , z dxn függvény, akkor új változókat vezetünk be: y1 = z, y2 = z0 = z(1) , y3 = z00 = z(2) ,...„ yn−1 = z(n−1) . Így megoldandó a     y1  F1 (x; y1 , y2 , . . . , yn ) y  F (x; y , y , . . . , y ) 1 2 n  d  2   2   ..  =  ..  dx  .   .     Fn (x; y1 , y2 , . . . , yn ) yn közönséges differenciálegyenlet rendszer (a Cauchy-átírás a példákon keresztül könyebben megérthet o˝ ). Kezdeti érték 0 0 0 feladatról beszélünk, ha az y(x0 ) = y1 , y2 , . . . , yn kezdeti értékb˝ol kiindulva kíváncsiak vagyunk y(x1 ) értékére. Peremérték feladatról beszélünk, ha y(x0 ) mellett el˝o van írva y(x1 ) is, ekkor azonban a megoldás létezéséhez általában meg kell adni a DE-ben egy szabad paramétert is. A DE jobb oldalához tartozó Jacobi mátrix elemeit következ o˝ képpen számíthatjuk: Ji j =

∂Fi (x; y1 , y2 , . . . , yn ) ∂y j

ahol i, j = 1 . . . n.

9.1. Kezdeti érték feladatok (IVP) A megoldó szintaktikája a következo˝ : [t,y] = megoldó(odefun,tspan,y0,options) ahol megoldó a numerikus módszert jelöli, odefun a DE jobb oldalát meghatározó függvényt, tspan vektor a megoldás intervallumát ([x_eleje x_vége]), y0 a kezdeti feltételek vektora és options a nem kötelez o˝ , esetleges további beállításokat tartalmazza. A numerikus megoldók az alábbiak lehetnek: NEM Merev feladatok: 29

9. FEJEZET. KÖZÖNSÉGES DIFFERENCIÁLEGYENLETEK (ODE) MEGOLDÁSA

ode45 ode23 ode113

Explicit, egylépéses, negyedrend˝u Runge-Kutta képlet. Általában ez lehet az els˝o próbálkozás. Szintén egylépéses Runge-Kutta típusú módszer. Enyhén merev feladatoknál gazdaságosabb lehet, mint az ode45. Változó rend˝u, többlépéses Adams-Bashforth-Moulton PECE megoldó. Ha nagyon pontos megoldásra van szükség vagy a DE nagyméret˝u (a jobb oldal kiszámítása drága), gazdaságosabb lehet, mint az ode45.

Merev (stiff) feladatok: ode15s

ode23s ode23t

Változó rend˝u többlépéses módszer. Ha az ode45 nem m˝uködik és gyanús, hogy ez a DE merevsége miatt van, ez lehet az els˝u nekifutás. Ezek is mind merev megoldók, a pontos dokumentáció a Helpben megtalálható.

ode23tb Az options változó fontosabb leheto˝ ségei: Hiba (pontosság) beállítása RelTol AbsTol NormControl

Relatív hiba. Abszolút hiba. A relatív hibát hibrid módon, a megoldás normájával számítja.

Megoldó kimenet OutputFcn OutputSel Refine

Stats

A felhasználó által beállítható kimeneti függvény (ilyen lehet pl. az odeplot függvény). A megoldásmátrix mely indexeit adja át a kimeneti fv.-nek. Ha a megoldó túl nagyokat lép és nem vagyunk megelégedve a megoldás felbontásával, ezzel átállíthatjuk az elmentett lépések s˝ur˝uségét, így a program maga nem lassul. Statisztika

Jacobi mátrix Jacobian JPattern

Vectorized

A Jacobi mátrixot explicit módon, függvény alakjában a felhasználó definiálja. A Jacobi mátrix nemnulla elemei helyére 1-et, máshova 0-t kell írni. így a megoldó a nulla elemeknek megfelel˝u deriváltakat meg sem próbálja kiszámítani, ami gyorsítja a programot. Vektorizálás. (Én nem igazán értem, miért jó. Ha valaki elmesélné, hálás lennék.)

Lépésköz InitialStep MaxStep

Javasolt kezdo˝ lépés. Maximális lépésköz. Ha a kiszámított megoldást túl durvának tartjuk és s˝uríteni szeretnénk a pontokat, a refine parancsot használjuk!

30


31

Esemény figyelés Ha valamilyen esemény bekövetkezését szeretnénk nyomon követni (pl. ütközés, valamilyen határérték elérése, stb.) ezt kapcsoljuk be. Ilyenkor az integrálás megállítható.

Events

9.1.1. Szabadsugár sebességeloszlása Egy kétdimenziós szabadsugár sebességeloszlásának számítása az alábbi kezdeti érték feladatra vezet: F F˙ = F¨ 2 ,

ahol F(0) = 1 e´ s

˙ F(0) = 0.

Bevezetve az y1 = F és az y2 = F˙ új koordinátákat kapjuk az els˝orend˝u DE rendszert: ! ! y˙ 1 y2 = √ y˙ 2 − y 1 y2 (A negatív el˝ojel˝u gyök kell y˙ 2 egyenletében, különben a megoldások a végtelenhez tartanak.) Az alábbi kis program megoldja a feladatot (szabadsugar.m) . Szabadsugar 1.2

1

0.8

F

0.6

0.4

0.2

0

−0.2

0

0.5

1

ξ

1.5

2

9.1. ábra. A szabadsugár feladat megoldása.

szadadsugar.m function szabadsugar [t,y] = ode45(@f,[0 2.5],[0; 1]); plot(t,y(:,2)); xlabel(’\xi’); ylabel(’F’); title(’Szabadsugar’); function dydt = f(t,y,mu) dydt = [ y(2) -sqrt(y(1)*y(2))];

2.5


32

9.1.2. Egy merev feladat: a Van der Pol egyenlet Balthazar van der Pool 1927-ben egy vákuumcsövet tartalmazó egyszer˝u áramkör vizsgálata közben id o˝ nként zajokat figyelt meg. Az áramkört modellezo˝ egyenlet x¨ + µ x2 − 1 x˙ + x = 0 valóban öngerjesztett rezgéseket mutat bizonyos µ paramétertartományokban (ami meglep o˝ volt, hiszen explicit id˝ofüggés nincs az egyenletekben, azaz senki sem ’rázza’ a rendszert). Hogy a numerikus megoldót használni tudjuk, el˝oször Cauchy átírással els˝orend˝u rendszerré alakítjuk az y1 = x és y2 = x˙ új változók segítségével, amivel az egyenlet az ! ! y2 y˙ 1 = y˙ 2 µ 1 − y21 y2 − y1 alakot ölti. A (van_der_Pol.m) példaprogramban átállítjuk a relatív hibát 10 −4 -re és az ode15s megoldót használjuk, mert a feladat merev (gyors-lassú mozgást végez a rendszer). van_der_Pol.m function van_der_Pol(mu) options = odeset(’RelTol’,1e-4); [t,y] = ode15s(@vdp,[0 3000],[2; 0],options,mu); plot(t,y(:,1)); xlabel(’t’); ylabel(’y_1’); title([’van der Pol egyenlet, \mu=’,num2str(mu)]); function dydt = vdp(t,y,mu) dydt = [ y(2) mu*(1-y(1)^2)*y(2)-y(1)];

van der Pol egyenlet, µ=1000 2.5

2

1.5

1

y1

0.5

0

−0.5

−1

−1.5

−2

−2.5

0

500

1000

1500 t

2000

2500

9.2. ábra. A van der Pol egyenlet megoldása.

3000


33

9.1.3. Eseménykezelés: vízcsepp alakja A csapból lassan kifolyó vízcsepp y(x) alakját a y00 1+y

+ 02 3/2

y0 x 1+

y02 1/2

= 2y00 (0) −

ρg y σ

egyenlet írja le, ahol σ a felületi feszültség, ρ a folyadék s˝ur˝usége. Az egyenletet numerikus célokra átírva kapjuk a !   y2 y˙ 1 3/2 =  y˙ 2 (A − By1 ) 1 + y22 −

y2 t

1 + y22

  

A (csepp.m) programban a numerikus megoldás során figyeljük y 2 értékét (ami a felületet leíró görbe meredeksége) és ha már majdnem függ˝oleges (itt y2 > 104 ), abbahagyjuk az integrálást. Ezt valósítja meg az events(t,x) eljárás. csepp.m function csepp global A B A=1; B=0.1343; options = odeset(’Events’,@events); [t,x] = ode45(@f_csepp,[1e-5 10],[0; 0],options); plot(t,x(:,1),’k’,-t,x(:,1),’k’); title([’Vizcsepp alak, A=’,num2str(A)]); xlabel(’x’); ylabel(’y’); %----------------------------------------------------function dxdt = f_csepp(t,x) global A B dxdt = [x(2); (A-B*x(1))*(1+x(2)^2)^(3/2)-x(2)/t*(1+x(2)^2)]; %-----------------------------------------------------function [value,isterminal,direction] = events(t,x) value = abs(x(2))-10^4; isterminal = 1; % megall az integralas direction = 0; % mindket iranyban

9.2. Perem érték feladatok (BVP) A peremérték feladatok esetén egyrészro˝ l ugyanúgy definiálni kell a DE jobb oldalát, mint a kezdeti érték feladatoknál, de két további eljárásra is szükség van: a peremfeltételeket definiáló függvényre és a megoldót inicializáló eljárásra. Nagyon fontos, hogy ’jó’ kiindulófüggvényt válasszunk az iterációhoz, mert ez jelent o˝ sen megkönnyíti a megoldást. BVP megoldó szintaktika: sol = bvp4c(odefun,bcfun,solinit),


34

Vízcsepp alak, A=1 3

2.5

y

2

1.5

1

0.5

0 −2.5

−2

−1.5

−1

−0.5

0 x

0.5

1

1.5

2

2.5

9.3. ábra. A vízcsepp alakja.

ahol odefun tartalmazza a DE jobb oldalát. Formája a szokásos dydx = odefun(x,y), ahol x skalár, y és dydx pedig oszlopvektor. bcfun definiálja a peremfeltételeket nullára rendezve. Formája res = bcfun(ya,yb), ahol ya és yb oszlopvektorok, melyek a megoldást reprezentálják a tartomány elején és végén. res egy oszlopvektor, mely a rezíduumokat tartalmazza. solinit inicializálja a megoldást x által kifeszített hálón y értékekkel. A peremfeltételeket a=solinit.x(1) és b=solinit.x(end) helyeken értelmezzük. Ez az eljárás opcionális, az inicializáló vektorok feltöltése másként is lehetséges (pl. adatfájlból). sol az eredményeket tartalmazó struktúra. Mivel a felhasznált numerikus megoldó mozgóhálós kollokációs módszer, az x-ben definiált kezdeti hálót a megoldó felülírta. Érdemes még a bvpinit és a deval eljárásokat használnunk; az els o˝ az inicializálást könnyíti meg, a második pedig a kapott eredmények kiértékelését. Fontos még tudnunk, hogy a peremfeltételek száma: N perem f eltételek = NODE + N paraméterek , ui. egyes esetekben szabad paramétereket kell hagyni a DE rendszerben a megoldás létezésének biztosítására. A következ˝o példákban extra paraméterek keresésével nem foglalkozunk (ld. Help: Matthieu egyenlet sajátértéke), de pl. ilyenre vezet periodikus megoldások keresése is, de ekkor a peremfeltétel egy ismeretlen x = T helyen van. Ilyenkor a T ismeretlen periódussal átskálázzuk az x változót és így már jól definiált a feladat: y0 (x) = T F(x, y) ahol y(0) = y(1), és még egy peremfeltételt el˝o kell írni T szabad paraméter miatt, ami általában a fázisra vonatkozik.

9.2.1. Couette áramlás sebességgradienssel Két síklap közötti lamináris, stacioner áramlást leírja az


35

∂2 v ∂p = µ 2, ∂x ∂y egyenlet, ahol a x a hosszirányú, y a keresztirányú koordináta. Ha a ∂p ∂x nyomásgradiens konstans, a fenti egyenlet zárt alakban megoldható. Két speciális eset a Poiseuille-áramlás (mindkét lap áll és a nyomásgradiens nemzérus, ekkor a sebességeloszlás parabolikus) és a Couette-áramlás (ekkor az egyik lap mozog, nincs nyomásgradiens, a sebességprofil pedig lineáris). Ha van nyomásgradiens és a fels o˝ lap is mozog, akkor is megoldható a feladat zárt alakban, a sebességprofil egy ’ferde tengely˝u’ parabola. Ha a két lap távolsága h = 1m, az mozgó fal sebessége v0 = 2m/s és az áramló közeg víz, akkor mivel Re = vh/ν ≈ 103 , az áramlás még a lamináris tartományban van. A probléma peremértékfeladatként megfogalmazva: y01 y02

= y2 1 ∂p , = ρ ν ∂x

a tapadás törvénye miatt a peremfeltételek pedig az y(0) = 0 és y(1) = v 0 alakot öltik. A Matlab program így néz ki : couette.m function couette global P v0 nu P

= -0.01; v0 = 2; nu = 0.001;

% Megoldas inicializalasa: solinit=bvpinit(linspace(0,1,10),@mat4init); % Megoldo hivasa: sol=bvp4c(@mat4ode,@mat4bc,solinit); % Eredmenyek interpolalasa: xint = linspace(0,1,100); Sxint = deval(sol,xint); % Szinezes csusztatofeszultseg szerint: tau_max=max(Sxint(2,:)); tau_min=min(Sxint(2,:)); for i=1:length(Sxint(1,:)) szin=(Sxint(2,i)-tau_min)/(tau_max-tau_min); plot(Sxint(1,i),xint(i),’+’,’MarkerEdgeColor’,[szin 0.5*szin 0.75-0.5*szin]); hold on; end; hold off; title([’Adatok: dp/dx=’,num2str(P),’ \nu=’,num2str(nu)],’FontSize’,14); xlabel(’aramlasi sebesseg, v’,’FontSize’,14); ylabel(’keresztemtszet, x’,’FontSize’,14); %--------------------------------------------function dydx=mat4ode(x,y) global P nu dydx=[ y(2); P/nu];


%---------------------------------------------function res=mat4bc(ya,yb) global v0 res=[ ya(1); yb(1)-v0]; %--------------------------------------------function yinit=mat4init(x) global v0 yinit = [ v0*x^2; 2*v0*x ];

Adatok: dp/dx=−0.01 ν=0.001 1

0.9

0.8

keresztemtszet, x

0.7

0.6

0.5

0.4

0.3

0.2

0.1

0

0

0.5

1

1.5

2

2.5

aramlasi sebesseg, v

9.4. ábra. Sebességeloszlás Couette áramlás esetén nyomásgradienssel, τ csúsztatófeszültség szerint színezve.

36

10. fejezet

Fourier analízis Fourier analízist és frekvenciaspektrumot az fft(y,N) paranccsal végezhetünk, ahol N 2 egész számú hatványa (mivel a Matlab gyors Fourier transzformációt alkalmaz). Inverz transzormáció az ifft paranccsal számítható. Például vegyünk egy 1000Hz-el mintavételezett jelet, mely egy 50 és egy 120 Hz-es komponenset tartalmaz, de el van rontva egy véletlenszer˝u, zérusátlagú zajjal. Végezzük el a Fourier-transzformációt és számítsuk ki a komplex amplitúdó abszolúrtékét (fft_demo.m ). fft_demo.m t = 0:0.001:0.6; x = sin(2*pi*50*t)+sin(2*pi*120*t); y = x + 2*randn(size(t)); n=8; Y = fft(y,2^n); Pyy = Y.* conj(Y) / 2^n; f = 1000*(0:2^(n-1))/2^n; hossz=min(length(Pyy),length(f)); subplot(2,1,1) plot(t,x), xlabel(’t’), ylabel(’x(t)’) grid on, axis([0 0.6 -2.5 2.5]) subplot(2,1,2) plot(f(1:hossz),Pyy(1:hossz)) xlabel(’f [Hz]’), ylabel(’amplitudo’) grid on

37

10. FEJEZET. FOURIER ANALÍZIS

38

2

x(t)

1 0 −1 −2 0

0.1

0.2

0.3 t

0.4

0.5

0.6

60

amplitudo

50 40 30 20 10 0

0

50

100

150

200

250 f [Hz]

300

350

10.1. ábra. Eredeti id˝ojel és a spektrum.

400

450

500

11. fejezet

Példaprogramok 11.1. Karakterisztikák módszere Hidraulikus tranziensek számítása esetén a csövekben az áramlást az 1D mozgásegyenlet és a kontinuitási egyenlet ∂v 1 ∂p + =S ∂t ρ ∂x

és

∂ρ ∂ρv + = 0. ∂t ∂x

Tegyük fel, hogy a s˝ur˝uség csak a nyomás függvénye, azaz ρ = ρ(p), így ∂ρ/∂t = dρ/dp × ∂p/∂t. Vezessük be a hangsebességet: dρ/dp = a−2 . Adjuk hozzá a mozgásegyenlethez a kontinitási egyenlet λ szorosát: ! ∂v ∂v λ ∂p 1 ∂p + + λρ + = S. ∂t a2 ∂t ∂x λρ2 ∂x Amennyiben az ismeretlen λ szorzót ± ρa -nak választjuk, az (x, t) síkon két speciális irányban ( dx dt = ±a) egy - egy közönséges differenciálegyenletet kapunk: dz d (p ± ρav) = ∓ρag ∓ ρaS. dt dx A numerikus módszer ezek után egyszer˝u. A fenti differenciálegyenletet egy egyenköz˝u rácson oldjuk meg, ahol az aktuális i-edik pontban a jellemz˝oket az i−1 b˝ol indított a meredekség˝u és az i+1-edik pontból indított −a meredekség˝u karakterisztika segítségével határozzuk meg. Két programot készítettünk: egy önálló megoldót kar.m és egy vezérl o˝ programot: kar_driver.m. A megoldó bemenete a cs˝o egyes paraméterei, a ’régi’ p és v nyomáseloszlás, ill. a peremfeltételek típusa és értéke (2 db). Ezek alapján kiszámítja és visszaadja a ∆t-vel késo˝ bbi cs˝oállapotot. A vezérl˝oprogram ezt az eljárást ismétli, ill. elo˝ készíti a számítást és megjeleníti az eredményeket.

39

11. FEJEZET. PÉLDAPROGRAMOK

40

kar.m function [puj,vuj]=kar(adatok,p,v,PF0_tipus,PF0_ertek,PFN_tipus,PFN_ertek) D L N a ro

= = = = =

adatok{1}(1); adatok{1}(2); adatok{1}(3); adatok{1}(4); adatok{1}(5);

% % % % %

csoatmero csohossz csodarabok szama (csomopontok szama: N+1) hangsebesseg suruseg

g = 9.81; dx = L/N; dt = dx/a; roa = ro*a; dzdx = adatok{2}; nu = 1e-6; lambda=0.02; for i=2:N ar = p(i-1)+ro*a*v(i-1); br = p(i+1)-ro*a*v(i+1); vuj(i) = (-2*dt*roa*g*dzdx(i)+ar-br)/roa/(2+dt*lambda/2/D*(abs(v(i+1))+abs(v(i-1)))); puj(i) = ( -dt*roa*g*dzdx(i)+ar)-roa*(1+dt*lambda/2/D* abs(v(i-1)))*vuj(i); end % Peremfeltetelek beallitasa. perem_1 = p(2)-ro*a*v(2) + dt*roa*g*dzdx(1); perem_N = p(N)+ro*a*v(N) - dt*roa*g*dzdx(N+1); szorzo1 = roa*(1+dt*lambda/2/D*abs(v(2))); szorzoN = roa*(1+dt*lambda/2/D*abs(v(N))); switch PF0_tipus case ’p’ puj(1) = PF0_ertek; vuj(1) = (puj(1)-perem_1)/szorzo1; case ’v’ vuj(1) = PF0_ertek; puj(1) = perem_1 + ro*a*vuj(1); otherwise error(’Ismeretlen peremfeltétel!’); end switch PFN_tipus case ’p’ puj(N+1) = PFN_ertek; vuj(N+1) = -(puj(N+1)-perem_N)/ro/a; case ’v’ vuj(N+1) = PFN_ertek; puj(N+1) = perem_N-szorzoN*vuj(N+1); otherwise error(’Ismeretlen peremfeltétel!’); end


41

kar_driver.m function kar_driver % adatok D=0.05; A=D^2*pi/4; L=10; N=10; a=1000; ro=1000; lambda=0.02; tmax=0.1; % inicializalas t=0; dt=L/N/a; x=0:L/N:L; dzdx=zeros(1:N); % cell{} tip. vektorral több adattipust be lehet csomagolni adatok{1}=[D L N a ro]; adatok{2}=dzdx; % Inicializálás (t=0-hoz tartozó állapot) p=2e5*ones(1,N+1); v=zeros(1,N+1); i=1; % indulamandula while t
5 4.5 4

p [bar]

3.5 3 10

2.5 8 2 6 1.5 4 1 0.1

2 0.09

0.08

0.07

0.06

0.05

0.04

0.03

0.02

0.01

0 0

t [s]

11.1. ábra. A kar_driver.m program grafikus kimenete.

x [m]


42

11.2. Nyíltfelszínu˝ csatornaáramlás Nyíltfelszín˝u stacioner csatornaáramlás esetén y+z+

v2 + h0 = a´ llandó, 2g

ahol y(x) a folyadékfelszín alakja, z a csatorna fenekének magassága, v az áramlási sebesség, h 0 pedig a veszteségmagasság. Ez utóbbi számítása a Chézy képlettel történik: dh0 v2 = 2 , dx C Rh

R1/6 e´ s C = h . n

A ahol Rh = K

Itt A(y) a nedvesített felület, K(y) a nedvesített kerület, n pedig az cso˝ érdességét˝ol függ˝o tényez˝o. Figyelembe véve a kontinuitást, azaz hogy Q = Av = a´ llando, ´ a fenti képlet kicsit megdolgozva a 2

Q dy i − A2 C2 Rh = 2 dx 1 − QA3 gB

alakot ölti, ahol i a csatorna lejtése, B a folyadéktükör szélessége és g a gravitációs gyorsulás. A feladatot kezedti érték feladatként fogalmazzuk meg, azaz el o˝ írunk egy kezdeti folyadékszintet és térfogatáramot. Ismert, hogy ha a térfogatáram kisebb, mint az egyenes felszínhez tartozó normál térfogatáram, a folyadékszint duzzad, ha pedig nagyobb, akkor gyorsul a mozgás és kialakulhat a rohanás jelensége. Ezt kezelend o˝ eseményfigyelést alkalmazunk, azaz a program automatikusan abbahagyja az integrálást, ha bekövetkezik az y = 0 esemény (rohanas(x,y)). Az alábbi program (nyiltfelszin.m) az egyszer˝uség kedvéért b szélesség˝u, téglalap keresztmetszet˝u csatornát feltételez . Futtassuk le a programot Q = 0.2, 0.25, 0.298, 0.299, 0.3[m 3/s] térfogatáramokkal!

folyadekfelszin, y [m]

Nyiltfelszinû csatornaáramlás: i=0.001 Q=0.298 csatorna feneke folyadék felszíne 1.5

1

0.5

aramlasi sebesseg, v [m/s]

0

0

100

200

300

400

500

600

700

800

900

1000

0

100

200

300

400

500

600

700

800

900

1000

1 0.98 0.96 0.94 0.92 0.9 0.88

csohossz, x [m]

11.2. ábra. A nyiltfelszin.m program grafikus kimenete.


43

nyiltfelszin.m function nyiltfelszin global B n g Q ii Q = 0.298; y0= 0.6; ii = 0.001; n = 0.01; B = 0.5; g = 9.81; L = 1000;

% % % % % % %

eloirt terfogataram szint a csatorna elején hidraulikai lejtes csoerdesseg csatorna szelesseg gravitacio csohossz

Rh = y0*B/(2*y0+B); C = Rh^(1/6)/n; Qn = y0*B*C*sqrt(ii*Rh); fprintf(’\n Nyiltfelszinu csatornaaramlas\n’); fprintf(’\n %g m-es csőeleji felszinhez tartozó térfogatáram: %g m^3/s’,y0,Qn); % Megoldo hivasa: options=odeset(’events’,@rohanas); [x,y]=ode45(@nyiltfelszin_de,[0 L],y0,options); % A csatorna alja és az ahhoz képesti folydékszint (rajzoláshoz): for i=1:length(x) yy(i)=ii*(1-x(i)/L)*L; yf(i)=yy(i)+y(i); v(i) =Q/y(i)/B; end figure(1) subplot(2,1,1), plot(x,yy,x,yf); title([’Nyiltfelszinű csatornaáramlás: i=’,num2str(ii),’ Q=’,num2str(Q)],’FontSize’,12); ylabel(’folyadekfelszin, y [m]’,’FontSize’,12) legend(’csatorna feneke’,’folyadék felszíne’) axis([0 L 0 yf(1)*1.2]) subplot(2,1,2), plot(x,v); xlabel(’csohossz, x [m]’,’FontSize’,12); ylabel(’aramlasi sebesseg, v [m/s]’,’FontSize’,12); %--------------------------------------------function dydx=nyiltfelszin_de(x,y) global B n g Q ii A=y(1)*B; Rh=y(1)*B/(2*y(1)+B); C=Rh^(1/6)/n; dydx=[ (ii-Q^2/A^2/C^2/Rh)/(1-Q^2*B/A^3/g)]; %--------------------------------------------function [val,ter,dir]=rohanas(x,y) val=y(1); % Az esemény definiciója ter=1; % Állítsuk meg az integrálást. dir=0; % Mindkét irányban (1: + -> -, 2: - -> +)


44

11.3. Függ˝olegesb˝ol vízszintesbe vezet˝o ívben mozgó anyagdugó Pneumatikus anyagszállítás témakörében felmerülo˝ probléma, hogy egy függo˝ legesb˝ol vízszintesbe vezet˝o 90o -os ívben hogyan mozog az anyagdugó. Tegyük fel, hogy az anyagdugó hossza nagyobb, mint az ív, azaz L > R π/2. Ekkor a mozgás 3 szakaszra osztható: • az anyagdugó tölti az ívet, • az anyagdugó az ívben mozog és • az anyagdugó ürül az ívb˝ol. A független változók az anyagdugó elejének (ill. végének) α szögelfordulása és az anyagdugó sebessége v a ,. helyzete z és gyorsulása aa . A mozgásegyenletek az egyes szakaszokban: 1. szakasz (telít˝odés): 0 ≤ α ≤  0   α =    z0 =

π 2

z F1 (α, z) = C31 + Π0 Π2

q

k dhs ρg2 p2

2 vg2 − Rz + p22 −

µ 2

(cos α − α) −

sin α−α 2

− αµRz2

2. szakasz (mozgás az ívben): z + R π2 ≤ L:  0  v =    a0 aa =     z0 = a

3. szakasz (ürülés): 0 ≤ α ≤

π 2

( 0 α = z0 =

aa F2 (aa , va ) va

z F3 (α, z)

Nincs értelme itt pontosan definiálni az egyes változókat és fizikai konstansokat, ezt az el o˝ adásjegyzete alapján mindenki megteheti. (Egyébként F 1 a legegyszer˝ubb mozgásegyenlet, ezért ezt írtam ki.) Azt viszont vegyük észre, hogy az egyenleteket egymás után szakaszosan kell megoldani éspedig úgy, hogy eseménykezeléssel meg kell állítani az integrálást és a végállapotból indítani a következo˝ szakasz megoldását. Ez az 1. és a 3. szakaszban viszonylag egyszer˝u, mert α értékét kell figyelni, ami az egyik független változó. Ellenben a 2. szakaszban az elmozdulást kell figyelni, amire (eredetileg) nincs egyenlet (csak va -ra és aa -ra), ezért csatolunk egy harmadik diff. egyenletet, mely az egyszer˝uen R a z = va kapcsolatból nyerhet˝o. Így már tudjuk követni az anyagdugó helyzetét is. A példaprogram megtalálható az interneten www.hds.bme.hu. A kimenet a 11.3 ábrán látható, az egyes szakaszok különböz o˝ szinekkel vannak feltüntetve.


45

s [m]

4

2

0

0

1

2

3

4

5

6

7

0

1

2

3

4

5

6

7

0

1

2

3

4

5

6

7

0

1

2

3

4

5

6

7

v [m/s]

3

2.5

2

a [m/s2]

2 0 −2 −4

Dp [Pa]

10.5 10 9.5 9

t [s]

11.3. ábra. A pneu_iv.m program grafikus kimenete.

Matlab Hogyan. c GPL. BME Hidrodinamikai Rendszerek Tanszék február 17. belső használatra

Recommend Documents