EÖTVÖS LORÁND TUDOMÁNYEGYETEM TERMÉSZETTUDOMÁNYI KAR
Bisztray Tamás
Geodetikus g¨ orb´ ek keres´ ese fel¨ uleteken BSc Szakdolgozat
Témavezető: Dr. Szeghy Dávid Geometria Tanszék
Budapest, 2015
Ko ¨sz o ¨netnyilv´ an´ ıt´ as Szeretném megköszönni Szeghy Dávidnak a témavezetés során nyújtott tanácsait és útmutatását, valamint, hogy rendszeres konzultációkkal segítette munkámat. Szeretnék köszönetet mondani továbbá Lócsi Leventének, amiért számíthattam a segítségére MATLAB-os programozási kérdésekben. Végül köszönettel tartozom szüleimnek a támogatásukért.
Tartalomjegyzék Bevezető . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 Geodetikus görbék, alapfogalmak . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 Középpont módszer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 Alkalmazás a síkban . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . 4 Térbeli kiterjesztés. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 Differenciálgeometriai megközelítés . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 Fixpont iteráció. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 Newton-Raphson módszer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 Gráfelméleti megközelítés . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 Szélességi keresés . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 Dijkstra algoritmus. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 Mohó legjobb út keresés . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 A⇤ algoritmus. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 Összefoglalás . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 Függelék
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
Bevezető A különböző felületeken való útkeresés a gyakorlati matematikában az egyik leggyakrabban felmerülő feladat. Ezen belül kimelt fontosságú a legrövidebb vagy minimális költségű út megtalálása. Sima felületeken ezt a kérdést a visszavezethetjük a geodetikus görbék keresésére, melyek lokálisan minimális hosszúak és ezen görbék családjából kerül ki a globálisan legrövidebb út is. Kontinensek közötti repülő útvonlak tervezésekor például a gép a célállomás és a felszállási ponton átmenő főkör rövidebbik szakaszán halad, s mint később látni fogjuk ez éppen egy geodetikus görbe. GPS útvonaltervezésben a legrövidebb út meghatározása két pont között az egyik legalapvetőbb probléma. Egy számítógépes játékban az NPC karakterek mozgatása szintén ebbe a kérdéskörbe tartozik. Merre érdemes a karaktert elindítani, esetleg akadályokat kell kikerülni vagy egy ívhossz szerint hosszabb utat érdemes preferálni egy rövidebb, de lassabb útvonallal szemben. Például megéri megkerülni egy tavat vagy gyorsabb azt átúszva eljutni a célhoz. Az alkalmazott matematikának több ága is foglalkozik e problémák megoldásával. A dolgozat célja, hogy több megoldást is mutasson geodetikus görbék megtalálására különböző területekről.
1
Geodetikus görbék, alapfogalmak Definíció: Az R3 tér egy M részhalmazát sima elemi felületnek nevezzük, ha van olyan R2 -beli U tartomány (nyílt összefüggő halmaz) és olyan C 1 osztályú r : U ! R3 leképezés, amelyre: • r(U ) = M •
Bármely u 2 U esetén a @1 r(u), @2 r(u) vektorok lineárisan függetlenek
•
r homeomorfizmus U és M között
Ekkor az r függvényt az M = r(U ) sima elemi felület egy paraméteres előállításának nevezzük. Definíció: Egy R3 -beli M halmazt sima felületnek mondunk, ha az M tetszőleges p T pontjának van olyan K nyílt környezete R3 -ben , hogy a K M metszet egy sima elemi felület. Definíció: Az R3 térben legyen adva egy M sima felület az r : U ! R3 paraméterezéssel. A y : I ! U síkgörbe áltan definiált felületi görbén a = r y : I ! M , t[ r(y(t)) leképezést értjük. A = r y felületi görbét az y görbe r szerinti képének nevezzük. Definíció: Az R3 térbeli folytonosan paraméterezett görbén egy : I ! R3 folytonos leképezését értjük. A (I) = { (t) | t✏I} ponthalmazt a görbe pályájának nevezzük. Számunkra azok a görbék lesznek érdekesek melyek egy adott felületen fekszenek. Tétel: Egy görbe ívhossza egyenlő a ráírt töröttvonalak hosszának határértékével. Legyen I intervallum egy felosztása az I = (t0 , ..., tm ) és : I ! R3 leképezés. A S görbét ezen felosztás alapján az m k=1 (tk 1 ) (tk ) töröttvonallal közelíthetjük, melynek hossza: m X k (tk ) (tk 1 )k k=1
Ekkor a görbe ívhossza a töröttvonalak hosszának határértéke: l=
lim
maxk (tk tk
1 )!0
(
m X k=1
k (tk )
(tk 1 )k)
Célunk a felületen lévő két pontot összekötő felületi görbék közül a legrövidebb megtalálása. Bizonyítható, hogy ez a legrövidebb görbe az úgynevezett geodetikus görbék közül kerül ki. A következőkben értelmezzük a geodetikus görbe fogalmát. 2
Definíció: Legyen M egy sima felület r : U ! M paraméterezéssel. Igazolható, hogy az u = (u1 , u2 ) paraméterezésű pontban az érintősík normálvektora: @1 r(u) ⇥ @2 r(u) Megjegyezhetjük, hogy ez 6= 0, mivel @1 r(u), @2 r(u) lineárisan függetlenek. Így a fenti vektor a felületre merőleges (normális) vektornak tekinthető. Definíció: Az M felület r paraméterezéséhez tartozó normális egységvektormezőn azt az N : U ! R3 leképezést értjük melyet az N (u1 , u2 ) =
@1 r(u1 , u2 ) ⇥ @2 r(u1 , u2 ) k@1 r(u1 , u2 ) ⇥ @2 r(u1 , u2 )k
egyenlet határoz meg tetszőleges (u1 , u2 ) 2 U esetén. A geodetikus görbékre általánosan igaz, hogy lokálisan a görbe két pontja között egy legrövidebb utat adnak. Ebből következik, hogy két pont között haladó legrövidebb út maga is geodetikus. Azonban egy tetszőleges geodetikus görbe globálisan nem feltétlenül a legrövidebb út két pont között. A geodetikusságra vonatkozó feltételeket a későbbiekben részletesebben tárgyaljuk. Igazolható a következő tétel ([4],8.11 állítás): Tétel: Legyenek P, Q 2 U rögzített pontok. Az r(P ) és r(Q) pontokat összekötő 00 = r y görbe akkor és csak akkor geodetikus ha,: 8t 2 I : (t) és N (y(t)) vektorok lineárisan öszefüggőek.
3
Középpont módszer A következőkben az [5] cikkben tárgyalt első módszert vizsgáljuk meg részletesebben, melynél adott felületen két előre rögzített pont között keresünk geodetikus utat. Az algoritmus lényege, hogy numerikus módszerek alkalmazásával tetszőleges kezdeti útra felhasználható. Legyen p kezdő és q végpont, valamint adott egy kezdeti út mely köztük vezet. Ennek javításval szereténk olyan utat kapni mely jól közelít egy, a két pont között húzódó geodetikus görbét. A problémát először síkban vizsgáljuk majd az itt alkalmazott módszert kiterjesztjük térbeli felületekre.
Alkalmazás síkban Legyen p,q2R2 . Tudjuk, hogy síkban két pont között vezető legrövidebb út az egyenes. A gyakorati alkalmazhatóság szempontjából fel kell tegyük, hogy minimum egy köztes pont adott. Legyen 0 = (p01 , ..., p0n ) n > 2 egy véges kezdeti pontsorozat, p0i 6= p0j és S S S p = p01 , q = p0n a végpontok. Az utat a p1 p2 p2 p3 ....... pn 1 pn , adja ahol p0i = (x0i , yi0 ) koordinátákkal írható fel. Legyen 0 m k. Kezdetben m = 0 és k értéke határozza meg, hány iterációt végzünk. Nézzük a következő algoritmust. •
m
útra 8 pi 1 pi szakszra kiszámoljuk a felezőpontot azaz
m pm i 1 +pi 2
i = 2, ..., n + m.
•
A kapott pontsorozathoz kezdőpontnak hozzáveszük p01 míg végpontnak a p0n pontot. m 0 Ezt az utat m+1 -nek nevezzük. m+1 = (p01 , pm 2 ...pn+m , pn )
•
m értékét egyel növeljük. Ha m < k az első pontra ugrunk egyébként kilépünk a ciklusból
A fenti jelöléseket használva legyen l( m ) = töröttvonal hossza. Állítás: 8m = 1, 2, 3, ... esetén `( m 1 ) `(
4
Pn+m i=2
m
).
pm i
pm i 1
az m-edik iterációs
1 pm i 1
pm i 1
pm i
1 pm i 2
pm i
1
Bizonyítás: Írjuk fel egy (m 1)-edik iterációs töröttvonal hosszát. Ez éppen l( m 1 ) = 1 1 Pn+m 1 m 1 Pn+m 1 P2 kpm pm 1 i i 1 k pi pm . i=2 j=1 i 1 . Ez felírható a következő formában: 2 i=2 Ez megfelel annak, hogy a szakaszok hossza helyett kétszer a félszakaszok hosszát összegez1 m 1 1 m m 1 m zük. A fenti ábrán ez például pm = pm felírást jelenti. i 2 pi 1 i 2 pi 1 + pi 1 pi 1 m 1 1 P P pm 1 kp i i 1 k Mivel j-től független a kifejezés, a szummákat felcserélhetjük: 2j=1 n+m . i=2 2 A háromszögegyenlőtlenségre hivatkozva szeretnénk belátni, hogy m 1 pm i 1 pi 1
+
1 m pm i 1 pi
=
1 pm i 1
1 pm i 2
pm i
+
2
1
1 pm i 1
pm i
2
pm i 1
miatt a teljes görbén elvégezve az iterációt annak hossza csökken. Nézzük az alábbi átalakításokat l(
m 1
)=
2 n+m X X 1 pm i j=1
i0
1
1 pm i 1
2
i=2
=
n+m X1
pm i
1
1 pm i 1
2
i=2
+
n+m X1
pm i
1
1 pm i 1
2
i=2
0
1 = i, i = i + 1 = 3, ...., n + m indexeltolást alkalmazzuk az első tagra n+m X 0
i =3
pm i0
1 1
pm i0 2
1 2
+
n+m X1 i=2
pm i
1
1 pm i 1
2
=
0
i -t visszacseréljük i-re és leválasztjuk az i = 2, m tagokat. " n+m 1 1 m 1 m 1 m 1 1 1 X pm p p p pm pm pm n+m 1 n+m 2 2 1 i 1 i 2 + + + i 2 2 2 i=3
5
1
1 pm i 1
2
#
(4egyenl˝ otlens´ eg)
z
kp m n+m
1 pm n+m
1
1
}|
pm n+m k
1 pm n+m
2
pm kp m 2 1 k { z }| pm 1 pm 2 1 + 2 2 n+m X
pm i
1
{
+
n+m X1
pm i
pm i 1 =
i=3
pm i 1 = l(
m
)
i=2
⌅
Állítás: Legyen egy véges pontsorozat által definiált út p kezdő és q végpont között. A középpont módszer iteratív alkalmazásával egy olyan görbesorozatot kapunk, mely egy egyenes szakaszhoz tart. Bizonyítás: Látható, hogy a minden iterációt követően egyel növekszik az utat alkotó pontok száma. A felezőpontok kiszámításakor egy n hosszú pontsorozat esetén n 1 köztes pontot számolunk ki, majd a kezdő és végpont hozzávételét követően n + 1 pont adódik. Így látható, hogy m iteráció elvégzésével egy k + m hosszú út keletkezik. Legyen m m m m m = (pm 1 , ..., pn+m ) az m-edik iteráció után kapott pontsorozat ahol pi = (xi , yi ). Ha például ki akarjuk számolni az első iterációből keletkező pontsorozat egy i-edik kop0 +p0 ordinátáját az triviálisan p1i = i 12 i . A kezdő és végpontot az iteráció során rögzítjük. m = 2-re ez p2i
=
p1i
1
+ p1i = 2
p0i
0 2 +pi 1
2
+ 2
p0i
0 1 +pi
2
=
p0i
2
+ 2p0i 22
1
+ p0i
Ez alapján a pm i koordináta mindig felírható a kezdeti koordináták egy lineáris kombinációjaként osztva 2m -el.
p01 p11 p21 p02 \ p12 p22 p03 \ p13 \ p23 p14 p24 p25
6
A fenti ábrán n = 3, m = 2, i = 3 esetben a pont koordinátáit visszavezetjük a kezdeti pontsorozat koordinátáira. Látható, i > n esetén a kezdeti koordináták között nincs i-edik tag amely szükséges lenne a pont visszafejtéséhez. Hasonló helyzet adódik n kezdőpont esetén ha i 5 m idexű koordiátát számolunk ezzel a módszerrel. Mivel igazak 0 0 p0 +p0 n a p = p01 = 1 2 1 és q = p0n = pn +p összefüggések és a kezdő valamint a végpont végig 2 fix, ezért minden oszlopot egészítsünk ki a következőképp: p .. .
p .. .
p .. .
···
p p p ··· p02 p12 p22 · · · q p13 p23 q p24 .. . q .. .. .. . . . . . . q
q
q
···
Ez így már egyszerűen visszavezethető a Pascal-háromszög egy sorára mely az együtthatókat fogja szolgáltatni. A pm i koordináta felel meg a háromszög csúcsának, erre egyet írunk. Ennek m 1-edik generációs ősei az előző formula alapján 1 1 súllyal kerülnek be míg az m 2-edik iterációból 1 2 1, ami a Pascal-háromszög harmadik sora. Így egy m-edik iterációból származó pont felírásához - a kezdeti koordináták lineáris kombinációjaként - a Pascal-háromszög m-edik sora adja az együtthatókat. A fentiekből felírható a m · p0i m + m1 · p0i m+1 + ... + mm 1 p0i 1 + m p0 0 m i pm = i 2m összefüggés 8m 2 {0, 1, 2, 3, ...} esetén. p0i := p (i = 0, 1, 2, ...), p0i := q (i = n + 1, n + 2, ...) Ezt beláthatjuk indukcióval is. Lemma: 8m 2 {0, 1, 2, 3, ...}-re pm i
m ✓ ◆ 1 X m m = m· p 2 j i m+j j=0
(i 2 Z)
p0i := p (i = 0, 1, 2, ...), p0i := q (i = n + 1, n + 2, ...). 7
Lemma Bizonyítás: Teljes indukció m szerint: P0 0 0 j=0 j · pi 0 m = 0 : pi = 20 (m
0+j
= p0i
1) ! m :
1 m pm i 1 + pi pm = i 2
1
1 = 2
1 2m
1
m X1 ✓
·
m
1 j
j=0
◆
· p0(i
1) (m 1)+j
+
1 2m
1
·
m X1 ✓
m
j=0
1 j
◆
· p0i
(m 1)+j
index eltolás: k = j + 1, 1 = m 2
m X1 ✓
m
j=0
1 j
◆
·
p0(i 1) (m 1)+j
+
1 = m 2
✓
m
1 0
1 = m 2
◆
· p0i
m X1 ✓
m j
j=0
m+0
+
✓ ◆ m · p0i 0
1
◆
m X1 ✓
· p0i
m
j=1
m+0
k
k=1
k-t visszacseréljük j-re: 1 = m 2
m ✓ X m
+
k=1
1 j
◆
m X1 ✓ j=1
+
m+j
m ✓ X m
m j
✓
m + j ◆
· p0i
j
1 1
◆
m+j
m ✓ ◆ 1 X m m = m· p 2 j i m+j j=0
◆ 1 · p0i 1
◆ 1 · p0i 1
· p0i
m+j
m+k
m+j
+
✓ ◆ m + · p0i m
!
!
=
✓
m m
m+m
=
◆ 1 · p0i 1 !
m+m
!
=
⌅
Legyen 0 = (p01 , ..., p0n ) véges pontsorozat a síkban. A p01 pontot toljuk el az origóba. Az összes többi pontot egy ugyanilyen irányú és hosszúságú vektorral szintén toljunk el. Ezután végezzünk el egy origó körüli forgatást minden ponton, úgy, hogy p0n pontosan az x-tengelyre essen. Így p01 = (0, 0) és pn1 = (t, 0) koordinátájú pontok lesznek, valamint könnyen látható, hogy az őket összekötő legrövidebb út pontosan az x-tengelyre esik. Az 8
=
!
=
iterációk és a fenti egybevágósági transzformáció felcserélhető, azaz az eredeti sorozat m-edik iterációjának transzformáltja megegyezik az eredeti transzformáltjának m-edik iterációjával. Továbbá a legrövidebb utak is egymásba mennek. Ezért elég a transzformáció utáni speciális helyzetű pontsorozatra kimondani az állítást. Legyen az m-edik iteráció hibája maxi2Z |yim |. Az yim számokra ugyanaz a képlet alkalmazható mint pm i -re azaz: m ✓ ◆ yim 1 + yim 1 X m m yi = = m· · yim m+j 2m 2 j j=0 yi0 = 0 ha (i = 1, 0, 1, 2, ...) vagy (i = n, n + 1, n + 2, ...) Rögzített m-re maxi2Z |yim |=maxi22,...,n+m 1 |yim |. Mivel az elforgatást követően a kezdő és végpont az x-tengelyen van, ebből adódóan az y koordinátájuk nulla, valamint az oszlopok kiegészítése is csupa nullákkal történik így nullánál nagyobb eltérés csak a köztes pontokon valósulhat meg. Az alábbi becslésnél kihasználjuk, hogy emiatt a kezdeti pontsorozatra csak a köztes (n 2) pontra teljesül, hogy |yi0 | > 0. Tehát tetszőleges i = 2, 3, ..., n + m 1 esetén |yim |
m ✓ ◆ 1 X m = m· · yi0 2 j j=0
m m/2 2m
m+j
m ✓ ◆ 1 X m m· · yi0 2 j j=0 M
· (n
z }| 2) · maxk=2,...,n
{
0 1 yk
m m/2 2m
m+j
· (n
m m/2 2m
2) · M
·
m X
yi0
j=0
m+j
m!1
!0 ⌅
Egy p01, ...p0n felosztás mellett, ahol x01 ...x0n az x koordináták vektora vegyük az p01 illetve p0 +p0 az p02 pontokat és végezzünk el egy iterációt azaz kiszámoljuk a 1 2 2 = p12 felezőpontot. Nyilvánvalóan ez a p01 p02 egyenesen lesz. Mivel a kezdőpont rögzített ezért minden iteráp0 +
0 p0 1 +p2
cióra p1 = p01 . A következő iterációban a 1 2 2 = p22 pont lesz a felosztás beli második pont mely szintén az eredeti p1 p2 egyenesre esik. Emiatt megfigyelhető, hogy a felosztás beli 2. pont az iteráció előrehaladtával tetszőlegesen közel kerülhet a kezdőponthoz. Ha kp1 pm 2 k < " elég kis "-ra akkor tekinthetjük őket egy pontnak. Ez a gyakorlati alkalmazás szempontjából igen hasznos, ugyanis számítási kapacitást terheli nagyon közeli 9
pontokra iterálni azonban egy idő után fölösleges mert a pontosságon többet ront, hogy kevesebb iterációt tudunk elvégezni mint amennyit javít. Ekkor az előzőhöz hasonló helyzet adódik most a pm 3 pont fog torlódni. Ezzel látható, hogy a feloszásbeli pontok a kezdő illetve végponthoz torlódnak. Mivel minden iterációban egyel nő a pontok száma ezért ha "-t úgy határozzuk meg, hogy csak nagy időközönként olvasztunk egybe az a számítás pontosságát nem befolyásolja. A végpontnál ugyanez az eljárás alkalmazható. Felmerül azonban, hogy a konvergencia miatt a felosztás beli pontok a köztes szakaszon nem fognak e megritkulni bizonyos számú lépés után. Lemma: Az iterációk során m ! 1 esetén bármely két szomszédos pont x koordinátája tetszőlegesen közel kerülhet egymáshoz, azaz 8" > 0 és 8i esetén xm xm i i+1 < " ha m elég nagy. ⇤ Bizonyítás: Legyen Mmax = maxi22,...,n x0i x0i+1 a legnagyobb távolság két szomszédos x koordináta között. Pm m Pm m Pm m 0 0 0 · x · x x0i+1 m+j ) i m+j i+1 m+j j=0 j j=0 j j=0 j · (xi m+j m xm x = = i i+1 2m 2m Pm
j=0
m j
· x0i
m+j
2m
·x0i+1
m+j
m m/2
· (n
2m
⇤ 1) · Mmax
(m!1)
!0
A fenti becslésnél kihasználjuk, hogy a kezdeti pontsorozat mivel n pontból áll, pontosan (n 1) szakaszból tevődik össze a töröttvonal. Így ha mindnek a hosszát felülről ⇤ becsüljük akkor pontosan (n 1) · Mmax adódik. A Pascal-háromszögből kapott együttm hatókat pedig m/2 -vel majoráljuk. Így beláttuk az állítást. ⌅
Ennek értelmében hiába torlódnak a felosztás beli pontok a kezdő illetve végponthoz, ennek mértéke egy iterációra névze nem túl nagy. Azonban minden lépés során eggyel több pont lesz, ami kompenzálja a jelenséget, így nem ritkulnak meg a pontok a köztes szakaszon. A konvergencia azonban megmarad. Az alábbi ábrákon látható, hogyan működik gyakorlatban az algoritmus. p = [0.1; 0.1], q = [0.9; 0.8] kezdő és végpontok között felvett tetszőleges számú -a példában 28- pontsorozaton végrehajott 10, 100 illetve 1000 iterációra és hogyan konvergál a két pontot összekötő egyeneshez. Megfigyelhető, hogy az iteráció előrehaladtával a kezdő és végpont körül sűrűbb a pontsorozat. 10
A feladat programozása során igen fontos, hogy a pontosorozat vektorának méretét előre lefoglaljuk ún. “preallocated” formában. Az algoritmus folyamán ugyanis mind az x mind az y koordináták vektora iterációnként eggyel nő. Ennek elmulasztásával a számítógépnek minden lépésnél le kell foglalni a memóriaterületet, egy eggyel nagyobb méretű vektor számára, melybe az új pontokat átmásolja. Ez lényegesen rontja az algoritmus futásidejét. Preallocation nélkül a futásidő 30 pontra 1000 iteráció esetén: Az eltelt idő 1.601453 másodperc. Ugyanez előre lefoglalt memóriával: Az eltelt idő 0.096031 másodperc.
(a) kezdeti út
(b) 10. iteráció
(c) 100. iteráció
(d) 1000. iteráció
Figure 1: Középpont módszer síkon
11
Térbeli kiterjesztés Vizsgáljuk meg hogyan tudnánk ezt a módszert alkalmazni felületeken. Legyen 0 r : U ! M sima felület. Adott 0 = (p01 , ..., p0n ⇢ U ). Világos, hogy r által 3 meghatározott pontok összekötéséből származó felezőpontok R -ben nem feltétlenül esnek a felületre. Ezen pontok pusztán felületre történő vetítése nem vezet eredményre hiszen lehet egy adott felezőpontban egy nagy függvényérték amivel rosszabb utat kaphatunk mint az eredeti. Legyen 0 m k. Kezdetben m = 0, k az iterációk száma. Nézzük a következő algoritmust: •
A paramétertartománybeli pontokat vetítsük fel a felületre azaz vegyük r Legyen az így kapott R3 -beli pontsorozat ⌘ = (z10 , ..., zn0 ).
•
8 zi 1 zi szakszra kiszámoljuk a felezőpontot azaz
•
Az így keletkezett R3 felezőpontokhoz meghatározzuk a legközelebbi felületi pontont.
•
A kapott pontokhoz kezdőpontnak hozzáveszük r p01 = z10 míg végpontnak a r p0n = m zn0 pontot. Ezt az utat ⌘ m+1 -nek nevezzük. ⌘ m+1 = (z10 , z2m ...zn+m , zn0 )
•
m értékét egyel növeljük. Ha m < k a második pontra ugrunk egyébként kilépünk a ciklusból
zi
1 +zi
2
0
-t.
i = 2, ..., n + m.
Míg a síkbeli algoritmus futásideje elenyésző még nagy pontsszám esetén is, felületeken már sokkal lassabb a geodetikus görbe meghatározása. Ennek oka, hogy az algoritmus harmadik pontjában szereplő lépéshez a MATLAB-ba beépített f minsearch függvényt használjuk. Ez jelentősen megnöveli a futásidőt. Minél több pontból áll a görbénk közelítése annál többre kell minimalizálni. A futásidő csökkentésére a fent említett "-os egybeolvasztás használható. Ennek használatára később visszatérünk, most vizsgáljuk meg, hogyan működik az algoritmus általános esetben.
12
(a) kezdeti út
(b) 20. iteráció
(c) 100. iteráció
(d) felülnézet
Figure 2: Középpont módszer függvénygrafikonon
A fenti ábrát vizsgálva felmerül a kérdés, hogy ez módszer képes e lokálisan a legrövidebb utat megtalálni azaz kikerülni nagyobb akadályokat, például egy mély gödröt vagy magas hegycsúcsot. Nyújtsuk meg a felületet, hogy a hegyek és a a gödrök sokkal nagyobbak legyenek. Így már világos, hogy a fenti példához hasonlóan nem vezethet út ezeken keresztül. Ugyanazzal a kezdeti pontsorozattal nézzük, milyen utat ad az algoritmus.
13
(a) felület
(b) felülnézet
Figure 3: Középpont módszer függvénygrafikonon
Nem megfelelő kezdeti pontsorozat esetén azonban előfordulhat, hogy a módszer nem konvergál. Vegyük példának az alábbi ábrát. Ez ugyanaz a felület mint a fenti példában, valamint a két végpont is megegyezik. Az egyetlen különbség a felvett kezdeti pontok. Itt ⇥-el az első és a második pontot jelöljük, hogy jól látható legyen miért akad meg a konvergencia.
14
(a) 200 iteráció után
(b) felülnézet
Figure 4: Középpont módszer függvénygrafikonon
Az első és a második koordinátát összekötő szakasz felezőpontja nyilván nincs a felületen. Ekkor az algoritmus szerint vesszük a hozzá legközelebb eső felületi pontot. De ez éppen egybeesik azzal a két ponttal melyből kiszámoltuk a felezőt. Ekkor az f minsearch kiválasztja az egyik pontot jelen esetben a másodikat. A következő iterációs lépésnél ugyanez játszódik le. Általánosságban ha két az utat alkotó pontra ez a jelenség teljesül, hogy az összekötő szakaszuk felezőpontjához eső legközelebbi pontok önmaguk, azon a szakaszon nem lesz konvergens az eljárás. Ezt könnyen kiküszöbölhetjük ha más kezdeti pontsorozatot választunk. Általánosságban éppen egy rosszat találni nem túl valószínű, a fenti pontsorozat a példa demonstrálásához lett megválasztva. 15
Térjünk vissza a futásidőre. A következő példában szerepő felületen a kezdeti pontsorozat 80 pontból áll. Ez azt jelenti, hogy 500 iterációs lépés esetén 580 ponból áll majd a geodetikus görbét közelítő töröttvonal. Mivel tudjuk, hogy a végpontokban sűrűbb a pontsorozat megfelelő számú lépés után, azokra elvégezni a számítást már nem érdemes. Ezért először minden 50-edik lépésben az első 10 és az utolsó 10 felosztásbeli pontot töröljük azaz összesen 20 pontot. Hasonlítsuk össze törléssel - Midepsz2() - és törlés nélkül - Mid() - az algoritmust 50, 100, 200 és 500 lépésre.
k´ ek narancs s´ arga lila 50 100 200 500 M id() 13.414823 mp 31.094625 mp 82.187663 mp 369.293678 mp M idepsz() 12.971356 mp 28.850147 mp 71.239822 mp 272.972471 mp (⇤) 186.385776 mp
16
(a) Mid()
(b) Midepsz()
(c) végpont
(d) 30 pont 40 lépésenként levágva, 500 iteráció (*)
Figure 5: Középpont módszer függvénygrafikonon
A levágás mértékét és azt, hogy milyen időközönként végezzük el befolyásolja a felület és annak meredeksége, a kezdeti pontsorozat mérete és az iterációk száma. Az első és a második ábrán a levágott és az eredeti algoritmus görbéi megegyeznek, így megállapíthatjuk, hogy az eljárás a futásidőt csökkenti és a pontosságot nem rontja. A harmadik ábrán megfigyelhető, hogy a levágást mivel elég finoman végeztük, nem tette pontatlanná a közelítést a végpontokban. Végül 15-15 pont levágása 40 iterációnként még gyorsabb futásidőt - 186.385776 másodperc - eredményez 500 iterációra és így sem látható különbség a pontosságot tekintve. 17
A módszer nem csak függvénygrafikonokon, hanem három dimenziós felületeken is alkalmazható, de itt általános esetben sem a minimalitást, sem azt nem tudjuk garantálni, hogy a görbesorozat egy geodetikus görbéhez konvergál. Egyszerű esetekben az algoritmus használható, például gömb esetén. Ehhez a paramétertartományt le kell cserélni az xy síkról a (0, ⇡) ⇥ (0, 2⇡) nyílt halmazra. Itt veszünk két pontot, majd köztük egy tetszőleges pontsorozatot és ezek képeit a gömbön. Ekkor az előző esettől eltérően nem a felezőpontokhoz legközelebb eső pont Descartes-koordinátáit számítjuk ki, hanem a felezőpont és egy polárkoordinátákkal megadott általános gömbi pont távolságfüggvényének a minimumát. Ehhez is az f minsearch függvényt használjuk. A polárkoordináták meghatározását követően a pontokat ismét a felületre vetítjük és így folytatjuk az iterációt. Gömbi koordináták esetében a Nelder-Mead féle módszer igen lassú konvergenciát eredményez, és nagy lépésszámra az f minsearch alapértelmezett opciókkal nem is tud elég számtítást elvégezni a minimum megkereséséhez. Ezért az options = struct(’MaxFunEvals’, 9999999) beállítást kell átadni a függvénynek, mellyel megemeljük az engedélyezett lépések számát.
18
(a) kezdeti út
(b) 10. iteráció
(c) 100. iteráció
(d) 1000. iteráció
Figure 6: Középpont módszer gömbön
Nézzük meg ennek egy gyakorlati alkalmazását melyet a bevezetőben említettünk. A földgömb mint tudjuk egy ellipszoid a = Egyenl´ıt˝ oisug´ ar(6, 378.1370km) b = P ol´ arissug´ ar(6, 356.7523km) paraméterekkel. Az (IUGG)[9] által használt közelítéssel a föld modellezhető egy gömbbel melynek sugara: 2a + b R1 = = 6, 371.009 kil´ om´ eter 3 19
Figure 7: Budapest-New York
Mivel a gömbfelületen két pont közötti legrövidebb út a rajtuk átmenő főkör rövidebbik szelete, ezért egy hosszabb repülő, vagy hajózási útvonal megtervezésekor ezt az útvonalat célszerű követni. Az egyéb szempontokat mint repülési magasság elérése, országhatárok, időjárási körülmények, melyek szintén befolyásolják egy ilyen útvonal megtervezését, ebben a dolgozatban nem tárgyaljuk. Számítsuk ki ennek segítségével, hogy egy Budapest - New York között haladó gép milyen útvonalat követne. A föld térképének gömbfelületen történő megjelenítéséhez a [10] referenciában szereplő segédprogram használható.
20
Differenciálgeometriai megközelítés Az első fejezetben definiáltuk a felület fogalmát. Definíció: Az M elemi felületnek az r paraméterezéséhez tartozó első főmennyiségein a gij (u1 , u2 ) =< @i r(u1 , u2 ), @j r(u1 , u2 ) > összefüggésekkel meghatározott gij : U ! R (i, j = 1, 2) függvényeket értjük. Az első főmennyiségek u✏U helyen vett függvényértékeiből képezzük a g11 (u) g12 (u) g21 (u) g22 (u)
G(u) =
!
szimmetrikus mátrixot. Vegyünk a U ⇢ R2 paramétertartományban egy y : I ! U ⇢ R2 sima görbét. Az y koordináta függvényei legyenek az yi : I ! R i = (1, 2) függvények. Ezekkel az y leképezés kifejezhető y(t) = (y1 (t), y2 (t)) alakban, ahol a t az I intervallum egy tetszőleges eleme. Egy [a, b] (a, b✏I, a < b) görbedarab ívhossza : v ˆb u 2 2 X uX 0 0 t l( [a, b]) = < gij (y(t))yi (t), yj (t) >dt i=1 j=1
a
Legyen a G(u) mátrix inverze H(u) ahol hij (u) (i, j = 1, 2) szám a H(u) szimmetrikus mátrix i-edik sorának j-edik eleme. Mivel ez a mátrix 2x2-es így könnyen felírható G(u) mátrixból. H(u) =
1 g11 g22
2 g12
"
g22 g21
g12 g11
#
A @1 r, @2 r, N : U ! R3 felületi vektormezők hármasát az M felület r paraméterezéséhez tartozó kísérő Gauss-bázisának mondjuk. Írjuk fel a @11 r, @12 r, @22 r 2 R3 vektorokat az u pont beli Gauss-bázisban: @11 r(u) =
1 11
· @1 r(u) +
2 11
· @2 r(u) + b11 (u) · N (u)
@12 r(u) =
1 12
· @1 r(u) +
2 12
· @2 r(u) + b12 (u) · N (u)
21
@22 r(u) =
1 22
2 22
· @1 r(u) +
· @2 r(u) + b22 (u) · N (u)
Definíció: A lij : U ! R (i, j = 1, 2) függvényeket az M felület r paraméterezéséhez tartozó Christoffel-féle szimbólumoknak nevezzük melyet az alábbi módon számítunk ki: 2
l ij
1 X kl = h (@i gjk + @j gjk 2 k=1
@k gij )
Tétel: A = r y görbe az M elemi felületnek egy geodetikus görbéje akkor és csak akkor, ha a koordináta függvényeire tetszőleges t 2 I esetén fennállnak az 00
yl (t) =
2 X 2 X
l ij
0
0
(y(t)) yi (t)yj (t) (#)
i=1 j=1
(l = 1, 2) összefüggések. Bizonyítás: ([4],8.4). Példa: Nézzük meg a Christoffel féle szimbólumok kiszámítását gömbre és írjuk fel rá a fenti differenciálegyenletet. A gömböt R3 -ban a követekező módon adhatjuk meg paraméteresen: r : (0, 2⇡) ⇥ (0, ⇡) ! R3 , r(u) = (cosu1 sinu2 , sinu1 sinu2 , cosu2 ) 2
3 sinu1 sinu2 cosu1 cosu2 6 7 r0 (u) = 4 cosu1 sinu2 sinu1 cosu2 5 0 sinu1 " # " # 2 sin u 0 1 0 1 2 G(u) = r0 (u)T r0 (u) = ) H(u) = sin2 u2 0 sin2 u2 0 1 Az első főmennyiségek deriváltjai:
•
@1 g11 = 0
•
@1 g12 = @1 g21 = 0
•
@1 g22 = 0 22
•
@2 g11 = 2sinu2 cosu2 = sin(2u2 )
•
@2 g12 = @2 g21 = 0
•
@2 g22 = 0
Számítsuk ki a Christoffel-féle szimbólumokat: 1 11
1 12
=
1 21
= 1 22
2 11
=
=
1 ⇥ 11 h (@1 g11 + @1 g11 2
1 ⇥ 11 h (@1 g21 + @2 g11 2
1 ⇥ 12 h (@1 g11 + @1 g11 2
2 12
=
2 21
=
⇤ @2 g21 ) = ctgu2
@1 g21 ) + h21 (@2 g12 + @1 g22
1 ⇥ 11 h (@2 g21 + @2 g12 2
@1 g22 ) + h21 (@2 g22 + @2 g22
@1 g11 ) + h22 (@1 g12 + @1 g21
1 ⇥ 12 h (@1 g21 + @2 g11 2 1 ⇥ 12 h (@2 g11 + @1 g12 2
⇤ @2 g11 ) = 0 ⇤ @2 g12 ) = ctgu2
@1 g12 ) + h21 (@1 g22 + @2 g21
1 ⇥ 11 h (@2 g11 + @1 g12 2 =
@1 g11 ) + h21 (@1 g12 + @1 g21
⇤ @2 g22 ) = 0
⇤ @2 g11 ) = sin(2u2 )
@1 g12 ) + h22 (@1 g22 + @2 g21 @1 g21 ) + h22 (@2 g12 + @1 g22
⇤ @2 g12 ) = 0 ⇤ @2 g21 ) = 0
⇤ 1 ⇥ 12 h (@2 g21 + @2 g12 @1 g22 ) + h22 (@2 g22 + @2 g22 @2 g22 ) = 0 2 Helyettesítsük be az egyenletbe. A felületi görbe r y, ahol y = (y1 , y2 ) 00 0 0 0 0 l = 1 eset : y1 (t) + 0 · y1 (t)2 + 2ctg(y2 (t)) · y1 (t)y2 (t) + 0 · y2 (t)2 = 0 2 22
=
+ 00
0
0
y1 (t) + 2ctg(y2 (t)) · y1 (t)y2 (t) = 0 23
00
0
0
0
0
l = 2 eset : y2 (t) + 2sin(y2 (t)) · cos(y2 (t)) · y1 (t)2 + 0 · y1 (t)y2 (t) + 0 · y2 (t)2 + 00
0
y2 (t) + sin(2 · y2 (t)) · y1 (t)2 = 0 Most megutatjuk, hogy (1), (2)-nek a főkörök megoldásai. Mivel a gömbi metrika forgásszimetrikus, azaz a gömb bármely tengelye körül forgatás izometria és az ilyen forgatás bármely főkört bármely főkörbe átviszi. Így ha egy példán be tudjuk lántni, hogy a főkör megoldása a geodetikus egyenleteknek akkor ebből következik, hogy bármely főkör az. Példa 1: xy síkbeli főkör. Legyen u1 2 (0, 2⇡), u2 2 (0, ⇡). A paramétertartományban vegyük a P ( ⇡2 , ⇡2 ), Q( 3⇡ , ⇡ ) pontokat. Az őket összekötő görbe paraméterezése y(t) = 2 2 (t, ⇡2 ) ( ⇡2 t 3⇡ ). 2 ⇡ y1 = t, y2 = 2 0 0 y1 = 1, y2 = 0 y1” = 0, y2” = 0 0 + 2ctg(0) · 1 · 0 = 0 X (1) ⇡ 0 + sin(2 · ) · 12 = 0 X (2) 2 Példa 2: yz síkbeli főkör. u1 2 ( ⇡2 , ⇡4 ), u2 2 ( ⇡2 , 3⇡ ), y(t) = ( ⇡2 , t) ( ⇡4 t 4 y1 = ⇡2 , y2 = t 0 0 y1 = 0, y2 = 1 y1” = 0, y2” = 0 0 + 2ctg(0) · 0 · 1 = 0 X (1)
3⇡ ). 4
0 + sin(2 · t) · 02 = 0 X (2)
Numerikus megoldások Adott tehát egy M 2 R3 -beli sima felület. A paramétertartományon vegyük a P, Q pontokat. Az r(P ) és r(Q) felületi pontok között egy olyan = r y felületi görbét keresünk amely geodetikus. Ez mint láttuk akkor és csak akkor teljesül, ha a görbe U -beli ősképe kielégíti a (#) differenciálegyenletet. Világos, hogy ez az egyenlet csak 24
a Christoffel-féle szimbólumokban függ a felülettől. A következőkben szereténk ennek a differenciálegyenletnek egy közelítő megoldását keresni, melyhez az [5] cikkben szereplő két módszert részletesebben kifejtük. A paramétertartományban meghatározott két pont közötti görbét és annak deriváltjait numerikus módszerekkel fejezzük ki. Ezek segítségével egy olyan algoritmust hajtunk végre, hogy az iterálás során adódó görbesorozat a differenciálegyenletünk egy megoldásához konvergál. Redukáljuk nullára a (#) differenciálegyenleteket. 00
yk (t) +
2 X 2 X i=1 j=1
k ij
0
0
(y(t)) · yi (t) · yj (t) = 0
(k = 1, 2)
Legyen N egy pozitív szám. Ezzel az I = [0, 1] intervallumot N egyenlő részre osztjuk és legyen két felosztásbeli pont távolsága ". Így egy felosztásbeli p-edik pont p · " alakban írható fel p = 0, ..., N . A keresett y : I ! U függvényt az y(p") = y p (p = 0, ..., N ) pontokat összekötő y0y1 töröttvonallal közelítjük. Legyen 4ykp = ykp+1 ykp
1
[
y1y2
(k = 1, 2, )
[
...
[
yN
1yN
(p = 1, ..., N
1) .
Szeretnénk a görbe y p pontjában vett deriváltját a centrális osztott differencia felhasználásával közelíteni: y p+1 ykp 1 0 yk (p") = k + O("2 ) 2" 00
yk (p") közelítését ennek segítségével nem tudjuk elvégezni mert a végpontokban nem tudunk centrális, csak féloldali differenciákat meghatározni. Ezért ykp+1 ykp "
ykp ykp "
1
2ykp + ykp 1 yk (p") = + O(" ) = + O("2 ) " "2 Az így kapott közelítéseket behelyettesítjuk a differenciálegyenletekbe. Rendezés után az alábbi egyenletrendszert kapjuk: 00
2
25
ykp+1
ykp+1
2ykp + ykp 2
1
2 1 X + · 8 i,j=1
k p ij (y )
· (yip+1
yip 1 ) · (yjp+1
yjp 1 ) = 0
(?)
(k = 1, 2) (p = 1, ...., N 1) . Továbbá a peremfeltételekből adódik még négy egyenlet: (y10 , y20 ) = P = (P1 , P2 ) (y1N , y2N ) = Q = (Q1 , Q2 ) Ez összesen 2N + 2 egyenlet az y10 , y20 , y11 , y21 , ..., y1N , y2N 2N +2 ismeretlenre nézve. Ezt a 2N +2 dimenziós vektort reprezentálhatjuk az alábbi 2(N + 1) méretű mátrixal is: k=1 k=2
"
y10 y11 ... y1N y20 y21 ... y2N
1 1
y1N y2N
#
Fixpont iteráció A (*) egyenletrendszert alakítsuk át fixpont keresési feladattá: y10 = P1 y20 = P2
ykp
ykp+1 + ykp = 2
1
2 1X + 8 i,j=1
p p k p ij (y )4yi 4yj
(k = 1, 2, p = 1...N
y1N , = Q1 y2N = Q2 Jelöljük F -el az a leképezést amely az x = (y10 , y20 , y11 , y21 , ..., y1N , y2N ) 2 R2N +2 26
1)
vektorhoz az egyenletrendszer jobb oldalán álló vektort rendeli. Ezzel az x = F (x) fixpont keresési feladathoz jutottunk. Ennek megoldásához felhasználhatjuk a Banachféle fixpont tételt. Definíció: Azokat a metrikus tereket, amelyekben minden Cauchy-sorozat konvergens, teljes metrikus térnek nevezzük. Igazolható, hogy RN teljes metrikus tér. Definíció: Legyen (X, d) egy metrikus tér, f : X ! X leképezés. Azt mondjuk, hogy f kontrakció, ha igaz a következő: 9 q 2 [0, 1) 8x1 , x2 2 X : d(f (x1 ), f (x2 )) q · d(x1 , x2 ) Nyilvánvaló, hogy a kontrakciók egyenletesen folytonosak. Tétel:(Banach-f e´le f ixpontt´ etel) Legyen X teljes metrikus tér és f : X ! X kontrakció. Ekkor 1.
Egyértelműen létezik f -nek fixpontja, azaz létezik egyetelen olyan x⇤ 2 X, amelyre f (x⇤ ) = x⇤ .
2.
Tetszőleges x0 2 X kezdőpont esetén az xn+1 := f (xn ) iteráció konvergens, és lim xn = x⇤ . n!1
3.
Érvényes az alábbi hibabecslés: d(x⇤ , xn )
qn 1
q
· d(x1 , x0 )
(n = 1, 2, 3...)
Térjünk vissza a fixpont keresési feladathoz. Az F iteratív alkalmazásával egy olyan R2N +2 -beli vektorokból álló x(i) sorozatot kapunk mely konvergál, egy a P és Q pont között vezető geodetikus görbe egy N + 1 pontból álló közelítéséhez. Az iterációt egy N x(0) = a = (a01 , a02 , a11 , a12 , ..., aN 1 , a2 ) kezdeti vektorral indítjuk. Az iteráció n-edik lépése x(n + 1) = F (x(n)), részletesebben: y10 (n + 1) = P1 y20 (n + 1) = P2 27
ykp (n+1)
2 ykp+1 (n) + ykp 1 (n) 1 X = + 2 8 i,j=1
p p k p ij (y (n))4yi (n)4yj (n)
(k = 1, 2, p = 1...N 1)
y1N (n + 1), = Q1 y2N (n + 1) = Q2 n = 0, 1, 2, 3... Vázlatosan: x(0) = a #F x(1) #F .. . #F x⇤ Ez a sorozat konvergens lesz, ha maxp=1...N |ap ap 1 | megfelelően kicsi. Szükséges még, hogy az egyes koordináták megfelelően közel legyenek F egy fixponjához, azaz a kezdeti paramétertartománybeli görbe, melyet a-val közelítünk, elég közel van egy olyan görbéhez mely kielégíti a differenciálegyenletet. Amennyiben x(n) konvergens a fixpont tétel miatt tudjuk, hogy ekkor a határértéke pontosan F fixpontja ami az eredeti differenciálegyenlet egy megoldását szolgáltatja. Az [5] cikkben tárgyalt példák azt mutatják, hogy a módszer jól alkalmazható síkban és tórusz felületén geodetikus görbék keresésére.
Newton -Raphson módszer A Newton-Raphson módszert először bemutatjuk egyismeretlenes egyenlet gyökének megkeresésére. Legyen f (x) = 0, ahol f : [a, b] ! R és f differenciálható. Az x0 kezdeti megoldásból úgy határozzuk meg x1 -et, hogy kiszámoljuk a függvény grafikonjához az (x0 , f (x0 )) pontban húzott érintő metszéspontját az x tengellyel. Ennek egyenlete: 28
0
f (x0 ) + f (x0 ) · (x x1 = x = x0
x0 ) = 0 f (x0 ) f 0 (x0 )
Hasonlóan nyerhető xn -ből xn+1 xn+1 = xn
f (xn ) f 0 (xn ) 0
Ez a módszer használható f : RN ! RN tipusú függvények esetében is, de ekkor f (xn )0 nel való osztás helyett az f (xn ) 2 RN ⇥N Jacobi-mátrix inverzével való szorzás szerepel. xn+1 = xn
h
0
f (xn )
i
1
· f (xn )
(n = 0, 1, 2, ...)
ahol x0 2 RN kezdeti megoldás. Mivel egy f : RN ! RN típusú f függvény esetén az f (x) = 0 vektoregyenlet egyenértékű az f1 (x1 , ..., xN ) = 0 .. . fN (x1 , ..., xN ) = 0 skalár egyenletrendszerrel, így eljutottunk a nemlineáris egyenletrendszerek NewtonRaphson módszeréhez. Ezt alkalmazzuk a (?) egyenletrendszer megoldására. Ehhez a peremfeltételeket is nullára kell redukálni: ykp+1
2ykp + ykp 2
(k = 1, 2, p = 1...N (y10 , y20 )
1
2 1 X + · 8 i,j=1
k p ij (y )
· (yip+1
yip 1 ) · (yjp+1
yjp 1 ) = 0
1). (y1N , y2N )
(P1 , P2 ) = (0, 0)
(Q1 , Q2 ) = (0, 0)
Azonban az így kapott F (y10 , y20 , y11 , y21 , ..., y1N , y2N ) = (0, ..., 0) egyenletre közvetlenül nem alkalmazható a Newton-Raphson módszer mivel a peremfeltételek állandósága miatt a Jacobi-mátrixban nulla sorok lépnek fel így nem létezik inverze. Ezért 29
a peremfeltételeknek megfelelő egyenleteket kiküszöböljük és csak a belső egyenleteknek megfelelő (2N 2) ⇥ (2N 2) méretű egyenletrendszert oldjuk meg. Az új egyenletrendszer : ykp+1
ykp+1
2 2ykp + Pk 1 X + · 2 8 i,j=1
2ykp + ykp 2 Qk
1
2 1 X + · 8 i,j=1
2ykp + ykp 2
1
k p ij (y )
p+1 k p ij (y )·(yi
2 1 X + · 8 i,j=1
k p ij (y )
· (yip+1
Pi ) · (yjp+1
Pj ) = 0 (p = 1)
yip 1 )·(yjp+1 yjp 1 ) = 0 (p = 2, ..., N 2)
· (Qi
yip 1 ) · (Qj
yjp 1 ) = 0 (p = N
(⇤⇤)
1)
(k = 1, 2) Látható, hogy az eredeti egyenletrendszerből a peremfeltételeknek megfelelő egyenleteket elhagytuk valamint a p = 1 és a p = N 1 egyenleteket a peremfeltételek figyelembevételével módosítottuk. Ezzel egy (2N 2) ⇥ (2N 2) méretű egyenletrendszert kaptunk. Jelölje H : R2N 2 ! R2N 2 azt a leképezést amely az (y11 , y21 , y12 , y22 , ..., y1N
1
, y2N
1
) 2 R2N
2
vektorhoz a (⇤⇤) egyenletrendszer bal oldán álló (2N 2) méretű vektort rendeli. Ezzel a Newton-Raphson módszer egy iterációja így írható fel: x(0) = s 2 R2N # ⇥ 0 ⇤ x(1) = x(0) H (x(0)) # ⇥ 0 ⇤ x(2) = x(1) H (x(1)) # .. . # x⇤ 30
2
1
1
· H(x(0)) · H(x(1))
Az iteráció egy lépésének számolásakor, az inverz mátrix kiszámítása helyett gyorsabb számítást tesz lehetővé ha egy lineáris egyenletrendszert oldunk meg mivel a ⇥ 0 ⇤ 1 z = x(i) H (x(i)) · H(x(i)) vektor kiszámítása ekvivalens a h
i 0 H (x(i)) · (x(i)
z) = H(x(i))
lineáris egyenletrendszer megoldásával. Az algoritmus konvergenciájához szükséges, hogy a H Jacobi mátrixa már nem szinguláris. Állítás: H determinánsa nem nulla. Bizonyítás: Tekintsük az eredeti problémának egy egyszerűbb esetét. Tegyük fel, hogy a felület lokálisan sík, tehát a Christoffel-féle szimbólumok azonosan nullák. kij ⌘ 0. Nézzük, ekkor az egyenletrendszert: yk2 ykp+1 Qk
2yk1 + Pk = 0 , yk2 2
2ykP + ykp 2
2yk1 =
Pk
p=1
1
2ykp + ykp 2
= 0 , yk2
2yk1 + ykp
1
=0
p = 2, ..., N
2
1
=0,
2ykp + ykp
1
=
Qk
p=N
1
k = 1, 2 A Jacobi-mátrix felírásához az ismeretleneket tartalmazó vektort most partícionáljuk k szerint.
31
" 2 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 4
k=1
k=2
z }| 1 2 y1 , y1 , ......, y1N 2
1 .. .
1
2
0
0
1 ..
···
···
.
1
0
2
2
, y1N 0 .. .
{ z }| 1 1 2 y2 , y2 , ......, y2N
, y2N
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
1
0
0
0
0
0
0
0
0
0
0
1
0
···
0 .. .
0
1
0
0
0
0
0
0
0
0
0
1
0
0
0
0
0
0
0
0
0
0
0 .. .
0
0
0
0
0
0
0
2
2
2
2
1 .. . 1
···
0
0 2 1
1 2
3
{ 1
2 7 76 76 76 76 76 76 76 76 76 76 76 76 76 76 76 76 76 76 76 76 74 5
# P1
3
7 7 7 7 7 7 0 7 7 7 Q1 7 7 P2 7 7 7 0 7 7 .. 7 . 7 7 7 0 5 Q2 0 .. .
Ennek az egyenletrendszernek a mátrixa ún. blokk-diagonális mátrix. Determinánsa egyenlő a blokkok determinánsának szorzatával. Ebből egy blokk az 2 3 2 1 0 0 0 0 ··· ··· 0 6 7 2 1 0 0 0 ··· ··· 0 7 6 1 6 7 6 0 1 2 1 0 0 ··· ··· 0 7 6 7 6 7 0 1 2 1 0 0 7 6 0 6 . .. 7 ... ... ... . An = 6 . 7 6 . 7 6 7 6 0 0 1 2 1 0 0 7 6 7 6 0 ··· ··· 0 0 1 2 1 0 7 6 7 6 7 0 0 1 2 1 5 4 0 ··· ··· 0 0 ··· ··· 0 0 0 0 1 2 n ⇥ n-es tridiagonális mátrix (n 2 N+ ). Nézzük a következő ciklust: •
f or i = 1 : n
1.
•
A mátrix i-edik sorának
•
end
i -szeresét i+1
hozzáadjuk az i + 1-edik sorhoz.
32
Példa: 2 2 6 6 0 6 6 0 4 0
4 ⇥ 4-es mátrix esetében a 2. és 3. 3 2 1 0 0 2 7 6 3 2 7 6 1 0 7 s3 := 3 ·s2 +s3 6 0 2 !6 1 2 1 7 5 4 0 0 1 2 0
iterációs lépés 1 3 2
0 0
0 1 4 3
1
0 0 1 2
3 7 7 7 7 5
2
2 6 6 0 !6 6 0 4 0
s4 := 34 ·s3 +s4
1 3 2
0 0
0 1
0 0 1
4 3
5 4
0
Ezzel az n⇥n-es tridiagonális részmátrixból felső háromszög mátrixot kapunk melynek főátlóbeli elemei rendre: 2, 32 , 43 , 54 , ..., n+1 . Egy mátrix determinánsa sorok konsn tansszorosának egymáshoz adásától nem változik. Egy felső háromszög mátrix determinánsa a főátlóbeli elemek szorzata. Végezzük el a szorzást: ( 2) ⇥ (
3 4 5 n+1 ) ⇥ ( ) ⇥ ( ) ⇥ ... ⇥ ( ) 2 3 4 n
Látható, hogy minden nevező a tőle jobbra álló számláló miatt kiegyszerűsödik így a determináns: detAn = ( 1)n (n + 1) A Jacobi mátrix determinánsa e két részmátrix determinánsának a szorzata azaz: detA2n = (n + 1)2 6= 0 Igazoltuk, hogy ha a felület lokálisan sík akkor alkalmazható a Newton-Raphson módszer. Be kell látni továbbá, hogy a Christoffel szimbólumok megjelenése nem rontja el az invertálhatóságot. A 4y p különbségek sűrű felosztás esetén kicsik. Ha a hármas szorzatot deriváljuk a szorzatban mindig lesz egy 4 különbség amely független a változótól. kij (y p ) minden pontban függ a felülettől és korlátos. Ezért a felosztást úgy kell megválasztani, hogy 4 elég kicsi legyen. Tudjuk, hogy a szinguláris mátrixok nullmértékű zárt halmazt alkotnak. Ha egy mátrix elemeit kicsivel változtatom akkor a determináns is csak kicsit változik. Mivel az eredeti determináns (n + 1)2 és az invertálható mátrixok nyílt halmazt alkotnak, azaz van olyan környezetük, hogy minden ott található mátrix invertálható. Tehát ha a mátrix elemeinek perturbációját veszem invertálható marad. Így megfelelő felosztás esetén H elemei csak kicsit változnak és ebből detH 6= 0 következik. ⌅
33
3 7 7 7 7 5
Gráfelméleti megközelítés Az eddigi módszerek alkalmazásakor minden feladatnál arra építettünk, hogy a felületen haladó geodetikus görbét egy töröttvonallal közelíthetjük. Ennek hibája a pontszám növelésével elméletileg tetszőlegesen kicsire csökkenthető, ugyanakkor a gyakorlati alkalmazásban a véges memóriakapacitás ennek határokat szab. Ebben a fejezetben arra keresünk megoldást, hogy magának a felületnek mutassunk egy olyan közelítését, mely lehetővé teszi gráfelméleti módszerek bevonását. A differenciálegyenletek megoldásánál ugyanis mind a Fixpont-módszer, mind a Newton-módszer esetében feltettük, hogy rendelkezésünkre áll egy olyan kezdeti pontsorozat mely elég közel van egy geodetikus görbéhez, így e helyett az iteráció során kapott pontok konvergálnak. Mivel a Középpont módszer - ugyan a gyakorlati szempontjából erre megfelelőnek bizonyult - konvergenciája nyitott probléma, így az esetleges hibák kiküszöbölése végett egy olyan algoritmusra van szükség melyet kezdeti út nélkül el tudunk indítani és a geodetikus görbének egy jó közelítését adja. Felületeket közelíthetünk egy a felületen elhelyezett rácsháló segítségével. Így a feladat visszavezethető a gráfelméletből ismert legrövidebb út keresési problémára. A következőkben megoldást mutatunk arra, hogyan készíthető el a felület ilyen módon való közelítése és miképp alkalmazhatók erre különböző gráfelméleti algoritmusok. Legyen G(V, E) gráf és legyen értelmezve w : E ! R éleken értelmezett súlyfüggvény. Vegyünk a paramétertartományon egy négyzetrács felosztást: ⇢ 1 2 V = 0, , , ..., 1 N N
2
⇢ I 2,
N 2 N+
E = {(u, v) ahol u e´s v szomsz´ edosak} Ennek pontjait felvetítjük a felületre. Ekkor a rácsbeli (u, v) élhez hozzárendeljük az kr(u) r(v)k euklideszi távolságot. A felület közelítésének pontosságát növeli, minél “sűrűbb” hálót illesztünk rá. Így egy síkgárf segítségével reprezentáltuk a felületet. Legyen G = (V, E) gráf és w : E ! R a rajta értelmezett súlyfüggvény. Ekkor s, t 2 V pontok közötti utat jelölje a (v1 , v2 , ..., vk ) pontsorozat ahol s = v1 és t = vk és {vi , vi+1 } 2 E, i = 1, ..., k 1, k 2 N . A legrövidebb út egy olyan pontsorozat, melyre
34
a következő összeg minimális k 1 X
w(vi , vi+1 )
i=1
Így a megfeleltetés kölcsönösen egyértelmű. Ha a paramétertartományban fekvő speciális rács gráfon meg tudunk határozni egy, a két pont között vezető minimális költségű utat, akkor ennek segítségével tudunk közelíteni a felületen levő két pont közötti minimális utakat. Először vizsgáljuk meg milyen ismert algoritmusok alkalmahatók ilyen típusú feladatokra. Szélességi keresés (Breadth First Search) Adott G = (V, E) gráf és válasszunk egy csúcsot ahonnan az algoritmust elindítjuk. Jelöljük az i-ediknek bejárt csúcsot vi -vel. A kiindulási csúcs így s = v0 . Legyen k annak a pontnak a sorszáma, melynek éppen a szomszédait vizsgáljuk. Ahhoz, hogy adott pontba vezető utat is meg tudjuk adni, minden ponthoz rendeljünk egy q változót melyre q(vx ) megmondja mely pontból léptünk vx -be. Az általános algoritmus a gráf minden pontjába mutat egy útvonalat de mi csak egy adott s ! t útra vagyunk kíváncsiak így legyen c´ el = t és k = 0. Nézzük az alábbi algoritmust: 1.
Van e vk -nak olyan szomszédja melyet még nem jártunk be? Igen: 2. pont. Nem: 3. pont.
2.
Választunk egyet a szomszédok közül, ezt jelöljük vi+1 -el és q(vi+1 ) = vk , i = i + 1. (Ha vi+1 = t ! Stop) 1. pontra lépünk
3.
Igaz e, hogy k = i Igen: Stop, Nem: 4. pont.
4.
k = k + 1, 1. pontra lépünk
Itt két index szerint vizsgáltuk a pontokat. Ha adott k-ra nem találunk több szomszédot akkor a feltérképezett pontok közül vesszük a legkisebb indexűt melynek még nem néztük a szomszédait, ez éppen vk+1 . Ha k = i akkor a térképen található utoljára bejárt ponton vagyunk melynek nincs már több szomszédja. A nála kisebb indexűekre már leellenőriztük, hogy van e szomszédjuk de nem volt, így nincs több fel nem térképezett pont. Értelemszerűen, ha idáig eljut az algoritmus, akkor éppen az utoljára bejárt pont volt a cél, máskülönben már megálltunk volna. Ez az algoritmus élsúlyozatlan esetben 35
működik, azaz a súlyozás tekinthető úgy, hogy w ⌘ 1. Ez egy felületet reprezentáló gráfon azt jelenti, hogy a felület lokálisan sík. Az algoritmus működési elve, hogy adott pontból kiindulva ún. “flood fill” azaz árasztásos feltöltés módszerrel térképezi fel a felületet. Ez lehetővé teszi azt is, hogy a lépések során a gráfot képező hálón elhelyezett akadályokat is ki tudjuk kerülni. Dijkstra algoritmus Egy hatékony módszer mely nemnegatív költségfüggvénnyel rendelkező gráfokon két pont között a legrövidebb utat megtalálja Dijsktra algoritmusa. Az előzőektől eltérően itt nemcsak azt tároljuk, mely pontból lépünk be, hanem, hogy mi egy adott pontba vezető út költsége. Az előzőkhöz hasonlóan G = (V, E) gráf, w : E ! R+ az éleken értelmezett költségfüggvény. Az eredeti algoritmus egy adott s pontból meghatározza a minimális út költségét a gráf összes pontjára, de az oda vezető utat nem állítja elő. Ezért minden ponthoz rendelünk egy q változót, valamint adott egy d() tömb melyben a csúcsokhoz az addig talált minimális út hosszát tároljuk. Legyen s = v0 gyökér, ahonnan az algoritmust indítjuk és t az a pont ahova el akarunk jutni. Ekkor d(s) = 0 és d(v) = 1, 8v 2 V s. Legyen S, T 2 V a gráf pontjainak egy felosztása, ahol S-ben tárojuk azokat a pontokat, melyeket az algoritmus során már tovább nem javítunk és T tartalmaz minden más pontot. {s} és legyen u0 := s
1.
S
2.
8(u0 , v) élre ahol v 2 T : ha(d(v) > d(u0 ) + w(u0 , v)) akkor d(v) Ha nincs ilyen (u0 , v) él: 3. pont
3.
A T -beli pontok közül legyen v ⇤ ahol d(v) érték minimális. ha(v ⇤ = t) akkor: Stop
4.
q(v ⇤ ) = u0 , S := S + v ⇤ , T := T
5.
Ha T üres: Stop, különben a 2. pontra lépünk
{s}, T
V
d(u0 ) + w(u0 , v).
v ⇤ , u0 := v ⇤
Ezzel a módszerrel már nem csak w ⌘ 1 költségfüggvénnyel rendelkező élhalmazon tudunk mozogni. Az élköltségek növelésével például reprezentálhatunk egy magas hegyet, vagy egy olyan útszakaszt, melyet szeretnénk elkerülni.
36
Mohó Legjobb Út Keresés (Greedy Best First Search) Az előző két algoritmus hátránya, hogy nem célirányosan a t pontba próbálnak eljutni, hanem minden pontot feltérképeznek és közben találnak rá t-re. Így végeredményben olyan csúcsokra is kiszámítjuk az odavezető legrövidebb utat, melyek az s ! t útnak nem alkotóelemei így fölösleges iterációkat végeztünk. Itt ahelyett, hogy a költségfüggvény figyelembevételével vennénk be a pontokat, azt vesszük be amely legközelebb van a célhoz. Legyen dist(u, t) = ku tk és S egy elsőbbségi tömb, melybe minden beszúrt ponthoz tartozik egy érték. A tömbből mindig csak a legfelső elemet tudjuk kivenni. A beszúrt elemek rendezése a pontokhoz tartozó értékek szerint történik növekvő sorendben. Két azonos értékű pont között a bekerülésük ideje dönt, tehát az új pont kerül feljebb. A tömbbe az S.put(pont, e´rt´ ek) parancsal szúrhatunk be és az S.get() paranccsal vehetünk ki elemeket. Miután egy elemet kivettünk az már nem szerepel a tömbben. T 2 V a még be nem járt pontok halmaza. Kezdetben S.put(s, 0), q(s) =üres. 1.
Ha S nem üres: u0 := S.get() különben: Stop
2.
Ha u0 = t : Stop
3.
8(u0 , v) élre ahol v 2 T :S.put(v, dist(v, t)), 8v-re: {q(v) = u0 , T
v} , 1. pont
Így pontosan abba az irányba indulunk el amerre a cél található, viszont nem vesszük figyelembe a költségfüggvényeket így ez megint leginkább sík esetében használható ahol nincsenek akadályok ott viszont nagyon eredményes. A⇤ algoritmus Egy olyan módszert szeretnénk amely ötvözi az eddigi algoritmusok erősségeit. Találjon egy minimális költségű utat és egyúttal ne végezzen annyi fölösleges iterációt. Erre használható az A⇤ algoritmus. Adott v csúcsra q(v) adja meg honnan léptünk be, és d(v) az addig talált legjobb út költségét a gyökérpontból. Legyen S elsőbbségi tömb az előző példához hasonlóan. Kezdetben q(s) = u¨res, d(s) = 0, d(v) = 1: 8v 2 V s, S.put(s, 0). 1.
Ha S nem üres: u0 := S.get() különben: Stop
2.
Ha u0 = t : Stop 37
3.
Ha van (u0 , v) él ahol v 2 T akkor ujkoltseg = d(u0 )+w(u0 , v), T := T v különben: 1.pont
4.
Ha ujkoltseg < d(v) akkor d(v)
5.
S.put(v, d(v) + dist(v, t)), q(v) = u0 , 3.pont
d(u0 ) + w(u0 , v)
Az A⇤ algoritmus tehát jól alkalmazható a számunkra szükséges felületetet reprezentáló gráfokon való útkereséshez. Egy futás után amikor az algoritmus a Stop állapotba ér és u0 6= t, akkor egyáltalán nem létezik s ! t út. Egyébként a végpontból visszafelé indulva vizsgáljuk, hogy honnan léptünk az adott pontba egészen amíg el nem jutunk s-ig: t, q(t), q(q(t)), ...., s. Ezzel a paramétertartományon kapunk egy pontsorozatot, ahol a geodetikus görbe - mely ebben az esetben a legrövidebb út is - közelíthető ezen pontok felületre vetítése és összekötése által.
38
Összefoglalás A dolgozatban ismertettünk több módszert is geodetikus görbék keresésére, és ezen belül volt, hogy a minimális hosszúságút is megtaláltuk. Az elsőként ismertetésre került Középpont módszer síkban megtalálja a legrövidebb utat, ezt be is bizonyítottuk. Függvénygrafikonok eseténben [1] cikk alapján konvergencia nyitott probléma de mint több példa is mutatta gyakorlatban jól használható. Szintén az [1] cikkben tárgyaltak szerint belátható egy ennél erősebb módszerről, hogy 3 dimenziós sokaságokon biztosan nem működik általános esetben. Egyszerű esetekben példál gömbön alkalmazható az algoritmus és ennek segítségével meg is találtuk azt a felületi görbét amely a két ponton áthaladó főkör rövidebbik íve. Ezt követően a geodetikus görbékre vonatkozó differenciálegyenleten kereszül kerestünk megoldást. Két pontot összekötő görbe paramétertartománybeli ősképét egy töröttvonalal közelítettük és a közelítési pontokban vett deriváltjait centrális és féloldali osztott differenciákkal fejeztük ki. Ezeket a differenciálegyenletbe helyettesítve, egy olyan egyenletrenszert kaptunk melyre ha pontonként teljesülnek az egyenletek, akkor mondhattuk, hogy ennek felületi képe egy geodetikus görbe közelítése. Ennek elérésére két módszer került ismertetésre. Az első a Banach-féle Fixpont tételre vezeti vissza a megoldás keresését, a második a Newton-Raphson módszerre támaszkodik. Mindkét módszernél fel kellett tenni bizonyos feltételeket mely a konvergencia biztosításához volt szükséges. Végezetül egy a függvénygrafikonok közelítésére szolgáló módszerrel gráfelméleti problémára vezettük vissza a feladatot és bemutatásra került, hogyan alkalmazhatók erre ismert gráfelméleti algoritmusok.
39
Függelék Pontok felvétele síkon
Középpont módszer síkon
40
Középpont módszer függvénygrafikonon "-os levágás nélkül.
41
Levágással
Példákban használt függvények
42
Felület ábrázolása
Gömb ábrázolása
43
Középpont módszer gömbön
Programban használt segédfüggvény
44
References [1] Department of Mathematical Sciences, University of Aberdeen, Edward Wright Building, Dunbar Street, Aberdeen AB9 2TY, U.K. - Notes on locally CAT(1) spaces [2] Stoyan Gisbert, Takó Galina - Numerikus módszerek 1 [3] Karátson János előadásai alapján, Kurits Tamás - Bevezetés a funkcionálanalízisbe [4] Verhóczky László - Klasszikus differenciálgeometria [5] Finding geodesics on surfaces - Jongmin Baek, Anand Deopurkar, Katherine Redfield [6] Katona Gyula Y., Recski András, Szabó Csaba - A Számítástudomány alapjai [7] http://theory.stanford.edu/~amitp/GameProgramming/AStarComparison.html [8] http://www.redblobgames.com/pathfinding/a-star/introduction.html [9] http://www.iugg.org/ [10] http://www.mathworks.com/matlabcentral/fileexchange/13823-3d-earth-example
45
