Physics | Hydrodynamics » Dr. Gausz Tamás - Néhány numerikus módszer mérnöki alkalmazásokhoz, áramlástani példa-feladattal

Dr. Gausz Tamás Néhány numerikus módszer mérnöki alkalmazásokhoz áramlástani példa-feladattal Budapest, 2020 Tartalomjegyzék Bevezetés 1 I. Bevezető ismeretek I.1 Relációk I.2 Direkt szorzat, szorzathalmaz I.3 Lineáris tér I.4 Lineáris altér, a tér bázisa I.5 Operátor, leképezés, funkcionál I.6 Normált tér, metrika I.7 Sorozatok, konvergencia, határérték, teljes tér I.8 Skaláris vagy belső szorzat I.9 Ortogonalitás I.10 Minimalizáló vektor, minimalizáló elem 2 2 3 4 6 8 9 11 12 13 14 II. Az alap-feladat és néhány megoldási módszer II.1 Az alap-feladat 15 15 II.2 A megoldás lehetőségei 16 II.3 A megoldás-függvény közelítése 17 III. A minta-feladat III.1 A modell-folyadék III.2 A Hamilton féle nabla-operátor és a Laplace operátor III.3 A minta-feladat III.4 A zárt alakú megoldás vizsgálata 19 19 19 21 25 IV. A közelítő megoldás függvénytere és az eltérés (hiba) függvény IV.1 A közelítő

megoldás függvénytere IV.2 Az eltérés (hiba) függvény a minta-feladat esetében 30 30 32 V. A kollokációs módszer V.1 A kollokációs módszer leírása V.2 A kollokációs módszer alkalmazása V.3 Kollokáció egyetlen bázis-függvénnyel, pontválasztás 33 33 33 36 VI. A legkisebb négyzetek módszere VI.1 A legkisebb négyzetek módszerének leírása VI.2 A legkisebb négyzetek módszerének alkalmazása VII. A Galjorkin eljárás VII.1 A Galjorkin eljárás leírása VII.2 A Galjorkin eljárás alkalmazása 39 39 39 43 43 43 Irodalomjegyzék 45 Mellékletek I. Kollokációs módszer, program-lista II. Legkisebb négyzetek módszere, program-lista 46 51 Bevezetés Az itt bemutatandó ismeretek bepillantást tesznek lehetővé néhány – általunk fontosnak gondolt – numerikus módszer felépítésébe és alkalmazásuk lehetőségeibe. Az elvi bemutatást áramlástani példával illusztráljuk A módszerek, természetesen más, illeszkedő

tudományterületre is alkalmazhatók. Feltétlenül rögzíteni kell azonban, hogy e tudományterületen, különösen az utóbbi évtizedekben, hatalmas mértékű fejlődés tapasztalható, és az itt bemutatott ismeretek e terület csak nagyon kis részét fedik le. E munka fő célja néhány, fontos elv bemutatása, a konkrét feladatmegoldásra a professzionális CFD, illetve multi-fizika programokat ajánlhatjuk! Hálás köszönetemet fejezem ki Prof. Dr Zobory István egyetemi tanár úrnak, aki a tanácsaival és megjegyzéseivel nagymértékben hozzájárult e munka színvonalának emeléséhez. A munka eredetileg előadást támogató segédanyagnak (bemutató) készült, ezért önmagában nem feltétlenül követhető – adott esetben a megfelelő szakirodalomban (ennek a tudományterületnek a szakirodalma rendkívül kiterjedt) történő utánolvasásra lehet szükség! Illetve az anyag mélyebb elsajátításához feltétlenül szükséges a saját munka. A

bemutató példákban, egyes helyeken használjuk a „SCILAB” matematikai programot (www.scilaborg) – ez egy szabad szoftver, egyébként meglehetően hasonlít a „MATLAB” kereskedelemben megvásárolható szoftverhez Használjuk még a „Yabasic” programozási nyelvet (www.yabasicde) – ez a SCILAB-hoz hasonlóan szabad szoftver. A jelöléseket az első előfordulások helyén definiáljuk. Itt hívjuk fel a figyelmet arra, hogy az x az y és a z általános halmaz-elemet jelent, az x, y, z -vel viszont a Descartes-féle koordinátarendszer három koordinátáját jelöljük Az anyagban viszonylag sok összefüggés, képlet fordul elő, ezek közül soknak nem szerepel a részletes levezetése – ezért is meg azért, mert hiba vagy elírás – gondos munka ellenére is – bárhol előfordulhat, a Tisztelt Érdeklődőnek alkalmazás előtt minden képletet önállóan ellenőriznie kell! 1 I. Vizsgáljuk meg a munka-környezetünket I. Bevezető ismeretek I.1

Relációk A relációkkal foglakozó, igen kiterjedt tudományterületről itt csak nagyon röviden szólunk: a véleményünk szerint legfontosabb relációkat mutatjuk csak be. Az A és a B halmazok közötti R binér (kétváltozós) relációt (meghatározott kapcsolattal) az [ a, b] ( a ∈ A, b ∈ B ) rendezett párnak valamely, meghatározott tulajdonsága jellemez. Bármely rendezett párra meg kell tudni állapítani, hogy rendelkezik-e a kérdéses, meghatározott tulajdonsággal, vagy sem. Azaz érvényes-e az aRb vagy sem. Elsőnek tekintsük az egyenlőséget (ekkor az aRb azt jelenti, hogy a = b ). Ez lineáris, binér, ekvivalencia reláció, azaz: ï Reflexív, vagyis teljesül az a = a ; ï Szimmetrikus, vagyis teljesül az a = b és a b = a is; ï Tranzitív, vagyis ha a = b és b = c , akkor a = c . Tekintsük másodszorra a kisebb vagy a nagyobb relációt (ekkor az aRb azt jelenti, hogy a < b vagy a > b ). Mindkét binér, lineáris reláció szigorú

rendezés, azaz: ï Irreflexív, vagyis az a < a nem teljesül és az a > a nem teljesül; ï Aszimmetrikus, vagyis ha a < b , akkor a b < a nem teljesül, és ha a > b , akkor a b > a nem teljesül; ï Tranzitív, vagyis ha a < b és b < c , akkor a < c , és ha a > b és b > c , akkor a > c . Tekintsük harmadszorra a kisebb-egyenlő vagy a nagyobb-egyenlő relációt (ekkor az aRb azt jelenti, hogy a ≤ b vagy a ≥ b ). Mindkét binér, lineáris reláció rendezés, azaz: ï Reflexív, vagyis az a ≤ a teljesül és az a ≥ a is teljesül; ï Antiszimmetrikus, vagyis ha a ≤ b , és b ≤ a , akkor a = b, és ha a ≥ b , és b ≥ a , akkor b = a ; ï Tranzitív, vagyis ha a ≤ b és b ≤ c , akkor a ≤ c , és ha a ≥ b és b ≥ c , akkor a ≥ c . 2 I. Vizsgáljuk meg a munka-környezetünket I.2 Direkt szorzat, szorzathalmaz Legyen A és B két, tetszőleges halmaz. Tekintsük minden olyan {a, b} párt, amelynél a ∈ A és b

∈ B . Ezt, az új halmazt az A és B halmaz szorzathalmazának nevezzük; jelölése A × B A szorzathalmazt előállító műveletet pedig direkt szorzatnak nevezzük Legyenek A1 , A2 Ai An tetszőleges halmazok és a1 ∈ A1 , a2 ∈ A2 ai ∈ Ai an ∈ An , akkor a szorzathalmazuk: A1 × A2 × × Ai × × An = {a1 , a2 ai an } . A szorzathalmazokat mi alapvetően a valós számok halmazának felhasználásával képezzük. Legyen ℜ a valós számok halmaza, ezek a számegyenesen helyet foglaló számok. Az ℜ× ℜ = ℜ2 például egy végtelen síkon helyet foglaló, összes ponttal szemléltethető. Az ℜ× ℜ× ℜ = ℜ3 halmaz pedig egy, végtelen, háromdimenziós tér összes pontjával szemléltethető n Végül az ℜ× ℜ× × ℜ = ∏ ℜ =ℜ n egy, n-dimenziós végtelen tér pontjaival i=1 szemléltethető. 3 I. Vizsgáljuk meg a munka-környezetünket I.3 Lineáris tér Az X nemüres halmazt lineáris térnek nevezzük a Φ számtest felett

(valós: Φ = ℜ vagy komplex: Φ = C ) ha: I. ∀ x, y ∈ X -hez hozzá van rendelve ( x + y ) ∈ X , amelyet az x és y elemek összegének nevezünk, és az összeadás művelete négy azonosságot teljesít: 1. x + y = y + x ∀ x, y ∈ X -re; 2. ( x + y ) + z = x + ( y + z ) ∀ x, y , z ∈ X -re; 3. ∃ ∅ : x + ∅ = x ∀ x ∈ X -re, vagyis létezik zérus elem; 4. ∀ x ∈ X -hez ha ∃ y ∈ X , hogy az x + y = ∅ , akkor létezik additív inverz, és y az x inverze: y = − x . II. ∀λ ∈ Φ és ∀ x ∈ X párhoz hozzá van rendelve ( λ x ) ∈ X , amelyet az x nek a λ számmal (skalárral) való szorzatának nevezünk, és a szorzásra teljesül az alábbi, négy tulajdonság: 1. (α β ) x = α ( β x ) ∀ x ∈ X -re és ∀α , β ∈ Φ -re; 2. α ( x + y ) = α x + α y ∀ x, y ∈ X -re és ∀α ∈ Φ -re; 3. (α + β ) x = α x + β x ∀ x ∈ X -re és ∀α , β ∈ Φ -re; 4. létezik egység elem: ∃ 1∈ Φ , hogy 1 x = x , ∀ x ∈ X

-re Megjegyzések: mi a következőkben a Φ = ℜ valós számtesttel dolgozunk; a valós számtesten értelmezett összeadás és szorzás kommutatív (felcserélhető), azaz ∀α , β ∈ℜ : α + β = β + α és αβ = βα <a kivonás és az osztás nem kommutatív>; asszociatív (csoportosítható), azaz ∀α , β , γ ∈ℜ : α + β + γ = α + ( β + γ ) = (α + β ) + γ és αβγ = α ( βγ ) = (αβ ) γ <a kivonás és az osztás nem asszociatív>; disztributív (tagolható) a valós számok szorzása az összeadásra nézve disztributív, azaz: ∀α , β , γ ∈ ℜ : α ( β + γ ) = α β + α γ ; 4 I. Vizsgáljuk meg a munka-környezetünket Példák: Tekintsük a valós szám n-eseket. Legyen e halmaz két, tetszőleges eleme: x = ( x1 , x2 xn ) ∈ℜn és y = ( y1 , y2 yn ) ∈ℜn Legyen az öszszeadás a következő: x + y = ( x1 + y1 , x2 + y2 xn + yn ) ∈ℜn , illetve a számmal történő szorzás α x = (α x1 ,α x2 α

xn ) ∈ℜn , α ∈ Φ . Ez, nyilvánvalóan lineáris tér A fenti példa, az n = 3 esetben éppen a hagyományos, háromdimenziós vektorteret is jelentheti, azzal a megjegyzéssel, hogy ezekben a vektorterekben a számmal történő szorzás mellett skaláris vagy belső, vektoriális és diadikus vagy külső szorzatot is értelmezünk. Tekintsük az [ a, b] zárt intervallum felett folytonos függvények halmazát: f , g ∈ C [ a, b ] . Értelmezzük az összeadást ( f + g )( x ) = f ( x ) + g ( x ) , x ∈ [ a, b ] és a számmal történő szorzást (α f )( x ) = α f ( x ) , α ∈ Φ, x ∈ [ a, b ] , f ∈ C [ a, b ] . Ez is lineáris tér Ellenpéldaként tekintsük a C [ a, b ] halmazból a pozitív értékeket felvevő függvények alkotta részhalmazt. Ez nyilvánvalóan nem alkot lineáris teret, hiszen például a mínusz eggyel történő szorzás kivezet ebből a térből. 5 I. Vizsgáljuk meg a munka-környezetünket I.4 Lineáris altér, a tér

bázisa Egy X lineáris térnek a nemüres M ( M ⊆ X ) részhalmaza lineáris altér, ha maga is lineáris tér az X térre definiált összeadással és skalárral való szorzásra nézve, vagyis, ha: ∀ x, y ∈ M ⇒ ( x + y ) ∈ M ; és ∀ x ∈ M és ∀α ∈ Φ ⇒ Példa: (α x ) ∈ M . A szokásos geometriai tér, mint vektortér lineáris altere lehet minden, az origóhoz illeszkedő síkoz rendelt vektortér. Legyen xi ∈ X és α i ∈ ℜ, i = 1, 2, n ; az xi elemek α i együtthatókkal képzett n lineáris kombinációja: ∑α i =1 i xi . Az xi elemeket lineárisan függetlennek nevezzük, ha velük a n ∑α i =1 i xi = ∅ csak úgy lehetséges, ha α i = 0, ∀i -re. Mivel egy lineáris tér altere maga is lineáris tér, a következőkben az „altér” megnevezést csak akkor írjuk ki, ha az szükséges. Ha egy lineáris térben tetszőlegesen választott n számú elem lineárisan független, de bármely ( n + 1) számú elem már

lineárisan összefüggő, akkor a tér n -dimenziós. Az {xi } elem-rendszer generátorrendszer, ha az elemeinek lineáris kombi- nációjával bármely x ∈ X felírható. Ha pedig az {xi } elem-rendszer elemei lineárisan függetlenek, akkor {xi } bázis A bázist alkotó elemeket báziselemeknek – például vektortér esetében bázisvektoroknak, függvénytér esetében bázisfüggvényeknek – nevezzük Kimondjuk, hogy: véges dimenziós lineáris terekben, az {xi } bázis felhaszn nálásával ∀ x ∈ X előállítható az x = ∑ α i xi lineáris kombinációval, és ez az elői =1 állítás egyértelmű! 6 I. Vizsgáljuk meg a munka-környezetünket Példák: A három dimenziós (3D) geometriai teret, mint az origóból felrakott, irányított egyenesdarabok halmazát a három, jól ismert, speciális választás esetében egymásra merőleges egységvektor ( i, j, k ) feszíti ki, és minden v helyvektor előállítható a v = α i + β j + γ k alakban (

α , β és γ ∈ℜ ). Ebben az esetben természe- tesnek vesszük, amit a korábban kimondott tétel bizonyít, hogy ti. bármely vektor bázis-előállítása egyértelmű. Tekintsük másik példaként azt, az alkalmasan választott {ϕi } , i = 1, 2 n bázisfüggvény-rendszert, amely bázisfüggvényekkel a megfelelő, n -dimenziós függvénytér előállítható. Vegyük igen egyszerű, konkrét példaként a [ −1,1] zárt interval- lum felett értelmezett, Csebisev polinomok első két tagját: T0 ( t ) = 1 ; T1 ( t ) = t , t ∈ [ −1,1] . Evvel a két függvénnyel bármely, a szóban forgó intervallum feletti egyenes darab generálható (előállítható). Mondjuk konkrétan az nyilvánvalóan az 5 t + 7 egyenes-darab 5 T0 ( t ) + 7 T1 ( t ) módon generálható. 7 I. Vizsgáljuk meg a munka-környezetünket I.5 Operátor, leképezés, funkcionál Tekintsük az X lineáris teret és az Y (lineáris) teret, illetve a T műveletet (operátort, operációt).

Legyen D ⊆ X és R ⊆ Y A T művelet (operátor, operáció) ∀ x ∈ D ⊆ X elemhez egyértelműen hozzárendeli az y = T x ∈ R ⊆ Y elemet: T : D R . A D halmaz a T operátor (értelmezési) tartománya, az R = {T x ∈ R ⊆ Y } halmaz pedig az operátor képtere. A T operátor egyértelmű, ha egy x ∈ D -hez egy y ∈ R tartozik, azonban különböző x -ekhez tartozhatnak azonos y -ok. A T művelet (operátor, operáció) bijekció, ha a hozzárendelés kölcsönösen egyértelmű, azaz ∀ x -hez különböző y tartozik. A T : D R operátor lineáris, ha T (α x + β y ) = α T x + β T y , itt: tetszőleges α , β ∈ℜ és x, y ∈ D . Ez a T összeg- és aránytartó tulajdonságát posztulálja Példák: az y = A x mátrix-vektor szorzat homogén, lineáris leképezés; lineáris operátor a T f ( t ) = d f ( t ) differenciálás [ f ( t ) differendt ciálható függvény!]; t lineáris operátor a T f ( t ) = ∫ f (τ ) dτ , t ,τ ∈ [ a, b]

integrálás a [ f ( t ) [ a, b] -ben Riemann szerint integrálható függvény!] A normált, lineáris téren értelmezett T : D R operátor kontrakciós, ha T x − T y ≤ α x - y , 0 < α < 1 . Itt hívjuk fel a figyelmet arra, hogy az igen fontos kontraktív leképezésekkel – a kontrakciós operátorokkal – illetve az ezekkel kapcsolatos ismeretekkel és tételekkel bővebben például [7] 5. fejezete foglalkozik A T : D ℜ funkcionál, azaz olyan leképezés, amely a valós (vagy a komplex) számtestre képez le. Példák: a T x = x norma, nemlineáris funkcionál ((itt az operátor a normaképzés); b a T x = ∫ x ( t ) dt határozott integrál, lineáris funkcionál (itt az opea rátort a határozott integrál jelenti). 8 I. Vizsgáljuk meg a munka-környezetünket I.6 Normált tér, metrika Az X lineáris tér normált tulajdonságú, ha ∀ x ∈ X -hez hozzá van rendelve pontosan egy (egy és csak egy) x nem negatív, valós szám az x normája,

úgy, hogy a következő tulajdonságok teljesülnek: ∀ x ∈ X : ∃ x ∈ℜ+ , hogy: ℜ+ = [ 0, ∞ ) 1./ x ≥ 0 ; illetve ha x = 0 , ez akkor és csak akkor igaz, ha x = ∅ ; 2./ α x = α x α ∈ℜ ; 3./ x + y ≤ x + y x, y ∈ X ; (háromszög egyenlőtlenség). Példák: tekintsük az x = ( x1 , x2 xi xn ) , xi ∈ℜ, i = 1, 2, n szám n-eseket (ezek lineáris teret alkotnak): n 1./ abszolút-érték-összeg norma: x 1 = ∑ xi ; i =1 2./ Euklidész-i norma: x 2 = n ∑x i =1 2 i ; 3./ koordináta abszolút érték maximum norma: x 3 = max xi i = 1, 2 n tekintsük az [ a, b] zárt intervallum felett (négyzetesen is) integrálható, valósértékű függvények lineáris terét: b 4./ abszolút érték norma: x 1 = ∫ x ( t ) dt ; a b 5./ Euklidész-i norma: x 2 = + ∫ x (t ) 2 dt . a Legyen X nem üres halmaz, és d : X × X ℜ+ , ∀ x, y ∈ X elempárhoz nemnegatív értékeket rendelő függvény. Ha ∀ x, y , z ∈ X -re igaz, hogy:

1./ d ( x, y ) = d ( y, x ) ; 2./ d ( x, y ) = 0 akkor és csak akkor igaz, ha x = y ; 3./ d ( x, y ) + d ( y, z ) ≥ d ( x, z ) ; akkor d metrika vagy távolságfüggvény. Ekkor d segítségével két elem (pont) távolsága definiálható. Metrikus tér: egy halmaz, és egy, a halmazon értelmezett metrika (távolságfüggvény) együttese: ( X , d ) . A metrikát a norma indukálja, illetve ezen a módon mérhetjük a távolságot egy halmaz két eleme között. 9 I. Vizsgáljuk meg a munka-környezetünket Példák: legyen az X a szokványos, három-dimenziós euklideszi tér, és legyen P1 ( x1 , y1 , z1 ) ∈ X és P2 ( x2 , y2 , z2 ) ∈ X két, tetszőleges pont. E két pont távolsága a korábban „2” számmal jelölt, euklideszi norma felhasználásával: d ( P1 , P2 ) = n ∑ ( x − x ) + ( y i =1 2 1 2 1 2 2 − y2 ) + ( z1 − z2 )  ;  (Könnyen belátható, hogy ez a távolságfüggvény teljesíti a három metrika tulajdonságot.)

legyen az X az [ a, b] zárt intervallum felett (négyzetesen is) integrálható, valósértékű függvények lineáris tere, és legyen x, y ∈ X . Akkor az x és y függvények távolsága a korábban „5” számmal jelölt, euklideszi norma felhasználásával: d ( x, y ) = + b ∫  x ( t ) − y ( t ) 2 dt . a (Könnyen belátható, hogy ez a távolságfüggvény is teljesíti a három metrika tulajdonságot.) Kitekintés: a távolság fogalma kiterjeszthető, definiálhatjuk például az A, B ⊂ X részhalmazok egymástól mért távolságát (ez a két részhalmaz egymáshoz legközelebbi elemének távolsága): d ( A, B ) = inf {d ( x, y ) : x ∈ A, y ∈ B} ; Definiálhatjuk továbbá egy kiválasztott elem ( x ∈ X ) és egy részhalmaz ( A ⊂ X ) távolságát (ez a kiválasztott elem és a részhalmaznak a kiválasztott elemhez legközelebb eső elemének a távolsága): d ( x, A ) = inf {d ( x, y ) : y ∈ A} ; Meghatározhatjuk végül egy

részhalmaz ( A ⊂ X ) átmérőjét (ez a részhalmaz egymástól legtávolabb eső két elemének a távolsága): diam ( A ) = sup {d ( x, y ) : x ∈ A, y ∈ A} ; 10 I. Vizsgáljuk meg a munka-környezetünket I.7 Sorozatok, konvergencia, határérték, teljes tér Tekintsük az ℜ valós számtest felett értelmezett X normált lineáris tér elemeinek valamilyen xi ∈ X , i = 1, 2 ∞ végtelen sorozatát: az előző pont értelmében: d ( xm , x n ) := x m − x n , vagyis az elemek különb- ségének normája az elemek távolságaként értelmezhető; ha ∃ x ∈ X , hogy x − x n 0 hacsak n ∞ , akkor a sorozat konvergens és határértéke az x ∈ X elem; jelölése: lim x n = x , n∞ az xi ∈ X , i = 1, 2 elemek alkotta sorozatot Cauchy sorozatnak vagy önmagában konvergens sorozatnak nevezzük, ha tetszőlegesen kicsi ε > 0 számhoz található olyan M = M (ε ) küszöbindex, hogy ha m, n ≥ M akkor teljesül, hogy x m − x n < ε , vagyis

x m − xn 0 , ha n ∞ ; a gyakorlati számításokban tehát – elegendően nagy m, n esetén akár x m , akár x n tekinthető a sorozat határértékének; fontos tétel, hogy az {xi } ℜ -beli vagy ℜn -beli sorozat akkor és csak akkor konvergens, ha {xi } Cauchy konvergens. A gyakorlati számításokban a Cauchy konvergencia alapvetően fontos, hiszen a numerikus, közelítő megoldás általában ilyen módon határozható meg: az x határérték helyett meg kell elégednünk az x m és x n meghatározásával. Az x m és az x n egymástól mért távolságát (jó esetben közeledését, a távolságuk csökkenését) kísérjük figyelemmel, amikor egy numerikus megoldás „konvergencia történet”-ét vizsgáljuk. Ki kell még mondani azt is, hogy lényegében minden gyakorlati számításban létezik egy legkisebb, de nem nulla távolság korlát, aminél közelebb az x m és az x n nem kerül (nem kerülhet) egymáshoz – a gyakorlati számolást akkor

érdemes, illetve célszerű befejezni, ha ezt a helyzetet elértük, vagy legalább elegendő mértékben megközelítettük. Az X normált teret teljesnek nevezzük, ha e térben minden önmagában konvergens sorozatnak van határértéke. A teljes, lineáris, normált teret Banach térnek nevezzük. Ez a tulajdonság azért fontos, mert ha nem lenne határérték, akkor nincs mihez konvergálni. Másrészt a határérték létezését gyakran fizikai alapon döntjük el – mivel a feladatunk általában a valóságban létező / működő esetre vonatkozik, ezért feltehető, hogy korrekt modell esetében van határérték. 11 I. Vizsgáljuk meg a munka-környezetünket I.8 Skalár vagy belső szorzat Tekintsünk az ℜ valós számtest feletti X valós lineáris tér elempárjain egy S ( x, y ) ∈ℜ , az x, y elempár kétváltozós függvényét. Ez – mivel a valós számtestre képez le – egy funkcionál Értelmezzük az x, y elempár skalár szorzatát (más

elnevezéssel: belső szorzatát) az alábbi módon: ∀ x, y ∈ X : ∃ S ( x, y ) ∈ℜ+ 1./ S ( x, x ) ≥ 0; és 2./ S ( x, y ) = S ( y, x ) ; ahol: ℜ+ = [ 0, ∞ ) S ( x, x ) = 0 ⇔ x = ∅; 3./ S ( x + y, z ) = S ( x, z ) + S ( y, z ) , x, y, z ∈ X ; 4./ S (α x, y ) = α S ( x, y ) , α ∈ ℜ Példák: tekintsük az x, y ∈ X , x = ( x1 , x2 xi xn ) és y = ( y1 , y2 yi yn ) rendezett valós szám n-eseket és legyen definíció szerint def n S ( x, y ) = ∑ xi yi függvény – ez skalár szorzat, mert teljesülnek rá a i =1 fent előírt feltételek (így értelmezzük például az n = 3 esetben a 3Ds vektorok skaláris szorzatát); tekintsük az x = x ( t ) és az y = y ( t ) , x, y ∈ X elemeket, legyenek ezek az [ a, b] zárt intervallum felett négyzetesen is integrálható függvények – definiáljuk a skalár szorzatukat: def b S ( x, y ) = ∫ x ( t ) y ( t ) dt ; a ez is skalár szorzat, mert teljesülnek rá a fent előírt

feltételek. Megjegyzések: az X lineáris teret, ha skalár szorzatot is definiálunk rajta, unitérnek nevezzük; legyen az x ∈ X normája: x = + S ( x, x ) , ezzel az X valós pre-Hilbert tér. 12 I. Vizsgáljuk meg a munka-környezetünket I.9 Ortogonalitás Az x, y ∈ X elemeket ortogonálisnak nevezzük, ha S ( x, y ) = 0 , vagyis, ha a skalár szorzatuk nulla. (Az X tér zérus-eleme a tér minden más elemére ortogonális, azaz merőleges) Megjegyzés: közismert és szemléletes példa, hogy a két- és háromdimenziós vektorok ortogonálisak vagy merőlegesek egymásra, ha a skalár szorzatuk nulla. De az ortogonalitás olyan, általános fogalom, amelyet rengeteg más esetre is – például függvényterek elemeire is – kiterjesztettek Példa: tekintsük a [ −1,1] zárt intervallum felett értelmezett, Csebisev polinomokat (az I.1 ábrán csak az első három polinomot tüntettük fel): T0 ( t ) = 1 ; T1 ( t ) = t ; T2 ( t ) = 2t 2 − 1 ;. Tn +1 ( t ) =

2Tn ( t ) − Tn −1 ( t ) ; ezek a polinomok páronként ortogonálisak, mert: 1 ∫ T ( t ) T ( t ) dt = 0 i j ha i ≠ j −1 I.1 ábra – Csebisev polinomok „n” darab ilyen polinom egy n-dimenziós függvényteret feszít ki (a polinomok alkotják a függvénytér bázisát – és ezen bázis elemek lineáris kombinációjával a tér minden eleme egyértelműen előállítható)! 13 I. Vizsgáljuk meg a munka-környezetünket I.10 Minimalizáló vektor, minimalizáló elem Kezdjük a kérdéskör vizsgálatát egy, a közvetlen szemlélet alapján belátható példával: Legyen H = ℜ3 háromdimenziós, Descartes féle koordinátarendszerben értelmezett vektortér, és legyen M = ℜ2 ⊂ ℜ3 kétdimenziós (sík) altér. Legyen a kiválasztott sík az " x, y " koordináta tengelyek által kifeszített sík. I.2 ábra – Minimalizáló vektor Az I.2 ábrán feltüntettünk egy r ∈ H térbeli vektort és ennek a kiválasztott síkra vett,

merőleges vetületét: m0 ∈ M . A geometriai szemlélet alapján közvetlenül belátható, hogy az m0 ∈ M , a kijelölt síkban minimalizáló vektora r ∈ H -nak, hiszen, minden más m ∈ M -re igaz, hogy: r − m > r − m0 . Ez tehát azt jelenti, hogy a kiválasztott altérben m0 van az r -hez legközelebb, a minimalizáló vektort kiszámítani pedig az S ( r − m0 , m0 ) = 0 skalár szorzat segítségével lehet. Adjuk meg a koordinátáival az I.2 ábrán vázolt vektort: r = i xr + j yr + k zr Erre az esetre rögtön belátható, hogy m0 = i xr + j yr és r − m0 = k zr (e két vektor skalár szorzata valóban nulla). Vagyis a minimalizáló vektort merőleges (ortogonális) vetítéssel kaptuk meg A későbbiekben bemutatandó Galjorkin eljárás lényege is ez: ennél a módszernél a közelítő megoldást úgy számoljuk ki, hogy a hibafüggvény vetülete annak a függvénytérnek a bázisfüggvényeire (amely térben a közelítő megoldást keressük)

zérus legyen. Azaz a hibafüggvény legyen minden bázisfüggvényre ortogonális – legyen a hibafüggvény és a bázisfüggvények skalár szorzata zérus! 14 II. Az alap-feladat és néhány megoldási módszer II. Az alap-feladat és néhány megoldási módszer II.1 Az alap-feladat Legyen: adott a feladat matematikai modellje; adott a tartomány ( G ), amely felett a megoldást keressük; adottak a tartomány peremén ( ∂G ) a peremfeltételek; ha szükségesek, akkor adottak a kezdeti feltételek. Mi a következő feladattal foglalkozunk: legyen az L lineáris operátor, értelmezési tartománya legyen U , és az operátor képtere legyen F . Vagyis az L : U F leképezést vizsgáljuk. Az alap-feladat tehát: Lu = f ahol: u = u ( x ) ∈ U , az ismeretlen függvény; x ∈ G ⊆ ℜn , a független változó(k); f ∈ F , adott függvény. A peremfeltételek pedig a ∂G = ∂G1 ∪ ∂G2 ∪ ∂G3 peremen: 1./ x ∈ ∂G1 : u := u0 ( x ) : Dirichlet

féle, vagy elsőfajú perem (peremérték feladat), a keresett függvény értéke adott a ∂G1 perem-részen; 2./ x ∈ ∂G2 : ∂u := g 0 ( x ) : Neumann féle, vagy másodfajú perem (pe∂n remérték feladat), a gradiens adott a ∂G2 perem-részen; 3./ x ∈ ∂G3 : az 1/ és 2/ együtt: Cauchy féle vagy harmadfajú vagy vegyes perem (peremérték feladat) Az x n változót foglal össze. Időben változó feladat esetén ezek közül az egyik lehet az idő (t). Ezért a következőkben, általánosságban csak peremértékekről beszélünk Az időre vonatkozó kezdeti érték általában tehát peremértékként értelmezhető Egyes, konkrét feladatok esetén azonban, amikor az úgy kézenfekvő, használhatjuk az időt, mint elnevezésével is megkülönböztetett, független változót és ezekben az esetekben kezdeti feltételeket is meg kell adunk. Egy feladat korrekt kitűzésű, ha: létezik a megoldása (egzisztencia); a megoldás egyértelmű (unicitás); a

megoldás stabil (bármely adat kicsiny megváltozása a megoldást is csak kicsivel változtatja meg). 15 II. Az alap-feladat és néhány megoldási módszer Megjegyzések: a megoldás létezéséről és egyértelműségéről gyakran fizikai oldalról igyekszünk megbizonyosodni, illetve pontosan ezért a konkrét feladatok numerikus megoldását valamilyen módon ellenőrizni – validálni – kell; egy feladat rendelkezhet akár egyféle peremmel is (pl. lehet a teljes ∂G perem elsőfajú, akkor ∂G = ∂G1 , illetve ebben az esetben ∂G2 = ∂G3 = ∅ ); a ∂G1 , ∂G2 és ∂G3 perem-részek páronként diszjunktak; a kezdeti- és peremfeltételek szabatos megadása rendkívül fontos: nem megfelelő értékek megadásával esetleg látványos, de tökéletesen eltérő, adott esetben teljesen használhatatlan megoldást kapunk (a validálás emiatt is fontos)! II.2 A megoldás lehetőségei Az Lu = f feladatot numerikusan kívánjuk megoldani. Erre alapvetően

kettő plusz egy út kínálkozik: a megoldás-függvény közelítése, approximációja – ebben a munkában ezzel a módszer-csoporttal foglalkozunk; az operátor átírása, amikoris a megoldást diszkrét pontokban közelítjük (ilyen átírás lehet pl. egy differenciálegyenlet differenciaegyenletté történő átírása) – ezzel a módszer-csoporttal ebben a munkában nem foglalkozunk; az első két módszer kombinálása (ilyen lehet például a végeselem módszer, ahol a tartományt is diszkretizáljuk, az operátort is átírjuk és a megoldás-függvényt is közelítjük) – ebben a munkában az ebbe a módszer csoportba sorolható feladattal sem foglalkozunk. Megjegyzendő, hogy a két utóbbi módszer-csoport rendkívül széles szakirodalommal rendelkezik – az érdeklődő Olvasónak e szakirodalom tanulmányozását javasoljuk. Tekintsük a megoldás-függvény approximációját. Az első út (a megoldásfüggvény közelítése) nagy vonalakban azt

jelenti, hogy a G tartomány felett, valamilyen módon megválasztott bázisfüggvények segítségével egy függvényteret feszítünk ki. A közelítő megoldást pedig a bázisfüggvények lineáris kombinációjaként keressük, úgy, hogy a közelítő függvény valamilyen értelemben a megoldás-függvényhez a legközelebb legyen. Ezt a kérdéskört a következőkben vizsgáljuk részletesebben 16 II. Az alap-feladat és néhány megoldási módszer A második lehetséges út választása azt jelenti, hogy a G tartományt diszkretizáljuk – például ráccsal fedjük le – és a megoldást ezekben a diszkrét (tehát véges számosságú) pontokban keressük. Ehhez az operátor megfelelő átírásán keresztül vezet az út. Ennek a rendkívül nagy területet felölelő módszer-csoportnak a tanulmányozására a vonatkozó szakirodalmat ajánljuk A harmadik út nyilvánvalóan az első két út kombinációja. Ebből a kérdéskörből jellemző példa a végeselem

módszer, amelyet szintén igen kiterjedt szakirodalom tárgyal. II.3 A megoldás-függvény közelítése Legyen: u ∈ U (a tényleges megoldás); és uɶ ∈ U n , U n ⊆ U (a közelítő megoldás); illetve ϕi , i = 1, 2 n az U n függvénytér bázisfüggvény-rendszere. n Az U n altér ∀uɶ ∈ U n eleme egyértelműen előállítható az uɶ = ∑ α i ϕi lineáris i =1 kombinációval. A feladatatunk az α i ∈ℜ együtthatók azon értékeinek megkeresése, amely értékekre az f − Luɶ = ε eltérés (hiba) a megadott peremfeltételek kielégítése mellett, (valamilyen értelemben) minimális lesz! A keresett α i együttható értékeket legtöbbször inhomogén, lineáris algebrai egyenletrendszer megoldásával kaphatjuk meg. A megoldás és a közelítő megoldás eltérése =+ ∫ ( u − uɶ ) 2 u − uɶ = S ( u − uɶ ) , ( u − uɶ )  = dx ; mi ezzel szemben az Lu − Luɶ = f − Luɶ eltérést vizsgáljuk, il- G letve ennek az

eltérésnek keressük a minimumát. Alapvető kérdés az, hogy ha a közelítő megoldás képét közelítjük a megoldás képéhez, akkor ebből következik-e, hogy a közelítő megoldás a megoldáshoz tart? A kérdésre válasz a vonatkozó szakirodalom tanulmányozásával kapható. Mérnöki oldalról közelítve a kérdést azt mondhatjuk, hogy ezért is rendkívül fontos a validáció! Elöljáróban bemutatunk néhány szempontot, amit ennél a módszer típusnál – lehetőség szerint - célszerű szem előtt tartani: a bázisfüggvények egyenként és lineáris kombinációjukkal is elégítsék ki a peremfeltételeket (ez legtöbbször csak iskola-példákban lehetséges); 17 II. Az alap-feladat és néhány megoldási módszer a bázisfüggvények legyenek páronként ortogonálisak S (ϕi , ϕ j ) = 0 ha i ≠ j ; a bázisfüggvények legyenek az operátor sajátfüggvényei: L ϕi := α ϕi , α ∈ℜ . Ezeket a szempontokat igen célszerű szem előtt

tartani, pontosan betartani azonban leginkább szintén iskola-példákban lehetséges. A megoldás-függvény közelítésére sok módszer létezik. Ebben a munkában ezek közül csak hárommal, a kollokációs módszerrel, a legkisebb négyzetek módszerével és a Galjorkin eljárással (Галёркин, Борис Григорьевич (1871 – 1945)) foglalkozunk 18 III. A minta-feladat III. A minta-feladat Az ebben a fejezetben vizsgálandó módszerek lényegének lehető alapos bemutatása érdekében egy minta-feladatot dolgoztunk ki, mely feladatot a bemutatott módszerekkel meg is oldunk. Hangsúlyozzuk, hogy a minta-feladat egy elméleti probléma, illetve annak vizsgálata – ebből a vizsgálatból gyakorlati következtetéseket csak kellő körültekintéssel szabad levonni. A következőkben, a vizsgálatainkban Descartes féle, derékszögű koordinátarendszert ( x, y, z ) alkalmazzuk – csak megjegyezzük, hogy más koordinátarendszerek is léteznek,

illetve adott esetben célszerűen alkalmazhatók is. III.1 A modell folyadék A minta-feladat egy, elméleti áramlástani probléma kitűzését és megoldását tartalmazza. Elméletinek nevezzük, mert a valóságos helyett modell folyadékot (modell közeget) alkalmazunk Ez a modell folyadék a valóságos közegeket alkotó részecskék helyett kontinuum, amely a teret matematikai értelemben folytonosan tölti ki. Így lehetséges, hogy a közeget, illetve annak áramlását jellemző függvények adhatók meg. Fontos lehet például a sűrűség mező ( ρ = ρ ( x, y, z, t ) ), a nyomás mező ( p = p ( x, y, z, t ) ), a sebesség mező ( c = c ( x, y, z, t ) ) és több, további függvény. A kontinuum modellen túlmenően azt is kikötjük, hogy a minta-feladatban szereplő közeg legyen ideális. Az ideális (tökéletes) közeg homogén (azaz egynemű), izotróp (azaz minden irányban egyformán viselkedik), összenyomhatatlan (állandó sűrűségű) és

súrlódásmentes (tökéletesen gördülékeny, vagyis csak nyomófeszültség keletkezik benne. III.2 A Hamilton féle nabla-operátor és a Laplace-operátor A Hamilton féle nabla-operátornak egy szimbolikus „ ∇ ”-val jelölt, vektoroperátort nevezünk, amelyet derékszögű koordinátarendszerben a következőképpen írhatunk: 19 III. A minta-feladat ∇=i ∂ ∂ ∂ +j +k . ∂x ∂y ∂z A nabla-operátort skalár-vektor függvényre alkalmazva a gradiens vektort kapjuk. Tekintsük például a nyomás gradiensét: grad p = ∇p ; ez derékszögű koordinátarendszerben a ∇p = i ∂p ∂p ∂p +j +k alakot ölti. A gradiens vektor ∂x ∂y ∂z megmutatja a nyomás (vagy más, skalár-vektor függvény) változásának irányát és mértékét. E művelet eredménye tehát egy vektor Ha a skalármennyiségeket nulla-indexes vagy nulladrendű tenzornak, a vektorokat egy-indexes vagy elsőrendű tenzornak tekintjük, akkor a nabla-operátor fenti típusú

alkalmazásával nulla-indexesből egy-indexes tenzort kapunk. A nabla-operátor vektor-vektor függvényre skaláris vagy belső szorzattal történő alkalmazásával a divergenciát kapjuk (ez skalármennyiség – ez a művelet tehát egy-indexesből nulla-indexes tenzorhoz vezet). Tekintsük például a  cx    ∂ ∂ ∂ sebességtér divergenciáját div c = ∇T c =  , ,  c y  . A skaláris vagy belső  ∂ x ∂ y ∂ z  c   z szorzást div c = a sor-oszlop ∂ cy kombináció szerint elvégezve kapjuk, hogy: ∂ cx ∂c + + z . (Ezt ismét derékszögű koordinátarendszerben határoz∂x ∂y ∂z tuk meg.) Állandó sűrűségű közeg esetén – a folytonosság törvénye szerint (az anyagmegmaradás elve alapján) ha a vizsgált áramlási térben nincs sem forrás sem nyelő – a sebességtér divergenciája zérus ( div c = 0 ). A nabla-operátor vektor-vektor függvényre vektori szorzattal történő

alkalmazásával a rotáció vektort kapjuk. Ebben az esetben tehát egy-indexes tenzorból indulunk és az eredmény is egy-indexes tenzor. Tekintsük ismét a sebességteret, ennek rotációja derékszögű koordináta rendszerben:  i  rot c = ∇ × c =  ∂  ∂x  cx   ∂ cz ∂ c y  −   j k  ∂y ∂z     ∂ ∂  =  ∂ cx − ∂ cz  ∂y ∂z   ∂z ∂x  cy cz   ∂ c ∂ c   y − x  ∂ x ∂ y  A rotáció fizikai értelemben a vizsgált pontban értelmezhető, a merevtestszerű forgást jellemző szögsebesség kétszerese. Vagyis a rotáció nem nulla, ha 20 III. A minta-feladat a modell folyadék – kontinuum – vizsgált pontbeli, matematikai értelemben vett részecskéje forog, tehát a szögsebessége nem nulla. Csak megjegyezzük, hogy a diadikus vagy külső szorzat a derivált tenzort eredményezi: D = c ∇T (egyindexesből kétindexes tenzort kapunk). A

derivált tenzort egy szimmetrikus és egy ferdén szimmetrikus tenzor összegére szokás bontani [ D = DS + D A , ahol DS = ( D + DT ) / 2 és D A = ( D − DT ) / 2 ]. DS -et, a fizikai tartalma alapján alakváltozási-sebesség tenzornak, D A -t pedig, szintén a fizikai tartalma alapján örvény tenzornak nevezzük [3]. A Laplace operátort ( ∆ ) úgy kapjuk, hogy a nabla operátort önmagával skalárisan megszorozzuk (ezért néha – véleményünk szerint nem túl szerencsés módon – nabla-négyzetnek is nevezik). Derékszögű koordinátarendszerben: ∆ = ∇T ∇ = ∂2 ∂2 ∂2 + + . ∂ x2 ∂ y 2 ∂ z 2 A Laplace operátor jelentősége igen nagy, sok-sok hőtani, áramlástani, illetve számos további, például mechanikai feladatban találkozunk vele. Fontos leszögezni, hogy létezik (található) olyan függvény-pár, amelyre az ∆x − ∆y > x − y , vagyis a Laplace operátor nem kontrakciós. III.3 A minta-feladat Tekintsünk egy

egység-négyzet keresztmetszetű (III.1 ábra), mindkét irányban végtelenbe nyúló csövet, melyben nincs alkotó irányú (hosszirányú) áramlás, de minden, az alkotóira merőleges keresztmetszetében azonos, örvényes, stacionárius (időtálló, időben állandó) síkáramlás jön létre. Válasszunk ki egy keresztmetszetet – ezt tüntettük fel a III.1 ábrán A továbbiakban, az ebben a keresztmetszetben létrejövő (a többi keresztmetszettel azonos) stacionárius síkáramlást vizsgáljuk (Természetesen, ez az áramlás magától nem jön létre, valamilyen módon létre kell hoznunk – pl. valahogyan megkeverjük a folyadékunkat.) Fizikai alapon elvárjuk, hogy ebben az áramlásban a ψ = áll. áramvonalak önmagukba záródó görbesereget alkossanak (III.4 ábra) Hasonlóképpen, fizikai belátás alapján elvárjuk, hogy a kialakuló áramlás jellemzői ciklikusan változzanak (Ez azt jelenti, hogy egy – tetszőleges – áramvonal tetszőleges

pontjába beérkező matematikai részecske (azaz egy, geometriai pont) áramlástani 21 III. A minta-feladat jellemzői (pl. sebesség, statikus nyomás stb) akárhány körüláramlás után is azonosak maradjanak - stacionárius áramlás – miközben e jellemzők egyébként pontról-pontra változnak.) A minta-feladatban legyen tehát a „G” tartomány az egység-négyzet, azaz: G = {0 ≤ x1 ≤ 1 ∧ 0 ≤ x2 ≤ 1} ; illetve vezessük be az x1 = x és az x2 = y jelölést, akkor: G = {0 ≤ x ≤ 1 ∧ 0 ≤ y ≤ 1} ; (III.1 ábra) A ∂G peremet (a „csőfalak”, azaz az egység-négyzet oldalai) vastag vonallal jelöltük. III.1 ábra – Számítási tartomány Az áramlástanban az örvényes síkáramlást az áramfüggvényre [ψ = ψ ( x, y ) ] felírt, Poisson féle, másodrendű, lineáris, elliptikus, parciális differenciálegyenlet írja le: ∂ 2ψ ∂ 2ψ ∆ψ = −ω , részletesen kiírva: + = −ω ( x, y ) ; 2 2 ∂x ∂y A fenti

matematikai modellben az „ ω ” az örvényességet jelenti. Az örvényesség, azaz a sebességmező rotációja fizikai jelentése szerint a helyi szögsebesség kétszerese Mi itt az örvényességet konkrétan megadjuk (előírjuk); a részletesen kiírt matematikai modell tehát: ∆ψ = ∂ 2ψ ∂ 2ψ + = −2  y (1 − y ) + x (1 − x )  ∂ x2 ∂ y2 (3.1) ahol az örvényesség: ω ( x, y ) = 2  y (1 − y ) + x (1 − x ) . A modellhez tartoznak még a peremfeltételek. Mivel a feladat kitűzése szerint a ∂G perem egyetlen, zárt áramvonal, az áramfüggvény értéke a peremen állandó. Ezt az állandó értéket megválaszthatjuk: legyen az áramfüggvény értéke a peremen azonosan nulla: ψ ( x, y ) = 0, ha { x, y} ∈ ∂G (ez elsőfajú, azaz Dirichlet féle perem). Vizsgáljuk meg az örvényességet. Fizikai értelemben – mint már említettük - az örvényesség a merevtest-szerű forgás szögsebességének kétszerese.

Csak megemlítjük, hogy a folyadékokban (deformálható közegekben) szögdefor- 22 III. A minta-feladat mációt (disztorziót), illetve szögdeformáció-sebességet is értelmezünk. Ezt, a szögdeformáció-sebességet és a hosszváltozási (dilatációs) sebességet tartalmazza a DS = ( D + DT ) / 2 , alakváltozási-sebesség tenzor. Ábrázoljuk az örvényességet a G tartomány, azaz az egység négyzet felett (III.2 ábra) Írjunk egy SCILAB eljárást (mi az 551 verziót használtuk, de a program-csomagot folyamatosan fejlesztik), ami a G tartomány feletti örvényességet felületként ábrázolja. (A konkrét eljárás az eredetileg folytonos felületet véges darabszámú rész-négyszögekkel közelítve rajzolja ki) Az eljárás program-listája: clear deff(omega=fv(x,y),omega=2*(y(1-y)+x(1-x))) x=[0:0.1:1]; y=x; fplot3d1(x,y,fv) - változók törlése; - az örvényesség-függvény; - a G kitűzése és felosztása; - „3d” rajz Ø III.2 ábra III.2

ábra – Az örvényesség Az örvényesség-függvény tehát egy, G feletti felület, ennek közelítő képe látható a III.2 ábrán Csak megjegyezzük, hogy pl a G finomabb felosztásával nyilván pontosabban közelítő felület állítható elő – nekünk, szemléltetésre ez a felosztási finomság megfelelő! A második megjegyzés technikai jellegű: az eljárás lefuttatása után a III.2 ábrán látható kép egy, grafikus ablakban jelenik meg. Ennek az ablaknak a fejléce számos lehetőséget kínál, ezek között van a megjelenő kép exportálásának lehetősége (is). A szóban forgó ábra fő része így készült. Végül felhívjuk a figyelmet arra, hogy az exportált képet kiegészítettük az „x”, az „y” és a felfele mutató „z” – itt ez konkrétan az ω ( x, y ) - tengely 23 III. A minta-feladat megadásával. Az eredeti (kiegészítetlen) ábrán az „x” tengely skálázása a láthatóság miatt önmagával párhuzamosan

eltolva, az ábrán jobbra lent látható Hasonló megfontolás alapján az „y” tengely skálázása eltolva, az ábrán balra lent látható. A program a „z” felosztását jelző skálázást is (automatikusan) eltolta A III.2 ábrára pillantva megállapítható, hogy az örvényesség a sarokpontokban nulla, illetve a (05,05) pontban maximális ( ω ( 05, 05) = 1 ); a sarokpontok kivételével pedig pozitív Ez – jobbsodrású koordináta rendszerben, felülnézetből tekintve – az óramutató járásával ellentétes forgásirányt jelent Illetve ez megfelel annak, hogy a G tartomány pozitív körüljárása során (a tartomány peremén egyensúlyozva) a bal kéz esik a tartomány belseje felé. A későbbiekben megmutatjuk, hogy a G tartományban kialakuló áramlás sebességeloszlása valóban eszerint alakul A feladat egy iskola-példa, úgy konstruáltuk, hogy létezzen zárt alakú megoldása: ψ ( x, y ) = x (1 − x ) y (1 − y ) + állandó Könnyen

belátható, hogy – a korábban megadott peremfeltétel miatt – a zárt alakú megoldás – ami egyébként a vizsgálni kívánt áramlás áramfüggvénye – kifejezésében szereplő állandó értéke zérus, tehát a zárt alakú megoldás a mi esetünkben: ψ ( x, y ) = x (1 − x ) y (1 − y ) (3.2) Azt, hogy ez valóban megoldás, azt a ∂ 2ψ ∂ 2ψ + = −2  y (1 − y ) + x (1 − x )  , a ∂ x2 ∂ y 2 minta-feladat szerinti differenciálegyenletbe történő visszahelyettesítéssel láthatjuk be. Differenciáljuk az áramfüggvényt parciálisan „x”-szerint: ∂ 2ψ ∂ψ = (1 − 2 x ) y (1 − y ) és = −2 y (1 − y ) . 2 ∂x ∂x Az „y” szerinti, második parciális derivált, értelemszerűen: ∂ 2ψ ∂y 2 = −2 x (1 − x ) A feladat kitűzése szerint a ∂G perem egyetlen, zárt áramvonal, az áramfüggvény értéke a peremen azonosan nulla – a fent megadott áramfüggvény teljesíti ezt a feltételt. Innen pedig

már nyilvánvaló, hogy a megadott áramfüggvény valóban zárt alakú megoldása a minta-feladatnak 24 III. A minta-feladat Az a tény, hogy ismerjük a kitűzött feladat zárt alakú megoldását azért fontos, mert így közvetlenül is ellenőrizhetjük a közelítő megoldásokat. Illetve a numerikus módszereknél igen fontos a validálásnak nevezett eljárás: ennek során megpróbáljuk a lehető legalaposabban igazolni, hogy a szóban forgó numerikus módszer helyes eredményre vezet-e. Az ismert megoldás igen alkalmas validálási lehetőséget jelent III.4 A zárt alakú megoldás vizsgálata Írjunk ismét SCILAB eljárást, amely a G tartomány felett az áramfüggvényt, (ψ ( x, y ) ) mint felületet ábrázolja: clear deff(aramfv=fv(x,y),aramfv=(y*(1-y)x(1-x))) x=[0:0.1:1]; y=x; fplot3d1(x,y,fv) - változók törlése; - az áramfüggvény; - a G kitűzése és felosztása; - „3d” rajz Ø III.3 ábra III.3 ábra – Az áramfüggvény A (3.1)

parciális differenciálegyenlet a minta-feladatunkban az áramfüggvény-felületet (III3 ábra) képezi le az örvényesség-felületre (III2) ábra Válasszunk ki egy tetszőleges { x, y} ∈ G koordináta párt Ez a koordináta pár pontosan egy áramfüggvény pontot kapcsol pontosan egy örvényesség ponthoz, vagy megfordítva pontosan egy örvényesség pontot kapcsol pontosan egy áramfüggvény ponthoz. Továbbá igaz az is, hogy bármely vagy az összes 25 III. A minta-feladat áramfüggvény-felület, illetve örvényesség felület pont-párhoz rendelhető koordináta pár. Ezzel tehát az áramfüggvény-felület és az örvényesség felület pontjai között egy kölcsönösen egyértelmű (bijektív) kapcsolatot adtunk meg. Végeredményben tehát (3.1) bijektív leképezés A III.3 ábráról leolvasható (vagy a (32) kifejezésből megállapítható), hogy a ψ ( x, y ) áramfüggvény értéke a ∂G -n, azaz a peremen – a feladat kitűzésének megfelelően

– azonosan nulla. Ezek szerint a ∂G valóban áramvonal A kialakuló áramlás további ψ = áll. áramvonalait megkapjuk, ha a III3 ábrán látható felületet az ( x, y ) síkkal párhuzamos síkkal metsszük Írjunk SCILAB eljárást, amely a G tartományon futó, például 5 áramvonalat ábrázolja: clear deff(aramfv=fv(x,y),aramfv=(y*(1-y)x(1-x))) x=[0:0.01:1]; y=x; fcontour(x,y,fv,5) - változók törlése; - az áramfüggvény; - a G kitűzése és felosztása; - 5 db. áramvonal Ø III4 ábra III.4 ábra – Az áramvonalak 26 III. A minta-feladat A III.4 ábrán látható öt, a G belsejébe eső áramvonal (az egyes áramvonalak mellé írt számok az áramfüggvény értékét jelentik; ψ max = ψ ( 0.5,05 ) = 00625 ), illetve az ábrát kiegészítettük a hatodik, ψ = 0 áramvonallal (vastag, piros vonal – ez a ∂G perem). Az ábra alapján megállapítható, hogy az áramvonalak a tartomány középpontjához tartva körhöz, a perem felé tartva az

(egység) négyzethez tartanak. Vizsgáljuk meg a kialakuló áramlás sebesség-eloszlását is. A sebesség-mező jelen esetben kétkomponensű vektor-mező, az egyes sebesség összetevők az áramfüggvény parciális deriválásával számíthatók: cx = ∂ψ = x (1 − x )(1 − 2 y ) ∂y és cy = − ∂ψ = − (1 − 2 x ) y (1 − y ) ; ∂x Írjunk SCILAB eljárást, amely a sebesség-mezőt ábrázolja: - változók törlése; - felosztás szám; - a G kitűzése és felosztása; clear n=17 for i=1:n x(i)=(i-1)/(n-1); y(i)=(i-1)/(n-1); end for i=1:n for j=1:n cx(i,j)=(1-2*y(j))x(i)(1-x(i)); cy(i,j)=(2*x(i)-1)y(j)(1-y(j)); end end champ(x,y,cx,cy) - sebesség összetevők számítása; - ábrázolás Ø III.5 ábra III.5 ábra – Sebesség-mező 27 III. A minta-feladat A III.5 ábrán látható sebesség vektorok (nyilak halmaza) mutatja a kialakuló áramlás irányát és a nyilak hosszúsága alapján a sebességek nagyságát is Leolvasható az

ábráról továbbá, a korábban már említett pozitív, felülnézetben az óramutató járásával ellentétes forgásirány is. Rajzoltunk egy további ábrát is (III.6 ábra), amelyen a sebesség-eloszlás néhány jellegzetessége látható: III.6 ábra – Kiválasztott sebesség eloszlások Tekintsük a G tartomány középpontjának a (0.5,05) pontot (ez az átlók metszéspontja) A sebesség komponensek kifejezése alapján megállapítjuk, hogy ebben a pontban cx = c y = 0 , azaz a sebesség itt nulla. A III6 ábrán feltüntetett sebesség-profilok ennek megfelelően alakulnak. Például a III.6 ábra alapján a sebesség-eloszlásban egyfajta ciklikusság ismerhető fel Az x = 05, y ∈ [ 0,1] és az x ∈ [0,1] , y = 05 vonalak G tartománybeli szakaszai felett a sebesség-eloszlás lineáris: cx ( x = 0.5, y ∈ [ 0,1]) = 025 (1 − 2 y ) cx ( x ∈ [ 0,1] , y = 0.5 ) = 0 és c y ( x = 0.5, y ∈ [ 0,1]) = 0 ; valamint és c y ( x ∈ [ 0,1] , y = 0.5 ) = −025

(1 − 2 x ) ; Vizsgáljuk meg a (0.5,05) középponttól az (1,1) pontig tartó fél főátló pontjaiban értelmezhető sebesség-eloszlást Ebben az esetben az y = x helyettesítéssel élhetünk, így: 28 III. A minta-feladat cx = x (1 − x )(1 − 2 x ) és c y = − (1 − 2 x ) x (1 − x ) ; Megállapítható, hogy a két sebesség összetevő csak az előjelében különbözik, vagyis az eredő sebesség merőleges a vizsgált fél főátlóra. Megállapítható az is, hogy a sebesség a fél-főátló mindkét végpontjában nulla. Kiszámolható, hogy a vizsgált sebesség-eloszlásnak a (0.789,0789) pontban maximuma van Ezek a megállapítások értelemszerűen (!) igazak a további három fél-főátlóra is (ez látható a III.6 ábrán is) Beláttuk, hogy a sebesség a ∂G sarokpontjaiban nulla. A nulla vektornak pedig bármilyen iránya lehet. Ezért matematikai szempontból nem probléma az, hogy az áramlás a sarokpontokban „sarkon fordul”. Ha

ideális helyett valóságos közeget vizsgálnánk, akkor – természetesen – bizonyos rész tartományokban az itt vizsgálttól lényegesen eltérő áramlási jellemzőket látnánk 29 IV. A közelítő megoldás függvénytere és a hiba-függvény IV. A közelítő megoldások függvénytere és az eltérés (hiba) függvény IV.1 A közelítő megoldások függvénytere A közelítő megoldások megkereséséhez definiálnunk kell azt a függvényteret, amelyben dolgozni fogunk. Legyen ez az U (a 31 összefüggéssel definiált minta-feladatban szereplő Laplace operátor értelmezési tartománya) n dimenziós altere: ( U n ) n ψɶ ∈ U n ⊆ U , ψɶ = ∑ α i φi i =1 Legyen n = 5 , vagyis a következőkben dolgozzunk egy, öt-dimenziós altérben. A közelítő áramfüggvényt (ψ közelítő megoldás) az alábbi öt bázisfüggvény ( φi , i = 1, 2 5 ) lineáris kombinációjaként keressük, a (41) szerinti formában: (4.1) 5 ψɶ = ∑ α iφi , α i

∈ ℜ , ahol: i =1 φ1 = sin (π x ) sin (π y ) φ2 = sin ( 3π x ) sin ( 3π y ) φ3 = sin ( 5π x ) sin ( 5π y ) φ4 = sin ( 7π x ) sin ( 7π y ) φ5 = sin (11π x ) sin (11π y ) ∆φ1 = −2π 2φ1 ∆φ2 = −18π 2φ2 ∆φ3 = −50π 2φ3 ∆φ4 = −98π 2φ4 ∆φ5 = −242π 2φ5 A numerikus, közelítő megoldás megkeresése ebben az esetben az α i együtthatók azon számértékének a meghatározását jelenti, amelyekkel a közelítő megoldás a pontos megoldáshoz – valamilyen értelemben – a legközelebb lesz. Alkalmazzuk a Laplace operátort a (4.1) összefüggésnél megadott bázisfüggvényekre A művelet eredménye az egyes bázisfüggvények mellett látható E művelet eredménye minden bázisfüggvény esetében a bázisfüggvény konstans-szorosa, ezek a bázisfüggvények tehát a Laplace operátor sajátfüggvényei. 30 IV. A közelítő megoldás függvénytere és a hiba-függvény Minden bázisfüggvényre igaz, hogy φi ( x, y )

= 0, i = 1, 2 5, ha { x, y} ∈ ∂G , vagyis a perempontokban minden bázisfüggvény értéke zérus. Ezek szerint minden bázisfüggvény önmagában kielégíti az előírt (a peremen zérus) peremfeltételt. A közelítő megoldás a bázisfüggvények véges lineáris kombinációja, ez a kombináció pedig – véges együtthatók esetén - automatikusan kielégíti a kirótt peremfeltételt. Ez tehát azt jelenti, hogy a közelítő megoldások keresésekor a peremfeltétellel nem kell foglalkoznunk Összetettebb feladatok esetén ez a kérdés ilyen egyszerűen nem intézhető el! IV.1 ábra – A bázisfüggvények A választott bázisfüggvények igen előnyös tulajdonsága még, hogy páron1 1  0 0  ként ortogonálisak: S (φi , φ j ) = ∫  ∫ φi ( x, y ) φ j ( x, y ) dx  dy = 0, ha i ≠ j . Az i = j esetben – négyzetgyökvonással kiegészítve - a tekintett bázisfüggvény S (φi , φi ) = 1 2  ∫0  ∫0 φi ( x, y ) dx

 dy = φi 1 normáját (az azonosan nulla függvénytől mért távolságát) kapjuk. Kiszámítva a bázisfüggvények normáját, azt találjuk, hogy: φi = 0.707, i = 1, 2 5 , vagyis minden bázisfüggvény normája azonos Érdemes felhívni a figyelmet arra, hogy a fentiekben definiált függvénytér bizonyos vonatkozásokban analógiát mutat a közismert, általánosan használt három-dimenziós Descartes-féle koordináta rendszerrel (geometriai tér), amit szintén páronként ortogonális egységvektorok feszítenek ki. Illetve amelyben 31 IV. A közelítő megoldás függvénytere és a hiba-függvény tetszőleges vektor a három bázis-vektor lineáris kombinációjaként állítható elő. Az általunk definiált függvénytérrel kapcsolatban megjegyezzük, hogy enynyi előnyös tulajdonsággal rendelkező bázisfüggvény rendszer leginkább iskolai példában fordul elő! IV.2 Az eltérés (hiba) függvény a minta-feladat esetében E fejezet elején

definiáltuk az eltérés (hiba) függvény általános alakját: f − Luɶ = ε . Írjuk fel ezt most részletesen, a minta-feladat esetében, a megoldáshoz definiált, öt-dimenziós függvénytér elemeinek felhasználásával:  5  ε = ω + ∆ψɶ = ω + ∆  ∑ α iφi  , α i ∈ℜ , illetve:  i =1  ε = ω − α1 2π 2φ1 − α 218π 2φ2 − α 3 50π 2φ3 − α 4 78π 2φ4 − α 5 242π 2φ5 (4.2) A feladatunk a fenti eltérés (hiba) függvényben szereplő α i ( i = 1, 2 5 ) együtthatók megválasztása, úgy, hogy a közelítő megoldás pontos megoldástól való eltérése – valamilyen értelemben – legyen kicsi (esetleg legkisebb). Ebben a munkában három módszert (kollokációs módszer, a legkisebb négyzetek módszere (LKN) és a Galjorkin eljárás) vizsgálunk meg. Mint megállapítottuk, jelen iskolapéldában a közelítő megoldás általunk definiált függvényterének bázisfüggvényei egyúttal a Laplace operátor

sajátfüggvényei is. Ezért a (41) alakban keresett közelítő megoldás minden (bármely) véges {α i } együttható-csoport esetén, a G tartomány peremén azonosan zérus értéket vesz fel. Ezért a közelítő megoldással számolt közelítő örvényesség legalább a perempontokban – a négy sarokpont kivételével – a feladatban megadott örvényességtől mindig különbözni fog! Ezt a hibát – ezen a módszer-csoporton belül, az általunk választott függvénytér alkalmazása esetén – nem tudjuk elkerülni. 32 V. A kollokációs módszer V. A kollokációs módszer V.1 A kollokációs módszer leírása n A közelítő megoldást – mint már leszögeztük – az uɶ = ∑ α i ϕi formában keresi =1 sük. Az egyszerűség kedvéért tekintsük azt a minta-feladatnak megfelelő esetet, amikor a bázisfüggvények egyenként és lineáris kombinációjukkal is kielégítik a peremfeltételeket Ennél a módszernél kiválasztunk n számú pontot (

xi ∈ G, i = 1, 2, n ), és az α i ∈ℜ együtthatók értékét úgy állapítjuk meg, hogy a hiba vagy eltérés-függvény a felvett pontokban nulla értéket vegyen fel. Vagyis:  n  f ( xi ) − Luɶ ( xi ) = f ( xi ) − L  ∑ α i ϕi ( xi )  = ε ( xi ) := 0 i = 1, 2, n ;   i =1 A fenti kifejezés n darab algebrai egyenletet jelent. Oldjuk meg a mintafeladatot ezzel a módszerrel V.2 A kollokációs módszer alkalmazása Ennél a módszernél tehát első lépésként megfelelő számú – itt, konkrétan öt – pontot kell választanunk. Legyenek ezek rendre az ( x1 , y1 ) , ( x2 , y2 ) , ( x3 , y3 ) , ( x4 , y4 ) és ( x5 , y5 ) pontok. A minta-feladat esetében, a (42) egyenlet szerint, az i –edik kollokációs pontra felírható: α1 2π 2φ1 ( xi , yi ) + α 218π 2φ2 ( xi , yi ) + α 3 50π 2φ3 ( xi , yi ) + +α 4 78π 2φ4 ( xi , yi ) + α 5 242π 2φ5 ( xi , yi ) = ω ( xi , yi ) i = 1, 2 5; Ez egy inhomogén, lineáris, algebrai

egyenlet-rendszer; itt öt egyenletünk van, öt ismeretlenhez. A fenti egyenlet-rendszer mátrixos formában:  2π 2φ1 ( x1 , y1 )    2π 2φ1 ( x3 , y3 )    2  2π φ1 ( x5 , y5 ) 50π 2φ3 ( x1 , y1 ) 242π 2φ5 ( x1 , y1 )  α1   ω ( x1 , y1 )       α 2  ω ( x2 , y2 )  50π 2φ3 ( x3 , y3 ) 242π 2φ5 ( x3 , y3 )  α 3  = ω ( x3 , y3 )  ;      α 4  ω ( x4 , y4 )  50π 2φ3 ( x5 , y5 ) 242π 2φ5 ( x5 , y5 )  α 5  ω ( x5 , y5 )  azaz: (5.1) Aα = ω ; 33 V. A kollokációs módszer Válasszunk tehát öt, valóban különböző pontot. A különbözőség – mint majd a későbbiekből kiderül – igen fontos. Legyenek a választott pontok az I Mellékletben található program-lista szerint (a programban a kollokációs pontok koordinátáit „xk”-val, illetve „yk”-val jelöltük): xk(1)=0.25;

yk(1)=013; xk(2)=0.58; yk(2)=096; xk(3)=0.67; yk(3)=061; xk(4)=0.86; yk(4)=054; xk(5)=0.95; yk(5)=043; A programot a SCILAB környezetben lefuttatva, az α i együtthatókra a következő eredményeket kaptuk: α1 = 0.1061, α 2 = 0000087, α 3 = −0000778, α 4 = −00005914, α 5 = 00002137; Első ránézésre az együtthatók értékénél rokonszenves, hogy az első együttható legalább százszor nagyobb, mint a többi együttható abszolút értéke. Ezek szerint – amint az várható is – a közelítő megoldásban a φ1 bázisfüggvénynek a legnagyobb a szerepe. A futtatás eredményeiből kiderül az is, hogy az eltérés (hiba) függvény értéke a kollokációs pontokban jó közelítéssel nulla (10-15 nagyságrendű vagy kisebb). Számítsuk ki a jelen közelítő megoldás esetében az eltérés (hiba) függvény önmagával vett skalár szorzatát: (( )( S ( ε , ε ) = S ω + ∆ψ , ω + ∆ψ )) ≅ 0.3967 Az 1.6, illetve az 18 pont szerint az eltérés

(hiba) függvény Euklidész-i normája (ezt nevezhetjük a hiba nagyságának) pedig: ε = + S ( ε , ε ) ≅ 0.6299 Illetve a későbbi összehasonlítások érdekében határozzuk meg az áramfüggvény és a közelítő áramfüggvény különbségének önmagával vett skalár szorzatát, továbbá az Euklidész-i normáját – ez a norma mutatja az áramfüggvény és a közelítő áramfüggvény távolságát: (( )( S ψ −ψ , ψ − ψ )) ≅ 0.0003949 és 34 ψ −ψ ≅ 0.01987 V. A kollokációs módszer Az (5.1) inhomogén, lineáris, öt-ismeretlenes, algebrai egyenletrendszerből – a konkrét kollokációs pontok megválasztása (kijelölése) – után kiszámítottuk az α i együtthatókat. Ezzel előállítható a (41) szerinti közelítő áramfüggvény – ezt az V.1 ábrán tüntettük fel: V.1 ábra – Kollokációs megoldás, közelítő áramvonalak Vizsgáljuk meg a közelítő áramvonalak alakulását! A zárt alakú megoldásból

kapott áramvonalak a III.4 ábrán láthatók – a közelítő áramvonalak ezekhez valamennyire hasonlóak, bár alaposabb szemrevételezés után kisebb eltérések azért láthatók A leglényegesebb eltérés a legbelső (majdnem kör) áramvonalon látszik: a pontos megoldás esetén az áramfüggvény értéke ezen az áramvonalon 0.0521, a közelítő megoldásnál ez a szám 0.0876 Külön kiszámoltuk, hogy míg a ψ ( 0.5, 05 ) = 00625 , addig a ψɶ ( 05, 05 ) = 0105 - ez csaknem 70%-os eltérés! Vagyis megállapítható, hogy a megoldásunk „szemre” ugyan nem túl rossz (az áramlás jellegét valamennyire mutatja), azonban számszerűen elfogadhatatlan! (Érdemes erre a CFD feladatok megoldásának szemre történő validálásánál emlékezni!) 35 V. A kollokációs módszer A program (pontosabban a program listája), amellyel a számítást végeztük, az I. Mellékletben olvasható Az elméleti összefüggések és a lista alapján a program futásképes

formában, SCILAB környezetben újraírható és futtatható, így a bemutatott eredmények ellenőrizhetők. V.3 Kollokáció egyetlen bázis-függvénnyel, pontválasztás Vegyük figyelembe, hogy a pontos megoldáshoz, a választott bázisfüggvények közül – szemrevételezés alapján, erre most módunk van – az első (φ1 ) van a legközelebb. Ezt általában nyilvánvalóan nem tehetjük meg; itt is csak azért lehetséges ez az összevetés, mert ismerjük a pontos megoldást (is). A levonható tanulság miatt hagyjuk el a második, harmadik, negyedik és ötödik bázisfüggvényt. Illetve válasszuk kollokációs pontnak a (05,05) pontot Mivel a ∆ψ ( 0.5, 05) = −1 , és ∆φ1 ( 05, 05 ) = −197392 , ezért az α1 = 00507 és a közelítő áramfüggvény: ψ = 0.0507 sin (π x ) sin (π y ) Rajzoltassuk ki az ebben az esetben meghatározható áramvonalakat (V.2 ábra): V.2 ábra – Kollokációs megoldás, közelítő áramvonalak, egy

bázisfüggvény 36 V. A kollokációs módszer Megállapítható, hogy ezek az áramvonalak az V.1 ábrán látható áramvonalaknál jobban hasonlítanak a III4 ábrán látható, pontos áramvonalakhoz Az eltérés (hiba) függvény normája ε = 0.29055 , illetve a pontos és a közelítő megoldás távolsága: ψ −ψ ≅ 0.008135 Mindkét norma jelentősen kisebb, mint az öt bázisfüggvényes kollokációs megoldással kapott normák – ez számszerűen is jelzi, hogy az egy együtthatós megoldás itt jobb. A tanulság: ha mód van rá, akkor a feladatunk megoldása előtt, a gépies megoldás helyett érdemes a rendelkezésünkre álló, műszaki-fizikai képet (illetve az esetleg további ismert, elvárható sajátosságokat) figyelembe venni. Térjünk vissza a kiinduláskor definiált öt-dimenziós függvénytérhez és válasszuk most az alábbi – első ránézésre talán jónak tűnő – kollokációs pontokat. A kiválasztott öt pontot az V.3 ábrán

tüntettük fel, a pontok koordinátái: ( x1 , y1 ) = (0.5, 05) ; ( x2 , y2 ) = (0.33, 033) ; ( x3 , y3 ) = (0.17, 017) ; ( x4 , y4 ) = (0.5, 017) ; ( x5 , y5 ) = (0.5, 033) V.3 ábra – Kollokációs pontok A program futtatásakor az alábbi üzenetet kaptuk: Warning : matrix is close to singular or badly scaled. rcond = 22010D-18 computing least squares solution. Ez azt jelenti, hogy a számítás ugyan ad valamilyen eredményt, de ez az eredmény teljesen használhatatlan! Az viszont megnyugtató, hogy a SCILAB ezt a problémát jelzi, illetve ennek nyomán a pontválasztást legalább addig módosíthatjuk, amíg ez a hibajelzés el nem tűnik. Végeredményben megállapíthatjuk, hogy a kollokációs módszer – a példánk szerint – nagyon érzékeny a megfelelő pont-választásra. Ezért e módszer alkalmazása csak kellő óvatossággal ajánlható. 37 VI. A legkisebb négyzetek módszere VI. A legkisebb négyzetek módszere (LKN) VI.1 A legkisebb négyzetek

módszerének leírása A legkisebb négyzetek módszere – mint az elnevezése is mutatja – egy igazán ígéretes eljárás. A módszernek több verziója is létezik, mi a következő változatot vizsgáljuk S ( ε , ε ) = F , F ⇒ min! (6.1) Az F funkcionál minimumát kell tehát meghatároznunk. A minimum szükséges feltételét a ∂F = 0, i = 1, 2 n egyenletrendszer adja meg. (Bizonyít∂α i ható, hogy ez az egyenletrendszer valóban az F funkcionál minimumát adja.) VI.2 A legkisebb négyzetek módszerének alkalmazása Határozzuk meg a minta-feladat közelítő megoldását (a közelítő áramfüggvényt) a legkisebb négyzetek módszerével! Azért, hogy a megoldás áttekinthetőbb legyen, vezessük be a következő, konkrét számokat jelentő összefüggéseket, illetve jelöléseket:  1  Bi , j := ∫  ∫ ∆φi ∆φ j dx  dy;  0 0 1  1  Bi ,ω := ∫  ∫ ∆φi ω dx  dy;  0 0 1 ahol i, j

= 1, 2⋯ 5; 1 1   Bω ,ω := ∫  ∫ ω 2 dx  dy; 0 0  (6.2) Írjuk fel részletesen a (6.1) funkcionált a (62)-ben definiált állandókkal: F = α12 B1,1 + α 22 B2,2 + α 32 B3,3 + α 42 B4,4 + α 52 B5,5 + Bω ,ω + + 2α1α 2 B1,2 + 2α1α 3 B1,3 + 2α1α 4 B1,4 + 2α1α 5 B1,5 + 2α1 B1,ω + + 2α 2α 3 B2,3 + 2α 2α 4 B2,4 + 2α 2α 5 B2,5 + 2α 2 B2,ω + + 2α 3α 4 B3,4 + 2α 3α 5 B3,5 + 2α 3 B3,ω + 2α 4α 5 B4,5 + 2α 4 B4,ω + 2α 5 B5,ω ; A legkisebb négyzete módszere szerint, az F funkcionál minimumát a ∂F ∂α i = 0, i = 1, 2 n inhomogén, lineáris, algebrai egyenletrendszer szolgáltatja; ebből számítjuk a közelítő megoldásban szereplő α i együtthatókat. 38 VI. A legkisebb négyzetek módszere Használjuk ki, hogy a közelítő megoldás függvényterét kifeszítő bázisfüggvények a Laplace operátor sajátfüggvényei és páronként ortogonálisak. Ezért Bi , j = 0, ha i ≠ j . Emiatt az

együttható mátrix diagonál mátrix lesz, azaz a főátlón kívüli összes eleme zérus Ez pedig azt jelenti, hogy az egyenletrendszer öt, szimultán egyenletre esik szét: B ∂F = 0 ⇒ 2α i Bi ,i + 2 Bi ,ω = 0 ⇒ α i = − i ,ω ; i = 1, 2 5. ∂α i Bi ,i Az iskolapéldában alkalmazott bázisfüggvények előnyös tulajdonságai következtében tehát az α i együtthatók számítása igen egyszerű lett. Ezért érdemes törekedni arra, hogy a bázisfüggvényeink tulajdonságai más esetekben is közelítsék a II.3-ban megfogalmazott tulajdonságokat, mivel ez a megoldás egyszerűsödéséhez vezet. Írjuk ki részletesen az első együttható kiszámítására szolgáló kifejezést: 1 1 1     φ ω dx dy ∫0  ∫0 1  ∫0  ∫0 φ1 ω dx  dy α1 = − = ; 1 1 1 1  2  2   2  2 2  −2π ∫0  ∫0 φ1 dx  dy 2π ∫0  ∫0 φ1 dx  dy ( (

−2π 2 ) 1 (6.3) ) Határozzuk meg az α1 számértékét. Kezdjük a számolást a nevezőben szereplő integrál kiszámításával, de számoljunk kicsit általánosabban: 1 1 1 2  1  x ∫0 sin ( k π x ) dx = − 2kπ cos ( k π x ) sin ( k π x ) 0 +  2  0 = 0.5, k = 1,3, 5, 7,11; Ezek szerint: 1 2  2    k π y sin ( )  ∫0   ∫ sin ( k π x )  dx  dy = 0.25, k = 1, 3,5, 7,11;  0  1 Az α i együtthatók számítására szolgáló kifejezések nevezőjében, (4.1) figyelembe vételével, a 025-öt az i=1 esetben 2π 2 -tel, az i=2 esetben 18π 2 -tel, az i=3 esetben 50π 2 -tel, az i=4 esetben 98π 2 -tel és az i=5 esetben 242π 2 -tel kell megszorozni. 39 VI. A legkisebb négyzetek módszere Ennek alapján írjuk ki részletesen az egyes együtthatók számítására szolgáló kifejezéseket: 1   ∫0  ∫0 φ1 ω dx  dy ; α1

= 2π 2 0.25 1 1   dx φ ω   dy 3 ∫0 ∫0  ; α3 = 2 50π 0.25 1 1   ∫0  ∫0 φ5 ω dx  dy ; α4 = 242π 2 0.25 1 1   ∫0 ∫0 φ2 ω dx  dy ; α2 = 18π 2 0.25 1 1   dx φ ω   dy 4 ∫0  ∫0  ; α4 = 2 98π 0.25 1 (6.4) 1   A (6.4) számlálóiban található Si = ∫  ∫ φi ω dx  dy kifejezés zárt alakban is ki0 0  1 számítható lenne, azonban soktagú, bonyolult kifejezést kapnánk. Integráljunk inkább numerikusan. Erre használható az alábbi SCILAB program: - változók törlése; - a G kitűzése; clear X=[0,0;1,1;1,0]; Y=[0,0;0,1;1,1]; function [fi1] = fiegy omega(x,y); fi1=sin(%pi*x)sin(%piy)(2(y(1-y)+x(1-x))); endfunction - az integrálandó függvény; format(e,12); S1=int2d(X,Y,fiegy omega); disp(S1); alfa1=S1/(2*%pi%pi0.25); disp(alfa1); - kiíratási számformátum; - integrálás; - α1

kiíratása. Amint az a későbbiekből kiderül, az α 4 és különösen az α 5 igen kis értékű, ezért az α i együtthatók számolására YABASIC programot (II. Melléklet) írtunk Az integrálásra az egyszerűség érdekében a „nyers erő” módszerét alkalmaztuk A G tartományt a VI1 ábrán látható módon ráccsal fedtük le, és meghatároztuk minden kis négyzet átlóinak metszéspontját: ∆x ∆y + p ∆x és yq = + q∆y, ahol: 2 2 p = 0,1, 2 (n − 1) és q = 0,1, 2 (n − 1) illetve: ∆x = ∆y = 1 n ; xp = Vezessük be a gi , p , q = φi ( x p , yq ) ω ( x p , yq ) mennyiséget – ez az i-edik bázis- függvény és az örvényesség értékének szorzata az ( x p , yq ) pontban. 40 VI. A legkisebb négyzetek módszere Ezzel a (6.4) kifejezésbeli számlálók (integrálok) közelítő számértéke:  n−1  n −1  Si ≅  ∑  ∑ gi , p , q   ∆x∆y .  p=0 q=0     Az Si számértékek

ismeretében a keresett együtthatók már meghatározhatók: VI.1 ábra – Tartomány felosztása Si αi S1 ≅ 0.328511 α1 ≅ 0.06657 S 2 ≅ 0.004056 α 2 ≅ 9.1317 ⋅10−5 S3 ≅ 0.000526 α 3 ≅ 4.2605 ⋅10−6 S 4 ≅ 0.000137 α 4 ≅ 5.6584 ⋅10−7 S5 ≅ 2.244 ⋅10−5 α 5 ≅ 3.7577 ⋅10−8 Az együtthatók értékeire pillantva megállapítható, hogy az α1 domináns, a további együtthatók legalább három nagyságrenddel kisebbek. A kollokációs módszernél adódó csökkenésnél jelentősen nagyobb a különbség! Amint azt már korábban említettük, az α 4 és különösen az α 5 valóban igen kis számértékkel bír. Az alábbi táblázatban, a második sorban az egy együtthatós közelítés szerepel; a további sorokban rendre egy-egy együtthatóval többet vettünk figyelembe. A teljes közelítő áramfüggvény a hatodik sorban látható: ε A közelítő áramfüggvény ψɶ ψɶ ψɶ ψɶ ψɶ = α1φ1 ; 0.2447511

0.2446168 0.2446146 0.2446144 0.2446144 = α1φ1 + α 2φ2 ; = α1φ1 + α 2φ2 + α 3φ3 ; = α1φ1 + α 2φ2 + α 3φ3 + α 4φ4 ; = α1φ1 + α 2φ2 + α 3φ3 + α 4φ4 + α 5φ5 ; 41 ψ −ψ 0.0017913 0.001790718 0.001790717 0.0017907165 0.0017907165 VI. A legkisebb négyzetek módszere A legfontosabb megállapítás, hogy közelítő megoldást döntően az első bázisfüggvény alkalmazásával kapjuk. A második bázisfüggvény bevonása kis mértékben javítja a megoldást. A harmadik bázisfüggvény bevonása már csak nagyon kis javulást hoz. Ha a negyedik bázis-függvénnyel is számolunk, akkor látunk még egy rendkívül kis pontosítást Az ötödik bázisfüggvény bevonása azonban – az általunk alkalmazott számítási pontosság mellett - már nem hoz kimutatható javulást – ennek alkalmazása tehát lényegében felesleges. Ezek szerint, a módszer ilyen formában történő alkalmazásával a közelítő áramfüggvény és a pontos

áramfüggvény távolsága nem csökkenthető egy korlátérték alá – megmaradó hibával kell számolnunk. VI.2 ábra – LKN megoldás, közelítő áramvonalak A VI.2 ábrán látható áramkép alig különbözik a III4 ábrán látható, pontos áramképtől. A számszerűséget illetően pedig megjegyezzük, hogy a közelítő áramfüggvény és a pontos áramfüggvény kiválasztott áramvonalainak számértéke kicsit kevesebb, mint 7%-kal tér el egymástól. A közelítés eredményei tehát, megfelelő óvatossággal, de számszerűen is használhatók. 42 VII. A Galjorkin eljárás VII. A Galjorkin eljárás VII.1 A Galjorkin eljárás leírása A Galjorkin eljárás elve a következő: S ( ε , ϕi ) = 0, i = 1, 2 n (7.1) A (7.1) formálisan n darab egyenletet jelent, n ismeretlenre Az egyenletekben a ϕi = ϕi ( x ) , x ∈ G ∪ ∂G az általunk definiált bázisfüggvény rendszer A bázisfüggvények választásához néhány szempontot a II3 pontban

adtunk meg A (7.1)-ben megadott skalár-szorzatok nulla értéke azt jelenti, hogy az ε = ε ( x ) eltérés (hiba) függvény minden bázisfüggvényre merőleges (ortogonális) lesz. Másképpen fogalmazva ez azt jelenti, hogy az eltérés (hiba) függvény vetülete az egyes bázisfüggvényekre nulla Ezek szerint abban a függvénytérben, ahol a közelítő megoldást keressük, az eltérés (hiba) függvény nullának látszik. A fentiek alapján kimondhatjuk – az I.10 pontban leírtak értelemszerű kiterjesztésével – hogy az ε = ε ( x ) eltérés (hiba) függvény minimalizáló elem (minimalizáló függvény) lesz VII.2 A Galjorkin eljárás alkalmazása Az eltérés (hiba) függvény konkrét, minta-feladatra vonatkozó alakját a IV.2 pont (4.2) egyenletével adtuk meg Ebben kihasználtuk, hogy az általunk választott bázisfüggvények a Laplace operátor sajátfüggvényei Számítsuk ki először a (7.1) szerinti skalár szorzatot az első bázisfüggvény

felhasználásával: S ( ε , ϕi ) = 0, írjuk fel a skalár szorzat kifejezését részletesen: 1  2 2 2 2 2 ∫0  ∫0 ω − α1 2π φ1 − α 218π φ2 − α 3 50π φ3 − α 4 78π φ4 − α5 242π φ5 φ1 dx  dy = 0 ; 1 ( ) 43 VII. A Galjorkin eljárás A bázisfüggvények ortogonálisak, ezért írható, hogy:  1  φ ω dx   dy 1 ∫0  ∫0  ⇒ α1 = ; 1 1     2π 2 ∫  ∫ φ12 dx  dy 0 0  1 1  2 ∫0  ∫0 ω − α1 2π φ1 φ1 dx  dy = 0 1 ( ) Az eredmény ugyanaz, mint a (6.3) kifejezés végeredménye! Végezzük el a skaláris szorzást rendre a φ2 , φ3 , φ4 és φ5 bázis-függvénnyel is, akkor végeredményként azt kapjuk, hogy az α i együtthatók azonosak a (6.4) kifejezéssel adott együtthatókkal. Ezek szerint – az általunk vizsgált, rendkívül kedvező példában – a legkisebb négyzetek módszerével kapott eredménnyel azonos

eredményre jutunk. A Galjorkin eljárás rendkívül fontos, mert igen széles körben (pl. számos CFD algoritmusban is) alkalmazzák. A közelítő megoldást meghatározó együtthatókat – lineáris feladat esetében – általában inhomogén, lineáris, algebrai egyenletrendszerből számoljuk. A mi, igen előnyös minta-feladatunkban, az LKN módszer alkalmazásánál és a Galjorkin eljárásnál ez az egyenletrendszer szimultán egyenletekre esett szét – az egyes együtthatókat külön-külön számíthattuk. Ez általában nincs így, de törekedni lehet (kell) arra, hogy a kapott egyenletrendszer együttható mátrixa ritka legyen (lehetőleg sok eleme legyen zérus - például tridiagonális vagy pentadiagonális mátrixot kapjunk). 44 Irodalomjegyzék Irodalomjegyzék [1] BREBBIA, C.A – CONNOR, JJ: Finite Element Techniques for Fluid Flow, Newness-Butterworths, London, 1976. [2] CAMPBELL S. – CHANCELIER, J-P – NIKOUKHAH, R: Modeling and Simulation in

Sciilab/Scicos, Springer, 2006. [3] GAUSZ, T.: Áramlástan, wwwdoksihu, 2012 [4] GÖRKE, L.: Halmazok, relációk, függvények, Tankönyvkiadó, 1969 [5] HOFFMANN, A.: Numerikus módszerek az áramlástan és a hőtan parciális differenciálegyenleteinek megoldására I, Tankönyvkiadó, 1978 [6] VAN KAN, J. – SEGAL, A: Numerik Partieller Differentialgleichungen für Ingenieure, B.G Teuner, Stuttgart, 1995 [7] [8] [9] MÁTÉ, L.: Funkcionálanalízis műszakiaknak, Műszaki Könyvkiadó, 1976 PETZ, D.: Bevezetés a lineáris analízisbe és az alkalmazásaiba, BME 2013 POPPER, GY. – KURUTZNÉ DR KOVÁCS MÁRTA: SZÁMÍTÁSTECHNIKA – A végeselem módszer matematikai alapjai. BME-MTI M319, Budapest, 1982 KARÁTSON, J. – HORVÁTH, R – IZSÁK, F: Parciális differenciálegyenletek numerikus módszerei számítógépes alkalmazásokkal, 2013, TÁMOP-412A/1-11/1 MSc Tananyagfejlesztés www.ybasicde – A „yabasic” nyelv honlapja [10] [11] 45 A

kollokációs módszer program-listája I. Melléklet: Kollokációs módszer, program-lista Az itt következő program (túl) egyszerű, de a célunk az, hogy könnyen érthető lépésekből épüljön fel. Messze nem használjuk ki a SCILAB környezet nyújtotta lehetőségeket Ha Tisztelt Olvasó úgy óhajtja, akkor alakítsa át a programot sokkal tömörebb, elegánsabb formába. A programot a SCILAB 5.51 verziójában fejlesztettük Jelenleg már a 61 verzió is elérhető – ebben is fut a program A Tisztelt Olvasó, ha óhajtja, futtassa a programot és tanulmányozza azokat az eredményeket is, amelyeket a program kijelez, de az V. pontban külön nem részleteztünk! 46 A kollokációs módszer program-listája clear disp ("----------------------------------------------------------------"); disp ("--K O L L O K A C I O S megoldas ---"); disp ("----------------------------------------------------------------");

elvalaszto="=====-------------------========-----------------------=====" //A megoldas-fuggveny es az orvenyesseg-fuggveny ------------------------function [psi]=megoldas(x, y); psi=x*(1-x)y(1-y); endfunction function [omegafv]=omega(x, y); omegafv=-2*(y(1-y)+x(1-x)); endfunction //Kovetkeznek a bazis-fuggvenyek ----------------------------------------function [fi1]=fiegy(x, y); fi1=sin(%pi*x)sin(%piy); endfunction function [fi2]=fiket(x, y); fi2=sin(%pi*x3)sin(%piy3); endfunction function [fi3]=fihar(x, y); fi3=sin(%pi*x5)sin(%piy5); endfunction function [fi4]=fineg(x, y); fi4=sin(%pi*x7)sin(%piy7); endfunction function [fi5]=fiot(x, y); fi5=sin(%pi*x11)sin(%piy11); endfunction //Kovetkeznek a bazis-fuggvenyek Laplace transzformaltjai ---------------function [Lfi1]=lapfiegy(x, y); Lfi1=-2*%pi%pisin(%pix)sin(%piy); endfunction function [Lfi2]=lapfiket(x, y); Lfi2=-18*%pi%pisin(%pix3)sin(%piy3); endfunction function [Lfi3]=lapfihar(x, y); Lfi3=-50*%pi%pisin(%pix5)sin(%piy5);

endfunction function [Lfi4]=lapfineg(x, y); Lfi4=-98*%pi%pisin(%pix7)sin(%piy7); endfunction function [Lfi5]=lapfiot(x, y); Lfi5=-242*%pi%pisin(%pix11)sin(%piy11); endfunction //Felosztjuk a G tartományt es numerikusan kiszamoljuk a megoldas-fuggvenyt -n=100; lep=1/n; xx=[0:lep:1]; yy=[0:lep:1]; for i=1:n+1; for j=1:n+1; zm(i,j)=megoldas(xx(i),yy(j)); end; end xset(window,0); xname(A tenyleges megoldas - 0. ablak); xtitle("A zart alaku megolda - az aramfuggveny"); plot3d1(xx,yy,zm); //Kirajzoltatjuk a megoldas-fuggvenyt --------------- 47 A kollokációs módszer program-listája xset(window,1); xname(A tenyleges megoldas - aramvonalak - 1. ablak); xtitle("A zart alaku megoldas - az aramvonalak"); contour(xx,yy,zm,5); //Aramvonalak a pontos megodas alapjan -----------//Megadjuk a kollokacios pontokat -------------------------------------xk(1)=0.25; yk(1)=013; xk(2)=0.58; yk(2)=096; xk(3)=0.67; yk(3)=061; xk(4)=0.86; yk(4)=054; xk(5)=0.95; yk(5)=043; //Az

egyutthato-matrix elemeinek szamitasa ----------------------------a=zeros(5,5); for j=1:5 a(j,1)=lapfiegy(xk(j),yk(j)); a(j,2)=lapfiket(xk(j),yk(j)); a(j,3)=lapfihar(xk(j),yk(j)); a(j,4)=lapfineg(xk(j),yk(j)); a(j,5)=lapfiot(xk(j),yk(j)); end disp ("Az egyutthato matrix:"); disp (a); disp (elvalaszto); //Az egyutthato matrix ellenorzese ------------------------------------b=inv(a); disp ("Az egyutthato matrix es inverzenek szorzata:"); disp (a*b); //A jobboldal es az alfa egyutthatok szamitasa ------------------------for i=1:5; jobboldal(i)=omega(xk(i),yk(i)); end alfa=ajobboldal; disp (elvalaszto); disp ("A kozelito megoldas egyutthatoi"); disp (alfa); //Az elteres szamitasa a kollokacios pontokban ------------------------for i=1:5 sv(i)=alfa(1)*lapfiegy(xk(i),yk(i))+alfa(2)lapfiket(xk(i),yk(i))+alfa(3)lapfihar(xk(i),yk(i))+alfa(4)lapfineg(xk(i),yk(i))+alfa(5)lapfiot(xk(i),yk(i)); sq(i)=omega(xk(i),yk(i))-sv(i); end disp (elvalaszto); disp ("A

kollokaciok erteke (kb. nullanak kell lennie):"); disp (sq); 48 A kollokációs módszer program-listája //Az orvenyesseg es a kozelito megoldasbol kapott orvenyesseg elterese ----elteres negyzetosszeg=0; for i=1:n+1 for j=1:n+1 kozelit(i,j)=alfa(1)*lapfiegy(xx(i),yy(j))+alfa(2)lapfiket(xx(i),yy(j))+alfa(3)lapfihar(xx(i),yy(j))+alfa(4)lapfineg(xx(i),yy(j))+alfa(5)lapfiot(xx(i),yy(j)); elteres(i,j)=omega(xx(i),yy(j))-kozelit(i,j); elteres negyzetosszeg=elteres negyzetosszeg+elteres(i,j)*elteres(i,j); end end epszilon epszilon=elteres negyzetosszeg*leplep; disp (elvalaszto); disp ("Az orvenyesseg es a kozelitesenek kulonbsege <kozelito> skalar szorzata:"); disp (epszilon epszilon); disp ("es a normaja:"); disp (sqrt(epszilon epszilon)); xset(window,2); xname(Az elteres fuggveny - 2. ablak); xtitle("LapHiba-fuggveny"); plot3d1(xx,yy,elteres); //Az elteres abrazolasa //A kozelito aramfuggveny szamitasa es a pontos aramfuggvenytol valo

elterese -----------------megoldas elteres negyzet=0; for i=1:n+1 for j=1:n+1 arfvkozelit(i,j)=alfa(1)*fiegy(xx(i),yy(j))+alfa(2)fiket(xx(i),yy(j))+alfa(3)fihar(xx(i),yy(j))+alfa(4)fineg(xx(i),yy(j))+alfa(5)fiot(xx(i),yy(j)); arfvelteres(i,j)=megoldas(xx(i),yy(j))-arfvkozelit(i,j); megoldas elteres negyzet=megoldas elteres negyzet+arfvelteres(i,j)*arfvelteres(i,j); end end megoldas elteres negyzet=megoldas elteres negyzet*leplep; disp ("Az aramfuggveny es a kozelitese kulonbsegenek <kozelito> skalarszorzata:"); disp (megoldas elteres negyzet); disp ("es a normaja:"); disp (sqrt(megoldas elteres negyzet)); xset(window,3); xname(Az elteres fuggveny - 3. ablak); xtitle("Hiba-fuggveny"); plot3d1(xx,yy,arfvelteres); //Az aramfuggveny elteres abrazolasa 49 A kollokációs módszer program-listája xset(window,4); xname(A kozelito aramvonalak - 4. ablak); xtitle("Kozelito aramvonalak"); contour(xx,yy,arfvkozelit,5); //A kozelito aramvonalak

kepe disp (elvalaszto); disp ("Az aramfv maximuma"); disp (megoldas(0.5,05)); disp ("Az kozelites maximuma"); disp (arfvkozelit(50,50)); disp ("*===== P R O G R A M V E G E ====="); 50 LKN program-listája II. Melléklet: Legkisebb négyzetek módszere, program-lista Az itt következő, YABASIC program egyszerű felépítésű, véleményünk szerint külön magyarázatot nem igényel. n=3333 : rem ez tobb, mint 11 millio pontot jelent dim x(n) : dim y(n) delta x=1/n : delta y=1/n for i=0 to n-1 x(i)=i/n+delta x/2 : y(i)=i/n+delta y/2 next i suq=0 for i=0 to n-1 for j=0 to n-1 suq=suq+sin(pi*x(i))sin(piy(j))2(y(j)(1-y(j))+x(i)(1-x(i))) next j next i S1=suq/n/n : print " S1 = ",S1 suq=0 for i=0 to n-1 for j=0 to n-1 suq=suq+sin(3*pix(i))sin(3piy(j))2(y(j)(1-y(j))+x(i)(1-x(i))) next j next i S2=suq/n/n : print " S2 = ",S2 suq=0 for i=0 to n-1 for j=0 to n-1 suq=suq+sin(5*pix(i))sin(5piy(j))2(y(j)(1-y(j))+x(i)(1-x(i))) next j next i

S3=suq/n/n : print " S3 = ",S3 51 LKN program-listája suq=0 for i=0 to n-1 for j=0 to n-1 suq=suq+sin(7*pix(i))sin(7piy(j))2(y(j)(1-y(j))+x(i)(1-x(i))) next j next i S4=suq/n/n : print " S4 = ",S4 suq=0 for i=0 to n-1 for j=0 to n-1 suq=suq+sin(11*pix(i))sin(11piy(j))2(y(j)(1-y(j))+x(i)(1-x(i))) next j next i S5=suq/n/n : print " S5 = ",S5 print : print "----------------------------------------------------------" : print alfa1=S1/(0.25*2pipi) : print "alfa1 = ",alfa1 alfa2=S2/(0.25*18pipi): print "alfa2 = ",alfa2 alfa3=S3/(0.25*50pipi): print "alfa3 = ",alfa3 alfa4=S4/(0.25*98pipi): print "alfa4 = ",alfa4 alfa5=S5/(0.25*242pipi): print "alfa5 = ",alfa5 print : print "----------------------------------------------------------" :print end 52