Szakdolgozat. Írta: Újvári Zsuzsanna. Matematika BSc Alkalmazott matematika szakirány

Kétváltozós integrálok numerikus közelítése a peremre való transzformálással Szakdolgozat Írta: Újvári Zsuzsanna Matematika BSc Alkalmazott matematika szakirány

Témavezet˝o: Dr. Gáspár Csaba egyetemi tanár Numerikus Analízis Tanszék

Eötvös Loránd Tudományegyetem Természettudományi Kar 2014

Köszönetnyilvánítás Els˝osorban szeretném megköszönni témavezet˝omnek, Dr. Gáspár Csabának, hogy felkeltette érdekl˝odésemet a téma iránt, hogy bizalommal fordulhattam hozzá kérdéseimmel, problémáimmal, s tanácsaival nagymértékben segítette a dolgozat elkészítését. Szeretném megköszönni továbbá Lócsi Leventének is a témámmal kapcsolatos programok, Matlabban történ˝o megírásában nyújtott segítségét, hasznos tanácsait, iránymutatását. Köszönettel tartozom családomnak, barátaimnak, hogy végig mellettem álltak, bíztattak és támogattak.

1

Tartalomjegyzék Bevezetés

5

1. Többváltozós függvények néhány integrálási módszere

6

1.1. Visszavezetés egyváltozós függvények integrálására . . . . . . . . . . . .

6

1.2. Az integrálási tartomány approximációja . . . . . . . . . . . . . . . . . .

7

2. Közelít˝o integrálás a tartomány peremére transzformálással

10

2.1. Szórt alappontú interpoláció a radiális bázisfüggvények módszerével . . .

10

2.2. A peremre való transzformálás módszere . . . . . . . . . . . . . . . . . .

11

3. Néhány numerikus példa

14

3.1. Integrálás négyzet tartományon . . . . . . . . . . . . . . . . . . . . . . .

14

3.2. Integrálás kör tartományon . . . . . . . . . . . . . . . . . . . . . . . . .

21

3.3. Integrálás am˝obaszer˝u tartományon . . . . . . . . . . . . . . . . . . . .

26

Összefoglalás

35

Irodalomjegyzék

36

Függelék

37

2

Ábrák jegyzéke 1.1. Tartomány közelítése téglalapokkal . . . . . . . . . . . . . . . . . . . . .

8

1.2. Tartomány közelítése háromszögekkel . . . . . . . . . . . . . . . . . . .

8

3.1. n db pont generálása a négyzeten . . . . . . . . . . . . . . . . . . . . . .

15

3.2. TPS-interpolációs függvény . . . . . . . . . . . . . . . . . . . . . . . .

15

3.3. A tesztfüggvény és a TPS-interpolációs függvény eltérése . . . . . . . . .

16

3.4. A tesztfüggvény és a TPS-interpolációs függvény eltérésének nagyítása .

16

3.5. A hiba értéke közelít a 0-hoz . . . . . . . . . . . . . . . . . . . . . . . .

17

3.6. A négyzet peremén felvett pontok . . . . . . . . . . . . . . . . . . . . .

18

3.7. A négyzet peremén felvett pontok növelése . . . . . . . . . . . . . . . .

20

3.8. n db pont generálása a körön belül . . . . . . . . . . . . . . . . . . . . .

21

3.9. TPS-interpolációs függvény . . . . . . . . . . . . . . . . . . . . . . . .

22


22

3.11. A hiba értéke közelít a 0-hoz . . . . . . . . . . . . . . . . . . . . . . . .

23

3.12. A kör peremén felvett pontok . . . . . . . . . . . . . . . . . . . . . . . .

24

3.13. A kör peremén felvett pontok növelése . . . . . . . . . . . . . . . . . . .

25

3.14. Am˝oba kirajzolása . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

26

3.15. Felbontás 10 ponttal . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

27

3.16. Felbontás 50 ponttal . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

27

3.17. Felbontás 150 ponttal . . . . . . . . . . . . . . . . . . . . . . . . . . . .

28

3.18. n db pont generálása az am˝obaszer˝u tartományon . . . . . . . . . . . . .

29

3.19. Tesztfüggvény . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

29

3.20. TPS-interpolációs függvény (1) . . . . . . . . . . . . . . . . . . . . . . .

30

3.21. TPS-interpolációs függvény (2) . . . . . . . . . . . . . . . . . . . . . . .

30


31

3.23. A tesztfüggvény és a TPS-interpolációs függvény eltérésének nagyítása .

31

3.24. Az am˝oba peremén felvett pontok . . . . . . . . . . . . . . . . . . . . .

32

3.25. Az am˝oba perempontjainak irányvektorai . . . . . . . . . . . . . . . . .

33

3

Táblázatok jegyzéke 3.1. Az alappontok és peremmenti pontok növekedésével kapott hiba . . . . .

17

3.2. Az integrál értékének változása a véletlen pontok növekedésével . . . . .

19

3.3. Az integrál értékének változása a perempontok növekedésével . . . . . .

20

3.4. Az alappontok és peremmenti pontok növekedésével kapott hiba . . . . .

23


24


25

3.7. Integrál értéke Riemann-integrálközelít˝o összeggel . . . . . . . . . . . .

28

3.8. Az alappontok és peremmenti pontok növekedésével kapott hiba . . . . .

32


34


34

4

Bevezetés A szakdolgozatom egy olyan módszert mutat be, amely lehet˝ové teszi, hogy a bonyolult tartományokon vett kétváltozós integrálokat is ki tudjuk integrálni, ami egyáltalán nem könny˝u. Itt egy olyan közelít˝o eljárást alkalmazunk, amihez nem kell a tartomány semmiféle sokszögfelbontása, hanem egy kétváltozós függvénynek a tartományon vett integrálját visszavezetjük egy alkalmas másik függvénynek a tartomány peremén vett integráljára, ami már lényegében egyváltozós, így az ismert kvadratúraformulák - értelemszer˝u változtatásokkal - már nehézség nélkül alkalmazhatók.

5

1. fejezet Többváltozós függvények néhány integrálási módszere Legyen egy Ω tartomány feletti integrál: ZZ I(f, Ω) =

f (x, y)dxdy Ω

Ennek egy közelítése: In (f, Ω) :=

n X

αi f (xi , yi ),

i=1

ahol az αi számok súlyok, az (x1 , y1 ), . . . , (xn , yn ) ∈ Ω pontok pedig integrálási alappontok. Egy ilyen összeget nevezünk kubatúraképletnek. Az egydimenziós esethez képest az a nagy különbség, hogy az integrációs tartományt is közelíteni kell. [2]

1.1.

Visszavezetés egyváltozós függvények integrálására

Legyen Ω ⊂ R2 zárt korlátos tartomány és f (x, y) folytonos ebben a tartományban. Legegyszer˝ubben úgy integrálhatunk, ha az alaptartomány olyan, hogy a következ˝o alakban átírhatjuk: ZZ Ω

Zb Zd(x) Zb f (x, y)dxdy = f (x, y)dxdy = g(x)dx, a c(x)

a

ahol c, d az x ∈ [a, b] folytonos függvényei, ekkor Zd(x) g(x) := f (x, y)dy c(x)

6

(1.1)

Adott a következ˝o képlet: Im (F ; 0, 1) :=

m X

(m) (m) αi F (ti )

Z1 ≈

i=1

F (t)dt 0

Ezt átírva g függvényre és [a, b] intervallumra: Zb g(x)dx ≈

m X

(m)

(m)

ai g(xi ) =: Im (g; a, b)

i=1

a

ebb˝ol (m) g(xi )

≈

n X

(m)

(n)

bij f (xi , yj )

j=1

Ezeket kombinálva kapjuk a következ˝ot: I(f, Ω) ≈ Im (g; a, b) =

m X

(m) (m) ai g(xi )

≈

m X

i=1

i=1

(m) ai

n X

(m)

(n)

bij f (xi , yj ) := Im,n (f, Ω)

j=1

Ha a tartomány a Q := {a ≤ x ≤ b, c ≤ y ≤ d} téglalap, az Im,n képlet: Zb Zd f (x, y)dxdy ≈ Im,n (f )Im,n (f ) := |Q|

I(f, Ω) := a

m X m X

(m)

(n)

(m)

(n)

f (xi , yj )αi αj ,

i=1 j=1

c

ahol |Q| := (b − a)(d − c). Az ilyen képletet tenzorszorzat formulának nevezzük. [1]

1.2.

Az integrálási tartomány approximációja

Ha a tartomány általánosabb és így a kiszámítandó integrál nem írható fel (1.1) alakban, akkor a kubatúraképlet számításához a tartományt fel kell darabolni és ezeket kell közelíteni. • Az egyik lehet˝oség, amikor egy téglalapos ráccsal borítjuk le a tartományt, h és k lépéstávolságra x és y irányban, és az eredeti Ω tartomány peremét a rács egyeneseivel közelítjük. Ezzel kapunk egy Ω-t közelít˝o Ω’ tartományt, melynek minden téglalapján valamilyen tenzorszorzat integrációt alkalmazhatunk.

7

1.1. ábra. Tartomány közelítése téglalapokkal

Ez egy durvább közelítés, amely járható lenne, de nem el˝onyös, mivel a peremen található téglalapok már csonkítottak, amelyeken nem alkalmazható a formula. • A másik lehet˝oség már finomabb és pontosabb az az, ha az Ω tartományt háromszögekre bontjuk fel.

1.2. ábra. Tartomány közelítése háromszögekkel

Jelöljük a háromszögeket ∆i -vel és számuk legyen n Ω0 :=

n [

∆i

i=1

Ezután minden ∆i felett lineárisan interpoláljuk az f függvényt. ∆i három csúcsában f vegye fel az fi1 , fi2 , fi3 értékeket. Ha ∆i nem elfajult háromszög, akkor ezen 3 érték egyértelm˝uen definiálja az erre illeszked˝o síkot. Ennek egyenlete legyen z = Fi (x, y). A sík és az x, y sík közötti hasáb I3 (∆i , f ) térfogata: 1 I3 (∆i , f ) = (fi1 + fi2 + fi3 )|∆i |, 3 8

(1.2)

ahol |∆i | a ∆i háromszög területe. A hasábok összeadásával kapjuk a következ˝o kubatúraképletet: n

1X I(Ω , f ) := (fi1 + fi2 + fi3 )|∆i | 3 i=1 0

Ezek után a kövezkez˝o közelítést alkalmazzuk: ZZ ZZ n ZZ X f (x, y)dxdy ≈ Fi (x, y)dxdy =: In (Ω0 , f ). f (x, y)dxdy ≈ Ω

i=1 ∆ i

Ω0

Az (1.2) képlet a trapézszabály általánosítása. Megemlíthetjük a középpont szabály általánosítását, melyet a következ˝o képlet ad meg: In (f ) =

n X

f (si , ti )|∆i |,

i=1

ahol (si , ti ) a ∆i súlypontja. Megjegyezzük, hogy mindkét képlet pontos az els˝o- illetve másodfokú (kétváltozós) polinomokon. [2]

9

2. fejezet Közelít˝o integrálás a tartomány peremére transzformálással 2.1.

Szórt alappontú interpoláció a radiális bázisfüggvények módszerével

Adottak: P1 , P2 , . . . , Pn ∈ R2 helyek (Pi = (xi , yi ), i = 1, . . . n) és az azokhoz rendelt f1 , f2 , . . . , fn ∈ R értékek (fi = f (Pi ), ahol f : R2 → R adott függvény), ekkor f (x, y) :=

n X

αj Φ(P − Pj ),

j=1

ahol Φ adott radiális bázisfüggvény, azaz a függvény az argumentum euklideszi normájának függvénye Φ(P ) = ϕ(kP k), ahol ϕ egy alkalmas egyváltozós függvény. Az ismeretlen α1 , α2 , . . . , αn együtthatók az interpolációs egyenletrendszerb˝ol számíthatóak ki: n X j=1

αj Φ(Pk − Pj ) =

n X

αj Φ(xk − xj , yk − yj ) = f (Pk ) = fk ,

j=1

A megoldás m˝uveletigénye: O(N 3 ). Az egyenletrendszer átírva: α1 · Φ(P1 − P1 ) + · · · + αn · Φ(P1 − Pn ) = f (P1 ) α1 · Φ(P2 − P1 ) + · · · + αn · Φ(P2 − Pn ) = f (P2 ) .. . α1 · Φ(Pn − P1 ) + · · · + αn · Φ(Pn − Pn ) = f (Pn )

10

k = 1, . . . , n

Másképp: 

Φ11 . . .  . ...  ..  Φn1 . . . {z |

Φ1n Φnn

A



  α1 f (P1 )  .   ..   ..  =  .    αn f (Pn ) }

   

Az A mátrixot interpolációs mátrixnak nevezzük. Egy gyakran alkalmazott radiális bázisfüggvény: Φ(x, y) =

p x2 + y 2 + c2 (MQ),

ahol c ∈ R adott skálázó paraméter. Egy másik népszer˝u radiális bázisfüggvény: Φ(x, y) := (x2 + y 2 ) log

p x2 + y 2 . Az

erre alapuló interpolációs technikát vékony lemez módszernek (thin plate splines) nevezzük. A kés˝obbiekben ezt fogjuk használni. Ekkor a mátrixelemek: ( 0, Φij = p ((xi − xj )2 + (yi − yj )2 ) log (xi − xj )2 + (yi − yj )2 ,

ha i = j ha i 6= j

Az így kapott A mátrix szimmetrikus és a diagonál elemei 0-k lesznek. Innen kapjuk meg az együtthatókat:  α1  .  ..  αn





 f (P1 )    ..  = A−1   .    f (Pn ) | {z } f

Fontos megjegyeznünk, hogy az egyenletrendszer akkor és csakis akkor oldható meg minden f konstans vektor esetén, ha A invertálható. [3]

2.2.

A peremre való transzformálás módszere

Legyen Ψ egy másik olyan radiális bázisfüggvény a Φ mellett, hogy ∆Ψ = Φ, ahol a ∆ a Laplace-operátor. Közvetlen számolással ellen˝orizhet˝o, hogy ha Ψ =

1 16

r4 log r −

Φ. Polárkoordinátákba átírva kapjuk: ∆Ψ(r) = 1r (rΨ0 (r))0 = Φ(r) 4 1 1 Ψ0 (r) = 16 4r3 log r + rr − 2r3 = 16 (4r3 log r − r3 ) rΨ0 (r) =

1 16

(rΨ0 (r))0 =

(4r4 log r − r4 ) 1 16

16r3 log r − 4r4 1r − 4r3 = r3 log r

∆Ψ(r) = 1r (rΨ0 (r))0 = r2 log r = Φ(r)

11

r4 2

, akkor ∆Ψ =

A

n P

αj Φ(Pk − Pj ) = fk megoldása után kapott α1 , . . . , αn együtthatókkal képezzük

j=1

az következ˝o függvényt: u(P ) :=

n P

αj Ψ(P − Pj ).

j=1

Erre a függvényre ∆u ≈ f teljesül, hiszen n n P P ∆u(P ) := αj ∆Ψ(P − Pj ) = αj Φ(P − Pj ) = f (P ). j=1

j=1

Ezt integrálva az adott tartományon azt kapjuk, hogy elég csak a peremen vett vonalintegrált számolni, ami lényegesen leegyszer˝usíti a számolást. [4][5] A peremmenti integrálok numerikus kiszámítása az ismert egyszer˝u kvadratúrafomulák (pl. trapézformula) értelemszer˝u általánosításával történhet: Z

Z

Z f (P )dΩ ≈ Ω

∆u(P )dΩ =

∂u dΓ, ∂n

Γ

Ω

ahol felhasználtuk a Gauss-féle divergencia-tételb˝ol nyomban adódó

R

∆u(P )dΩ =

Ω

R Γ

∂u dΓ ∂n

egyenl˝oséget. ∂u A ∂n kiszámítása az alábbi formulával történik: n

X ∂Ψ ∂u (P ) = αj (P − Pj ), ∂n ∂n j=1

(2.1)

ahol:

∂Ψ ∂Ψ ∂Ψ = · nx + · ny , ∂n ∂x ∂y és az nx és ny az adott pontba húzott kifelé mutató irányvektor koordinátái, azaz n = (nx , ny ). A Ψ függvény (derékszög˝u koordinátákat használva) az alábbi alakú: 1 Ψ(x, y) = 16

2 2 2 p (x + y ) 2 2 2 (x + y ) log (x2 + y 2 ) − 2

A Ψ függvény x szerinti deriváltja most már közvetlenül számolható: ∂Ψ 1 = x(x2 + y 2 )(2 log(x2 + y 2 ) − 1) ∂x 16 A Ψ függvény y szerinti deriváltja hasonlóan adódik: ∂Ψ 1 = y(x2 + y 2 )(2 log(x2 + y 2 ) − 1) ∂y 16 Ezek után (2.1) kifejezést a konkrét αj -k ismeretében már egyszer˝uen ki lehet értékelni.

12

R ∂u ∂u peremen felvett értékeit kiszámítva, a ∂n dΓ peremintegrál közelít˝o kiszámíA ∂n Γ

tása már történhet a trapézformula alkalmazásával: Z Γ

M −1 X ∂u dΓ ≈ ∂n j=0

∂u (Qj ) ∂n

∂u + ∂n (Qj+1 ) · kQj+1 − Qj k, 2

ahol Q0 , Q1 , . . . , QM −1 már a Γ peremen elhelyezett pontok, és QM := Q0 .

13

3. fejezet Néhány numerikus példa 3.1.

Integrálás négyzet tartományon

Legyen adott f (x, y) = sin(πx) sin(πy) függvény és a radiális bázisfüggvény pedig a TPS (Thin Plate Spline), azaz Φ(r) = r2 log r. Az Ω tartomány legyen az egységnégyzet: Ω := [0, 1] × [0, 1]. R Az f (x, y)dxdy integrál pontos értéke egyszer˝uen kiszámítható: Ω



Z1 Z1

Z f (x, y)dxdy =

sin(πx) sin(πy)dxdy =  0

Ω

=

0

2

Z1

sin(πx)dx = 0

1 !2 2 2 1 4 = − cos(πx) = 2 ≈ 0, 4053 π π π 0

Matlab program segítségével kiszámoltuk az f (x, y) függvényünk közelítését az adott radiális bázisfüggvény segítségével. El˝oször is generálunk n (n = 500) pontot a [0, 1] × [0, 1] négyzeten, és úgy, hogy bármely 2 pont távolsága legalább ε legyen (ebben az esetben ε = 0.02).

14

3.1. ábra. n db pont generálása a négyzeten

Itt láthatjuk a bázisfüggvénnyel számított közelítést:

3.2. ábra. TPS-interpolációs függvény

15

Most pedig nézzük meg, hogy mennyi a tesztfüggvényt˝ol való eltérése, azaz a hiba:

3.3. ábra. A tesztfüggvény és a TPS-interpolációs függvény eltérése

A hiba értéke igen kicsi, 0, 67%-körül van. Megfigyelhet˝o, hogy a sarkokban van a legnagyobb eltérés:

3.4. ábra. A tesztfüggvény és a TPS-interpolációs függvény eltérésének nagyítása

16

Nézzük meg a következ˝o táblázatban hogyan változik a hiba értéke a tartományban elhelyezett alappontok (n) és a peremmenti pontok (m) növekedésével:

3.1. táblázat. Az alappontok és peremmenti pontok növekedésével kapott hiba

3.5. ábra. A hiba értéke közelít a 0-hoz

17

A [0, 1] × [0, 1] négyzet peremén vegyünk fel pontokat

3.6. ábra. A négyzet peremén felvett pontok

Ebben az esetben az irányvektorokat könny˝u meghatározni: • a lefelé mutató irányvektor: n = (0, −1) • a jobbra mutató irányvektor: n = (1, 0) • a felfelé mutató irányvektor: n = (0, 1) • a balra mutató irányvektor: n = (−1, 0)

18

A Matlab programmal számolt eredményt a következ˝o táblázat foglalja össze és megmutatja hogyan közelít az integrál számított értéke, ha az alappontok számát növeljük, a peremen felvett pontok száma fix (m = 200):

3.2. táblázat. Az integrál értékének változása a véletlen pontok növekedésével

Megvizsgáltuk továbbá, hogyan közelítünk az integrálhoz, ha a peremen felvett pontok számát növeljük egy fix véletlen pontszám mellett (n = 200):

19

3.7. ábra. A négyzet peremén felvett pontok növelése

3.3. táblázat. Az integrál értékének változása a perempontok növekedésével

Tapasztalataink szerint 320 − 640 véletlen pont között (fix perempont mellett) már a pontos integrálhoz közeli értéket adott a program, és 1280-ra 10 esetben is 0, 4053 lett az 20

eredmény. Minél nagyobb a véletlen pontok és a peremen felvett pontok száma, annál jobb a közelítésünk.

3.2.

Integrálás kör tartományon

p Legyen Ω az egységkör, és f (x, y) = 1 − x2 − y 2 függvény, mely egy félgömböt ír le. R Az f (x, y)dxdy integrál pontos értéke a félgömb térfogata: 23 π ≈ 2, 0944 Ω

Matlab programmal kiszámoljuk az f (x, y) függvényünk közelítését a TPS radiális bázisfüggvény segítségével. Generálunk n (n = 200) számot az egység sugarú körön belül, és megadjuk, hogy bármely 2 pont távolsága legalább ε legyen (itt ε = 0.05).

3.8. ábra. n db pont generálása a körön belül

21

Itt láthatjuk a bázisfüggvénnyel számított közelítést:

3.9. ábra. TPS-interpolációs függvény

Nézzük meg, hogy mennyi itt a hiba, azaz a tesztfüggvényt˝ol való eltérés:


A hiba értéke 0, 74%-körül van. 22

Nézzük meg hogyan változik a hiba értéke az alappontok (n) és a peremmenti pontok (m) növekedésével, amit a következ˝o táblázat szemléltet:


3.11. ábra. A hiba értéke közelít a 0-hoz

23

Az egység sugarú kör peremén felvett pontok

3.12. ábra. A kör peremén felvett pontok

Itt az irányvektorok koordinátái megegyeznek a peremen felvett pontok koordinátáival, azaz a j. pont x és y koordinátája a következ˝o: x(j) = cos j·2π , y(j) = sin j·2π , ahol m m az m a pontok száma. A programmal számolt eredményeket a következ˝o táblázat szemlélteti, hogy miként változik (közelít) az integrál értéke, ha az alappontok számát növelem, és a kör peremén felvett pontok száma m = 50:


24

Valamint megnézzük hogyan változik az integrál értéke, ha a peremen felvett pontok számát növelem s a véletlen pontok száma fix (n = 200):

3.13. ábra. A kör peremén felvett pontok növelése


Összességében az integrál számított értékei, mind a véletlen pontok növekedésekor, mind pedig a peremen felvett pontok növekedésekor nagyon jól közelítenek a pontos integrálhoz, azaz a félgömb térfogátához.

25

3.3.

Integrálás am˝obaszeru˝ tartományon

Itt legyen a tesztfüggvény f (x, y) = cos( x2 ) cos( y2 ) függvény, melynek integrálját közelítjük a már korábban is használt radiális bázisfüggvénnyel: Φ(r) = r2 log r. A tartomány most egy szabálytalan am˝obaszer˝u tartomány. A perem paraméterezése: x(t) = (esin(t) sin2 (2t) + ecos(t) cos2 (2t)) cos(t) t ∈ [0, 2π] y(t) = (esin(t) sin2 (2t) + ecos(t) cos2 (2t)) sin(t) Az am˝oba alakját az alábbi ábra szemlélteti:

3.14. ábra. Am˝oba kirajzolása

A pontos integrál értéke ebben az esetben nem számolható ki, mint korábban, így itt a "pontos" értéket egy Riemann-integrálközelít˝o összeggel közelítjük: Z X f (x, y)dxdy ≈ f (xk , yj ) · h2 , Ω

k,j=1 m (xk ,xj )∈Ω

ahol h a rács lépésköze. A relatív hiba kiszámításakor a pontos integrál-érték helyett ezzel a Reimann-integrálközelít˝o összeggel számoltunk. A [−π, π] × [−π, π] négyzetet s˝ur˝u ráccsal borítjuk le, minél nagyobb a felbontás, annál jobb a közelítésünk. Nézzünk meg néhány képet a felbontásról:

26

3.15. ábra. Felbontás 10 ponttal


27


Most nézzük meg, hogy mennyi a Riemann-integrálközelít˝o összeggel számított érték, ha a felbontás növelésével jobban és jobban közelítünk a "pontos" integrálhoz:

3.7. táblázat. Integrál értéke Riemann-integrálközelít˝o összeggel

Feltételezhetjük, hogy az f (x, y) függvényünk integrálja 5, 0645. Matlab program segítségével számoljuk ki az f (x, y) függvényünk közelítését a már ismert módon. El˝oször is generálunk n (n = 1500) pontot a [−π, π] × [−π, π] négyzeten, 28

viszont csak azokkal a pontokkal számolunk tovább, amelyek az am˝obaszer˝u tartományon belül esnek (piros szín˝uek),bármely 2 pont távolsága legalább ε legyen (ebben az esetben ε = 0.02):

3.18. ábra. n db pont generálása az am˝obaszer˝u tartományon

A tesztfüggvénnyel számított közelítés:

3.19. ábra. Tesztfüggvény

29

Itt látható az am˝oba feletti rész:

3.20. ábra. TPS-interpolációs függvény (1)

Egy másik szögb˝ol készített kép:

3.21. ábra. TPS-interpolációs függvény (2)

30

Ezen a képen látható a hiba, azaz a bázisfüggvényt˝ol való eltérés. A hiba mértéke az ábra léptékében igen kicsi, azaz a TPS-interpoláció pontossága nagyon jó!


A hiba értéke az am˝oba körvonalán a legnagyobb.

3.23. ábra. A tesztfüggvény és a TPS-interpolációs függvény eltérésének nagyítása

31

Itt is nézzük meg hogyan változik a hiba értéke az alappontok (n) és perempontok (m) növekedésével. Mivel a [−π, π] × [−π, π] négyzetet szórjuk meg pontokkal, ezért az eddigiekhez képest nagy pontszámmal kell indulnunk, mert csak azokkal számolunk tovább, amelyek az am˝obán belül esnek, s ha kevés pontot adunk meg, akkor lehetséges, hogy 1 pont sem esik belülre.


Az am˝oba peremén veszünk fel pontokat, s itt már nem ekvidisztáns a felbontás:

3.24. ábra. Az am˝oba peremén felvett pontok

32

Ebben az esetben az irányvektorokat a következ˝oképpen határozzuk meg: (y 0 (t), −x0 (t)) n(t) = p , ((x0 (t))2 + (y 0 (t))2 )2 ahol

x0 (t) = esin(t) cos2 (t) sin2 (2t) + esin(t) (− sin(t)) sin2 (2t) + esin(t) cos(t)4 sin(2t) cos(2t)+ +ecos(t) (− sin(t)) cos(t) cos2 (2t)+ecos(t) (− sin(t)) cos2 (2t)+ecos(t) cos(t)4 cos(2t)(− sin(2t)),

y 0 (t) = esin(t) cos(t) sin(t) sin2 (2t)+esin(t) cos(t) sin2 (2t)+esin(t) sin(t)4 sin(2t) cos(2t)+ +ecos(t) (− sin2 (t)) cos2 (2t) + ecos(t) cos(t) cos2 (2t) + ecos(t) sin(t)4 cos(2t)(− sin(2t))

A perempontokhoz az irányvektorokat a következ˝o ábra jól szemlélteti:

3.25. ábra. Az am˝oba perempontjainak irányvektorai

33

A Matlab programmal számolt eredményeket mutatja a következ˝o táblázat, melyben megnéztük hogyan közelít az integrál értéke, ha az alappontok számát növelem, a perempontok száma fix (m = 200):


Vizsgáljuk meg azt is, hogy mennyire változik az integrál értéke, ha a peremen fevett pontok számát növelem egy fix (n = 1000) véletlen pontszám mellett:


Az am˝obaszer˝u tartományon vett f (x, y) függvényünk integrálja jól közelít a Riemannintegrálközelít˝o összeggel számított értékhez.

34

Összefoglalás Olyan módszer bemutatása volt a cél, ami a kétváltozós általános tartományon vett integrált visszavezeti egy peremmenti integrál kiszámításra, ami numerikusan sokkal egyszer˝ubb, mint a kétváltozós függvény integrálása, mivel az már egy egyváltozós integrál. Ezzel a módszerrel megkerüljük a tartomány sokszögekre való felbontását s ehelyett szórt alappontú interpolációs technikát alkalmazunk. Három különböz˝o tartományon vizsgáltunk különböz˝o függvényeknek az integrálját s ezek közelítéseit vizsgáltuk. Az els˝o tartomány az egységnégyzet volt, s itt az f (x, y) = sin(πx) sin(πy) függvény integrálját közelítettük. A kapott eredmények jól közelítették a pontos integrál értékét, mind az egységnégyzetet megszórt véletlenpontok, mind a négyzet peremén fel vett pontok növekedésekor, s adott értékek mellett már a pontos integrált kaptuk. A második tartomány az egységsugarú kör volt, amin a félgömb térfogatát, azaz az p f (x, y) = 1 − x2 − y 2 integrálját számoltuk ki s hasonlítottuk össze a pontos térfogattal. Itt is jól közelített a számolt eredményünk a pontos térfogathoz. A harmadik egy am˝obaszer˝u tartomány volt. Itt az f (x, y) = cos( x2 ) cos( y2 ) függvény integrálját közelítettük az adott tartományon. Itt a véletlen pontok számát kell˝oen nagynak kellett vennünk, mivel csak azokkal a pontokkal foglalkoztunk a további számolás során, amelyek az am˝obaszer˝u tartományon belül estek. Az f (x, y) függvényünk pontos integrálját nem számíthattuk ki, mint korábban, így ezt egy Riemann-integrálközelít˝o összeggel kaptuk meg, s ezt felhasználva kaptuk meg az integrál számított értékének relatív hibáját.

35

Irodalomjegyzék [1] Stoyan Gisbert, Takó Galina, Numerikus módszerek III., Typotex, 1997. [2] Stoyan Gisbert, Numerikus matematika mérnököknek és programozóknak, Typotex, 2007. [3] Franke. R, Scattered Data Interpolation: Test of Some Methods, Mathematics of Computation, Vol. 38, No. 157, 1982. [4] Partridge P.W., Brebbia C.A., Computer Implementation of the BEM Dual Reciprocity Method for the Solution of general Field Equations, Communication in Applied Numerical Methods, Vol. 6, 83-92., 1990. [5] Tang W., Transforming Domain into Boundary Integrals in BEM, A Generalized Approach Lecture Notes in Engineering (ed. by C.A. Brebbia, S.A. Orszag), Vol. 35, SpringerVerlag, Berlin, Heidelberg, New York, London, Paris, Tokyo, 1988

36

Függelék Itt találhatók azok a programok, melyek segítségével közelítettem az integrál pontos értékét, mellyel szakdolgozatomban foglalkoztam.

negyzet.m f u n c t i o n n e g y z e t ( n ,m) % m e g h í v á s a : n e g y z e t ( n ,m) , a h o l % n − a v é l e t l e n p o n t o k száma % m − a n é g y z e t peremén f e l v e t t p o n t o k száma a = 0; b = 1; %−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−− % n db p o n t g e n e r á l á s a a [ 0 , 1 ] x [ 0 , 1 ] n é g y s z ö g b e n a % p o n t o k _ g e n e r a l a s a .m f á j l b a n t ö r t é n i k ; k é t p o n t % t á v o l s á g a l e g a l á b b eps l e s z [x , y] = pontok_generalasa (a , b , n ) ; plot (x , y , ’∗ ’) t i t l e ( ’ n db . p o n t g e n e r á l á s a ’ ) pause %−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−− % T e s z t f ü g g v é n y d e f i n i á l á s a , e z t a t e s z t f v .m f á j l b a í r j u k be f = tesztfv (x , y ); %−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−− % A TPS−i n t e r p o l á c i ó s m á t r i x e l o˝ á l l í t á s a A = zeros (n , n ) ; for i = 1: n for j = 1: n i f i == j A( i , j ) = 0 ; else r = norm ( [ x ( i ) y ( i ) ] − [ x ( j ) y ( j ) ] ) ; A( i , j ) = r ∗ r ∗ l o g ( r ) ; end 37

end end %−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−− % TPS−i n t e r p o l á c i ó e g y ü t t h a t ó i n a k s z á m í t á s a a l f a = A\ f ; %−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−− % Az i n t e r p o l á c i ó s f ü g g v é n y m e g j e l e n í t é s e f e l b o n t a s = 40; xf = l i n s p a c e ( a , b , f e l b o n t a s ) ; yf = l i n s p a c e ( a , b , f e l b o n t a s ) ; [ x i , y i ] = m e s h g r i d ( xf , y f ) ; zi=zeros ( felbontas ) ; for i = 1: felbontas for j = 1: felbontas for k = 1: n r = norm ( [ x i ( i , j ) y i ( i , j ) ] − [ x ( k ) y ( k ) ] ) ; i f r ~=0 z i ( i , j )= z i ( i , j )+ a l f a ( k )∗ r ∗ r ∗ log ( r ) ; end end end end s u r f ( xi , yi , z i ) ; t i t l e ( ’ TPS−i n t e r p o l á c i ó s f ü g g v é n y ’ ) pause %−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−− % Az i n t e r p o l á c i ó h i b á j á n a k m e g j e l e n í t é s e f i = t e s z t f v ( xi , y i ) ; s u r f ( xi , yi , f i − z i ) ; a x i s ( [ 0 1 0 1 −1 1 ] ) t i t l e ( ’A t e s z t f v . é s a TPS−i n t e r p o l á c i ó s f ü g g v é n y e l t é r é s e ’ ) pause %−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−

38

% A peremen f e l v e t t p o n t o k d e f i n i á l á s a x1 = z e r o s ( 1 ,m− 1 ) ; xx1 = l i n s p a c e ( a , b ,m ) ; y1 = z e r o s ( 1 ,m− 1 ) ; f o r i = 1 :m−1 x1 ( i ) = ( xx1 ( i ) + xx1 ( i + 1 ) ) / 2 ; end %−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−− y2 = z e r o s ( 1 ,m− 1 ) ; x2 = z e r o s ( 1 ,m−1)+1; yy2 = l i n s p a c e ( a , b ,m ) ; f o r i = 1 :m−1 y2 ( i ) = ( yy2 ( i ) + yy2 ( i + 1 ) ) / 2 ; end %−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−− x3 = z e r o s ( 1 ,m− 1 ) ; xx3 = l i n s p a c e ( b , a ,m ) ; y3 = z e r o s ( 1 ,m−1)+1; f o r i = 1 :m−1 x3 ( i ) = ( xx3 ( i ) + xx3 ( i + 1 ) ) / 2 ; end %−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−− y4 = z e r o s ( 1 ,m− 1 ) ; x4 = z e r o s ( 1 ,m− 1 ) ; yy4 = l i n s p a c e ( b , a ,m ) ; f o r i = 1 :m−1 y4 ( i ) = ( yy4 ( i ) + yy4 ( i + 1 ) ) / 2 ; end x k o r = [ x1 x2 x3 x4 ] ; y k o r = [ y1 y2 y3 y4 ] ; r = [0 1 1 0 0 ] ; t = [0 0 1 1 0 ] ; p l o t ( r , t , xkor , ykor , ’ ∗ r ’ )

39

%−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−− % Irányvektorok n1_x = z e r o s ( 1 ,m− 1 ) ; n1_y = z e r o s ( 1 ,m−1) −1; n2_x = z e r o s ( 1 ,m−1)+1; n2_y = z e r o s ( 1 ,m− 1 ) ; n3_x = z e r o s ( 1 ,m− 1 ) ; n3_y = z e r o s ( 1 ,m−1)+1; n4_x = z e r o s ( 1 ,m−1) −1; n4_y = z e r o s ( 1 ,m− 1 ) ; nx = [ n1_x n2_x n3_x n4_x ] ; ny = [ n1_y n2_y n3_y n4_y ] ; %−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−− M = 4 ∗ (m− 1 ) ; h = x1 ( 2 ) − x1 ( 1 ) ; dx = z e r o s ( s i z e ( x k o r ) ) ; dy = z e r o s ( s i z e ( x k o r ) ) ; dpdn = z e r o s ( s i z e ( x k o r ) ) ; f o r i = 1 :M for j = 1: n r = norm ( [ x k o r ( i ) y k o r ( i ) ] − [ x ( j ) y ( j ) ] ) ; dx ( i ) = dx ( i ) + ( a l f a ( j ) ∗ ( 1 / 1 6 ) ∗ ∗ ( ( x k o r ( i )−x ( j ) ) ∗ ( r ^ 2 ) ∗ ( 2 ∗ l o g ( r ^ 2 ) − 1 ) ) ) ; dy ( i ) = dy ( i ) + ( a l f a ( j ) ∗ ( 1 / 1 6 ) ∗ ∗ ( ( y k o r ( i )−y ( j ) ) ∗ ( r ^ 2 ) ∗ ( 2 ∗ l o g ( r ^ 2 ) − 1 ) ) ) ; end dpdn ( i ) = ( dx ( i ) ∗ nx ( i ) ) + ( dy ( i ) ∗ ny ( i ) ) ; end i n t e g r a l = h ∗ sum ( dpdn ) ; display ( integral )

40

kor.m f u n c t i o n k o r ( n ,m) % m e g h í v á s a : k o r ( n ,m) , a h o l % n − a v é l e t l e n p o n t o k száma % m − a n é g y z e t peremén f e l v e t t p o n t o k száma a = −1; b = 1 ; %−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−− % Kör k i r a j z o l á s a t = l i n s p a c e ( 0 , 2∗ p i , 1 0 0 ) ; p l o t ( cos ( t ) , s i n ( t ) ) axis equal pause h o l d on %−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−− % n db p o n t g e n e r á l á s a a z e g y s é g k ö r b e n a % p o n t o k _ g e n e r a l a s a _ k o r .m f á j l b a n t ö r t é n i k ; k é t p o n t % t á v o l s á g a l e g a l á b b eps l e s z [x , y] = pontok_generalasa_kor (a , b , n ) ; p l o t ( x , y , ’ g ∗ ’) t i t l e ( ’ n db . p o n t g e n e r á l á s a ’ ) pause hold off %−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−− % T e s z t f ü g g v é n y d e f i n i á l á s a , e z t a t e s z t f v _ k o r .m % f á j l b a í r j u k be f = tesztfv_kor (x , y ); %−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−− % A TPS−i n t e r p o l á c i ó s m á t r i x e l o˝ á l l í t á s a A = zeros (n , n ) ; for i = 1: n for j = 1: n i f i == j A( i , j ) = 0 ; else 41

r = norm ( [ x ( i ) y ( i ) ] − [ x ( j ) y ( j ) ] ) ; A( i , j ) = r ∗ r ∗ l o g ( r ) ; end end end %−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−− % TPS−i n t e r p o l á c i ó e g y ü t t h a t ó i n a k s z á m í t á s a a l f a = A\ f ; %−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−− % Az i n t e r p o l á c i ó s f ü g g v é n y m e g j e l e n í t é s e f e l b o n t a s = 40; rf = linspace (0 ,0.9999 , felbontas ) ; f i f = l i n s p a c e (0 ,2∗ pi , f e l b o n t a s ) ; [ r , f i ] = meshgrid ( rf , f i f ) ; [ xi , y i ] = p o l 2 c a r t ( f i , r ) ; zi=zeros ( felbontas ) ; for i = 1: felbontas for j = 1: felbontas for k = 1: n r = norm ( [ x i ( i , j ) y i ( i , j ) ] − [ x ( k ) y ( k ) ] ) ; i f r ~=0 z i ( i , j )= z i ( i , j )+ a l f a ( k )∗ r ∗ r ∗ log ( r ) ; end end i f x i ( i , j )^2+ y i ( i , j )^2 > 1 ; z i ( i , j ) = NaN ; end end end s u r f ( xi , yi , z i ) ; axis equal t i t l e ( ’ TPS−i n t e r p o l á c i ó s f ü g g v é n y ’ ) pause %−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−

42

% Az i n t e r p o l á c i ó h i b á j á n a k m e g j e l e n í t é s e f i = t e s z t f v _ k o r ( xi , y i ) ; s u r f ( xi , yi , f i − z i ) ; a x i s ([ −1 1 −1 1 −1 1 ] ) t i t l e ( ’A t e s z t f v . é s a TPS−i n t e r p o l á c i ó s f ü g g v é n y e l t é r é s e ’ ) pause %−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−− %K ö r v o n a l o n l é v o˝ p o n t o k x1 = z e r o s ( 1 ,m ) ; y1 = z e r o s ( 1 ,m ) ; f o r i = 1 :m x1 ( i ) = c o s ( i ∗ ( ( 2 ∗ p i ) / m ) ) ; end f o r j = 1 :m y1 ( j ) = s i n ( j ∗ ( ( 2 ∗ p i ) / m ) ) ; end x k o r = x1 ; y k o r = y1 ; p l o t ( c o s ( t ) , s i n ( t ) , xkor , ykor , ’ r ∗ ’ ) a x i s ( [ − 1 . 1 1 . 1 −1.1 1 . 1 ] ) %−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−− %I r á n y v e k t o r o k nx = x1 ; ny = y1 ; %−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−− h = norm ( [ x1 ( 2 ) y1 ( 2 ) ] − [ x1 ( 1 ) y1 ( 1 ) ] ) ; dx = z e r o s ( s i z e ( x k o r ) ) ; dy = z e r o s ( s i z e ( x k o r ) ) ; dpdn = z e r o s ( s i z e ( x k o r ) ) ; f o r i = 1 :m for j = 1: n r = norm ( [ x k o r ( i ) y k o r ( i ) ] − [ x ( j ) y ( j ) ] ) ;

43

dx ( i ) = dx ( i ) + a l f a ( j ) ∗ ( 1 / 1 6 ) ∗ ∗ ( ( x k o r ( i )−x ( j ) ) ∗ ( r ^ 2 ) ∗ ( 2 ∗ l o g ( r ^ 2 ) − 1 ) ) ; dy ( i ) = dy ( i ) + a l f a ( j ) ∗ ( 1 / 1 6 ) ∗ ∗ ( ( y k o r ( i )−y ( j ) ) ∗ ( r ^ 2 ) ∗ ( 2 ∗ l o g ( r ^ 2 ) − 1 ) ) ; end dpdn ( i ) = ( dx ( i ) ∗ nx ( i ) ) + ( dy ( i ) ∗ ny ( i ) ) ; end i n t e g r a l = h ∗ sum ( dpdn ) ; display ( integral )

44

amoba.m f u n c t i o n amoba ( n ,m) % m e g h í v á s a : amoba ( n ,m) , a h o l % n − a v é l e t l e n p o n t o k száma % m − a n é g y z e t peremén f e l v e t t p o n t o k száma a = −p i ; b = p i ; %−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−− % Am˝oba k i r a j z o l á s a t = l i n s p a c e ( 0 , 2∗ p i , 1 0 0 ) ; r = exp ( s i n ( t ) ) . ∗ ( ( s i n ( 2 ∗ t ) ) . ^ 2 ) + exp ( c o s ( t ) ) . ∗ ( ( c o s ( 2 ∗ t ) ) . ^ 2 ) ; xv = r . ∗ c o s ( t ) ; yv = r . ∗ s i n ( t ) ; p l o t ( xv , yv ) a x i s ([ −4 4 −4 4 ] ) pause %−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−− % n db p o n t g e n e r á l á s a a [− p i , p i ] x[− p i , p i ] n é g y s z ö g b e n , % amely a p o n t o k _ g e n e r a l a s a _ a m o b a .m f á j l b a n t ö r t é n i k ; % két pont t á v o l s á g a l e g a l á b b eps l e s z [x , y ] = pontok_generalasa_amoba ( a , b , n ) ; i n = i n p o l y g o n ( x , y , xv , yv ) ; p l o t ( xv , yv , x ( i n ) , y ( i n ) , ’ r . ’ , x ( ~ i n ) , y ( ~ i n ) , ’ b . ’ ) pause xx = x ( i n ) ; yy = y ( i n ) ; kx = x ( ~ i n ) ; ky = y ( ~ i n ) ; [ m e r e t ,mm] = s i z e ( xx ) ; %h o l d o f f %−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−− % T e s z t f ü g g v é n y d e f i n i á l á s a , e z t a t e s z t f v _ a m o b a .m % f á j l b a í r j u k be f = t e s z t f v _ a m o b a ( xx , yy ) ; %−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−− % A TPS−i n t e r p o l á c i ó s m á t r i x e l o˝ á l l í t á s a 45

A = z e r o s ( meret , meret ) ; f o r i = 1: meret f o r j = 1: meret i f i == j A( i , j ) = 0 ; else r = norm ( [ xx ( i ) yy ( i ) ] − [ xx ( j ) yy ( j ) ] ) ; A( i , j ) = r ∗ r ∗ l o g ( r ) ; end end end %−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−− % TPS−i n t e r p o l á c i ó e g y ü t t h a t ó i n a k s z á m í t á s a a l f a = A\ f ; %−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−− % Az i n t e r p o l á c i ó s f ü g g v é n y m e g j e l e n í t é s e f e l b o n t a s = 40; xf = l i n s p a c e ( a , b , f e l b o n t a s ) ; yf = l i n s p a c e ( a , b , f e l b o n t a s ) ; [ x i , y i ] = m e s h g r i d ( xf , y f ) ; zi = zeros ( felbontas , felbontas ) ; for i = 1: felbontas for j = 1: felbontas f o r k = 1: meret r = norm ( [ x i ( i , j ) y i ( i , j ) ] − [ xx ( k ) yy ( k ) ] ) ; i f r > 0.0001 zi ( i , j ) = zi ( i , j ) + a l f a ( k )∗ r ∗ r ∗ log ( r ) ; end end end end i n = i n p o l y g o n ( x i , y i , xv , yv ) ; s u r f ( xi , yi , i n . ∗ z i ) ; t i t l e ( ’ TPS−i n t e r p o l á c i ó s f ü g g v é n y ’ )

46

pause %−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−− % Az i n t e r p o l á c i ó h i b á j á n a k m e g j e l e n í t é s e f i = t e s z t f v _ a m o b a ( xi , y i ) ; s u r f ( xi , yi , i n . ∗ ( f i − z i ) ) ; a x i s ([ −4 4 −4 4 0 1 ] ) t i t l e ( ’A t e s z t f ü g g v é n y é s a TPS f ü g g v é n y e l t é r é s e ’ ) pause s u r f ( xi , yi , i n . ∗ ( f i − z i ) ) ; t i t l e ( ’A t e s z t f ü g g v é n y é s a TPS f ü g g v é n y e l t é r é s e ’ ) pause %−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−− % % A peremen f e l v e t t p o n t o k d e f i n i á l á s a t = l i n s p a c e ( 0 , 2∗ p i , m ) ; r = exp ( s i n ( t ) ) . ∗ ( ( s i n ( 2 ∗ t ) ) . ^ 2 ) + exp ( c o s ( t ) ) . ∗ ( ( c o s ( 2 ∗ t ) ) . ^ 2 ) ; x1 = r . ∗ c o s ( t ) ; y1 = r . ∗ s i n ( t ) ; p l o t ( x1 , y1 , ’ ∗ g ’ ) pause %−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−− %I r á n y v e k t o r o k ex = exp ( s i n ( t ) ) . ∗ ( c o s ( t ) . ^ 2 ) . ∗ ( ( s i n ( 2 ∗ t ) ) . ^ 2 ) + + exp ( s i n ( t ) ) . ∗ ( − s i n ( t ) ) . ∗ ( ( s i n ( 2 ∗ t ) ) . ^ 2 ) + + exp ( s i n ( t ) ) . ∗ ( c o s ( t ) ) . ∗ ( 4 ∗ ( s i n ( 2 ∗ t ) . ∗ ( c o s ( 2 ∗ t ) ) ) ) + + exp ( c o s ( t ) ) . ∗ ( − s i n ( t ) ) . ∗ ( c o s ( t ) ) . ∗ ( ( c o s ( 2 ∗ t ) ) . ^ 2 ) + + exp ( c o s ( t ) ) . ∗ ( − s i n ( t ) ) . ∗ ( ( c o s ( 2 ∗ t ) ) . ^ 2 ) + + exp ( c o s ( t ) ) . ∗ ( c o s ( t ) ) . ∗ ( 4 ∗ ( c o s ( 2 ∗ t ) . ∗ ( − s i n ( 2 ∗ t ) ) ) ) ; ey = exp ( s i n ( t ) ) . ∗ ( c o s ( t ) ) . ∗ ( s i n ( t ) ) . ∗ ( ( s i n ( 2 ∗ t ) ) . ^ 2 ) + + exp ( s i n ( t ) ) . ∗ ( c o s ( t ) ) . ∗ ( ( s i n ( 2 ∗ t ) ) . ^ 2 ) + + exp ( s i n ( t ) ) . ∗ ( s i n ( t ) ) . ∗ ( 4 ∗ ( s i n ( 2 ∗ t ) . ∗ ( c o s ( 2 ∗ t ) ) ) ) + + exp ( c o s ( t ) ) . ∗ ( − s i n ( t ) . ^ 2 ) . ∗ ( ( c o s ( 2 ∗ t ) ) . ^ 2 ) + + exp ( c o s ( t ) ) . ∗ ( c o s ( t ) ) . ∗ ( ( c o s ( 2 ∗ t ) ) . ^ 2 ) + + exp ( c o s ( t ) ) . ∗ ( s i n ( t ) ) . ∗ ( − 4 ∗ ( c o s ( 2 ∗ t ) . ∗ s i n ( 2 ∗ t ) ) ) ;

47

nx = z e r o s ( 1 ,m ) ; ny = z e r o s ( 1 ,m ) ; n e v e z o = s q r t ( ex . ^ 2 + ey . ^ 2 ) ; f o r i = 1 :m nx ( i ) = ey ( i ) / n e v e z o ( i ) ; ny ( i ) = −ex ( i ) / n e v e z o ( i ) ; end h o l d on q u i v e r ( x1 , y1 , nx , ny ) pause %−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−− dx = z e r o s ( 1 ,m ) ; dy = z e r o s ( 1 ,m ) ; dpdn = z e r o s ( 1 ,m ) ; f o r i = 1 :m−1 h = norm ( [ x1 ( i + 1) y1 ( i + 1 ) ] − [ x1 ( i ) y1 ( i ) ] ) ; f o r j = 1: meret r = norm ( [ x1 ( i ) y1 ( i ) ] − [ xx ( j ) yy ( j ) ] ) ; dx ( i ) = dx ( i ) + ( a l f a ( j ) ∗ ( 1 / 1 6 ) ∗ ∗ ( ( x1 ( i )−xx ( j ) ) ∗ ( r ^ 2 ) ∗ ( 2 ∗ l o g ( r ^ 2 ) − 1 ) ) ) ; dy ( i ) = dy ( i ) + ( a l f a ( j ) ∗ ( 1 / 1 6 ) ∗ ∗ ( ( y1 ( i )−yy ( j ) ) ∗ ( r ^ 2 ) ∗ ( 2 ∗ l o g ( r ^ 2 ) − 1 ) ) ) ; end dpdn ( i ) = ( ( dx ( i ) ∗ nx ( i ) ) + ( dy ( i ) ∗ ny ( i ) ) ) ∗ h ; end i n t e g r a l = sum ( dpdn ) ; display ( integral )

48

riemann.m f u n c t i o n riemann ( f e l b o n t ) x = l i n s p a c e (− p i , p i , f e l b o n t ) ; y = l i n s p a c e (− p i , p i , f e l b o n t ) ; [X, Y] = m e s h g r i d ( x , y ) ; Z = c o s (X . / 2 ) . ∗ c o s (Y . / 2 ) ; s u r f (X, Y, Z ) pause t = l i n s p a c e ( 0 , 2∗ p i , 1 0 0 ) ; r = exp ( s i n ( t ) ) . ∗ ( ( s i n ( 2 ∗ t ) ) . ^ 2 ) + exp ( c o s ( t ) ) . ∗ ( ( c o s ( 2 ∗ t ) ) . ^ 2 ) ; xv = r . ∗ c o s ( t ) ; yv = r . ∗ s i n ( t ) ; i n = i n p o l y g o n (X, Y, xv , yv ) ; p l o t ( xv , yv , X( i n ) ,Y( i n ) , ’ r . ’ ,X( ~ i n ) ,Y( ~ i n ) , ’ g . ’ ) pause xk = X( i n ) ; yk = Y( i n ) ; h = 2∗ p i / ( f e l b o n t − 1 ) ; [ s o r o k , o s z l o p o k ] = s i z e ( xk ) ; i n t e g r a l = 0; f o r j =1: sorok i n t e g r a l = i n t e g r a l + t e s z t f v _ a m o b a ( xk ( j ) , yk ( j ) ) ∗ h ^ 2 ; end display ( integral )

49

Szakdolgozat. Írta: Újvári Zsuzsanna. Matematika BSc Alkalmazott matematika szakirány

Recommend Documents