16. Nemlineáris egyenletek megoldása I

1

16. Nemlineáris egyenletek megoldása I. Eddig lényegében lineáris egyenletrendszerek megoldásával foglalkoztunk. De sokszor felvetődik az f ( x) = 0

(16.1)

egyenlet egy (vagy esetleg több) gyökének keresése, ahol az f ( x) ∈ C[a, b] egyváltozós függvény. Minden olyan x∗ érték, amelyre f ( x∗ ) = 0 , a (16.1) egyenlet gyöke vagy f ( x) zérushelye. A gyök az x∗ helyen m -edrendű, ha f ( x) = ( x − x∗ ) m g ( x), g ( x∗ ) ≠ 0 alakban írható. Azzal az esettel foglalkozunk, amikor a megoldás közelítése valamilyen numerikus módszer segítségével végezhető.

16.1. A gyököt tartalmazó intervallum Ha a függvény előjelet vált: f ( a) f (b) < 0 , akkor a folytonosság miatt legalább 1 gyök található [a, b] -ben. Ha létezik f ( x ) első deriváltja is, és előjeltartó [a, b] -ben, akkor csak 1 gyök van. Ha a függvény nem monoton, akkor [a, b] -t célszerű olyan részintervallumokra bontani, ahol az intervallum két végpontja között előjelváltás van. Ílymódon f ( x) páratlan gyökeit el tudjuk különíteni. Deriválható függvény esetén a páros gyököket kereshetjük f ′( x) gyökeiként, mert ekkor a párosakat páratlanná tettük, de kereshetjük f ( x) / f ′( x) gyökeit is, amelyek mind egyszeresek.

16.2. Fixpont iteráció Egy lehetséges eljárás, hogy megpróbáljuk az f ( x) = 0 egyenletet fixpont-egyenletté alakítani: x = F ( x) .

(16.2)

Példa: legyen a megoldandó egyenlet: x 2 − sin( x) = 0 . Ekkor próbálkozhatunk az xk +1 = sin( xk ) iterációval. Fixpont-egyenletet mindig tudunk készíteni, hiszen x = x + cf ( x) is ilyen, ahol c nemzérus állandó, de választhatunk valamely c ( x ) függvényt is olymódon, hogy az iteráció konvergencia-tulajdonságai javuljanak. A fixpont létezéséről szól a

16.2.1 Brouwer fixp onttétel1 Ha F ( x) folytonos [a, b] -ben és F : [a, b] → [a, b] , akkor létezik fixpontja. Bizonyítás. Legyen g ( x) = x − F ( x) , ekkor g (a) ≤ 0 és g (b) ≥ 0 , amiből g (a) g (b) ≤ 0 . Ha itt egyenlőségjel érvényes, akkor már van egy gyök, ha pedig a < jel érvényes, akkor a folytonosság miatt kell léteznie gyöknek [a, b] -ben. ■

16.2.2 Tétel, kontra kció Ha F : S → kontrakció.

folytonosan differenciálható az S zárt intervallumon és F ′( x) < 1, ∀x ∈ S , akkor F

Bizonyítás. A Lagrange középérték tétel alapján x, y ∈ S -re ∃ζ : F ( x) − F ( y ) = F ′(ζ )( x − y ) . Térjünk át az abszolút értékre és használjuk fel, hogy ∃ F ′( x) maximuma S -ben:

1

A tétel többdimenziós megfogalmazása: ha a folytonos F ( x) függvény a gömböt önmagába képezi le, akkor van fixpontja.

2

Hegedűs: Numerikus Analízis II.

F ( x) − F ( y ) ≤ max F ′( x) x − y , x, y ∈ S x∈S

tehát F kontrakció q = max x∈S F ′( x) < 1 kontrakciós állandóval. ■

16.2.3 Következmé ny Ha F ( x) kontrakció, akkor a Banach fixponttétel szerint csak egy gyök van és a kontrakciós állandó ismeretében a közelítés pontosságát is becsülni tudjuk. Visszatérve a fenti példához: a kapott iteráció biztosan konvergens abban a tartományban, ahol cos( x) ′ sin( x) = < 1 . Látjuk, ha x = 0 vagy x = π , ezzel a kifejezéssel baj van, mert a formula 2 sin( x)

(

)

kiértékelésekor 0-val kéne osztani. De például x = π / 4 esetén a kifejezés értéke 1/ 8 , ami már jobbnak tűnik. Ha megrajzoljuk az x 2 parabolát és a sin( x) függvény képét, látjuk, hogy két nemnegatív gyök van, az egyik a zérus, a másik pedig közel x = π / 4 -hez, tehát remélhető, hogy az iteráció π / 4 -gyel indítva konvergens. De az is látszik, hogyha nagyon kicsi pozitív értékkel indítjuk az iterációt, akkor sem kapjuk meg a zérus gyököt, mert az iteráció mindig elvisz a nagyobbik gyök irányába. Ha azonban az x = arcsin( x 2 ) iterációt készítjük, könnyen meggyőződhetünk arról, hogy kis pozitív x -re zérushoz tart. Ha azonban x = 1 -gyel indítunk, akkor először a π / 2 értéket kapjuk, majd komplex számokat, mivel az argumentum nagyobb 1-nél. A tanulság: ügyelnünk kell, a kapott függvény hova képez le, és a leképezés tartományában megmaradnak-e a konvergencia tulajdonságok, illetve azt a gyököt kapjuk-e, amit szeretnénk meghatározni. Ha a megoldandó egyenletben több helyen is szerepel x , akkor több x = F ( x) kifejezés is készíthető. Például szerepeljen két helyen, ekkor meg lehet mutatni: a kapott két iteráció egy adott helyen egyszerre nem lehet konvergens. Legyen ugyanis f ( x1 , x2 ) = 0 a megoldandó egyenlet, ahol a két előfordulást x1 és x2 -vel azonosítjuk. Legyen Fi ( x) az az iterációs függvény, amelyet xi kifejezésével kaptunk. Ez azt jelenti, hogy

f ( F1 ( x), x) = 0

és f ( x, F2 ( x)) = 0.

(16.3)

Legyen α egyszeres gyök: f (α ,α ) = 0 . Ha a (16.3)-ben szereplő kifejezéseket deriváljuk x szerint és helyettesítjük x = α -t, az eredmény: f1′ (α ,α ) F1′ (α ) + f 2′ (α , a ) = 0, f1′ (α ,α ) + f 2′ (α ,α ) F2′ (α ) = 0, ahol f alsó indexe azt jelöli, melyik hely szerint deriváltunk. Ahhoz, hogy egyszeres gyök mellett nemzérus megoldást kapjunk fi′ (α ,α ) -ra kell, hogy a kapott rendszer determinánsa 0 legyen: F1′ (α )

1

1

F2′ (α )

= F1′ (α ) F2′ (α ) − 1 = 0,

ahonnan F1′ (α ) = 1/ F2′ (α ) .

(16.4)

Hacsak nem 1 abszolút értékűek a deriváltak, a gyök közelében az egyik iterációs függvény konvergens, a másik meg divergens lesz. Hasonlóan lehet vizsgálni azt az esetet, amikor x 2-nél

3 többször fordul elő. De ekkor a helyzet rosszabb, az is lehetséges, hogy egyik iterációs függvény sem konvergens. Emiatt azt célszerű tennünk, hogy az x -ek előfordulását két csoportba osztjuk és x -et az egyik csoportból teljesen kifejezük. Például 3 x 2 − 2 x + exp(2.2 x) + 1 = 0 -nél az iterációs függvényre kereshetjük a másodfokú polinom gyökeit úgy, hogy a konstans tag helyére exp(2.2 x) + 1 -et írunk.

16.3. A konvergencia-sebesség Legyen az xn sorozat konvergens, lim n→∞ xn = x* . Jelölje ε = xn − x* az n -edik hibát. Ekkor, ha létezik c állandó és p ≥ 1 szám úgy, hogy p

ε n +1 ≤ c ε n , n = 0,1,…,

(16.5)

akkor az xn sorozat konvergenciája p -edrendű. Ha •

p = 1 , akkor a konvergencia lineáris vagy elsőrendű,

•

1 < p < 2 , akkor a konvergencia szuperlineáris,

•

p = 2 , akkor kvadratikus vagy másodrendű,

•

p = 3 , akkor köbös, vagy harmadrendű.

A p szám jellemzi az iterációs módszer konvergenciájának sebességét. Ha például p = 2 , akkor ez nagyjából azt jelenti, hogy lépésenként az értékes jegyek száma megduplázódik. A fixpont iteráció nem rendelkezik ezzel a sebességgel. Megmutatjuk, hogy p = 1 , azaz a konvergenciája elsőrendű, amennyiben F ′( x* ) ≠ 0 . Ugyanis ε n +1 = xn +1 − x* = F ( xn ) − F ( x* ) ≤ q xn − x* = q ε n .

(16.6)

Ha F ′( x* ) = 0 , akkor a konvergencia magasabb rendű. Erre vonatkozik a következő

16.3.1 Tétel Legyen F valós függvény: F ( S ) ⊂ S ⊂ , S zárt. Tegyük fel, F ∈ C m ( S ) és F ( k ) ( x* ) = 0 , k = 1, 2, , m − 1 . Ekkor az F által meghatározott iteráció konvergencia-sebessége p = m -edrendű. Bizonyítás. Az x* körüli Taylor-polinom m -edrendű maradéktaggal F ( x) = F ( x* ) + F ′( x* )( x − x* ) + … +

F ( m ) (ξ x ) F ( m −1) ( x* ) ( x − x* ) m −1 + ( x − x* ) m (m − 1)! m!

ahol a feltevés szerint az első, második, … , m − 1 -edik derivált eltűnik. Helyettesítsük x = xn -et, vegyük figyelembe, hogy x* = F ( x* ) és xn +1 = F ( xn ) , ezzel xn +1 − x* =

F ( m ) (ξ x ) ( xn − x* ) m , m!

ahonnan ε n +1 =

F ( m ) (ξ x ) m!

xn − x*

m

≤

Mm εn m!

m

, n = 0,1,…

ahol M k = max F ( k ) ( x ) . Innen látható, a konvergencia m -edrendű. ■ x∈S

4


16.4. Newton-iteráció (Newton-Raphson módszer) és a szelőmódszer Ha a függvény első deriváltja létezik a gyök környezetében, akkor a gyököt közelíthetjük úgy, hogy az xn pontban a függvényhez húzott érintő metszéspontját vesszük az x tengellyel. Ez ugyanaz, mint amikor az xn körüli elsőfokú Taylor-polinomot zérussá tesszük és xn +1 -re megoldjuk:

0 = f ( xn ) + f ′( xn )( xn +1 − xn ), innen a Newton-Raphson módszer iterációs formulája: xn +1 = xn −

f ( xn ) = F ( xn ) . f ′( xn )

(16.7)

A szelőmódszert ebből úgy nyerjük, hogy a derivált helyére az utolsó két pontra felírt osztott differenciát tesszük: xn +1 = xn −

f ( xn )( xn − xn −1 ) f ( xn ) xn −1 − f ( xn −1 ) xn = = F ( xn −1 , xn ) , f ( xn ) − f ( xn −1 ) f ( xn ) − f ( xn −1 )

(16.8)

tehát ez az iterációs függvény két pontra támaszkodik. A módszer előnye a Newton-módszerrel szemben, hogy nem kell hozzá a derivált, amit néha eléggé körülményes kiszámítani. Hátránya pedig a kisebb konvergencia-sebesség.

16.4.1 Tétel, a sze lőmódszer hibája Legyen f ( x) ∈ C 2 [ xn −1 , xn , x∗ ] , ekkor a szelőmódszernél az n + 1 -edik iterált hibájára fennáll ε n +1 = ε nε n −1

f ′′(ξ ) , ξ ,η ∈ [ x∗ , xn −1 , xn ] , 2 f ′(η )

(16.9)

ahol x∗ a zérushely és [ x∗ , xn −1 , xn ] az adott pontok által lefedett intervallum. Bizonyítás. Az állítást (16.8)-ból osztott differenciák segítségével származtatjuk. Kihasználjuk, hogy f ( x∗ ) = 0 : ε n +1 = xn +1 − x∗ = xn − x∗ − ( xn − x∗ )

 f ( xn ) − f ( x∗ ) f [ x∗ , xn ]  1 ε 1 = −  = n xn − x∗ f [ xn −1 , xn ] f [ xn −1 , xn ]  

 f [ xn −1 , xn ] − f [ x∗ , xn ] ε n −1  f [ x∗ , xn −1 , xn ] ε ε = εn  =  n n 1 − f [ xn −1 , xn ] xn −1 − x*  f [ xn −1 , xn ]  és innen az osztott differenciák és a deriváltak között érvényes összefüggés (14.1.1 Következmény) segítségével kapjuk az eredményt. ■

16.4.2 Következmén y A Newton-módszerre vonatkozó eredményt az xn −1 → xn határátmenettel kapjuk: ε n +1 = ε n2

f ′′(ξ ) , ξ ∈ [ xn , x∗ ] , 2 f ′( xn )

Látjuk, ha van konvergencia, akkor az másodrendű, feltéve, hogy f ′( x∗ ) ≠ 0 .

(16.10)

5

16.4.3 Tétel, mono ton konvergencia Legyen f ∈ C 2 [ a, b], f ( x∗ ) = 0, x∗ ∈ [ a, b] , az f ′( x), f ′′( x) deriváltak ne váltsanak előjelet [a, b] ben, továbbá az x0 ∈ [ a, b] kezdőpontra teljesüljön f ( x0 ) f ′′( x0 ) > 0 . Ekkor a Newton-módszer konvergens és az általa készített xn sorozat monoton módon tart az x∗ zérushelyhez. Bizonyítás. A Newton-módszer (16.10) formulája szerint az összes iterált a gyöktől vagy jobbra, vagy balra helyezkedik el, mert f ′′ / f ′ előjele állandó. A (16.7) formulából x∗ -ot levonva

ε n +1 = ε n − f ( xn ) / f ′( xn )

(16.11)

következik. Az f ( x0 ) f ′′( x0 ) > 0 feltétel miatt f ′′( x0 ) / f ′( x0 ) és f ( x0 ) / f ′( x0 ) előjele megegyezik. Emiatt, ha (16.10)-ben ε1 > 0 , akkor f ′′( x0 ) / f ′( x0 ) pozitív és (16.11)-ben ε 0 kissebítve van és az összes további lépésben ε n > 0 kisebbedik. Hasonlóan kapjuk, hogy ε1 < 0 esetén az összes további ε n < 0 nagyobbodik, tehát az ε n -ek vagy felülről vagy alulról monoton módon tartanak 0-hoz. ■ Következmény. A (16.9) formula mutatja, hogyha a szelőmódszert úgy indítjuk, hogy x0 , x1 ∈ [a, b] , ε 0 ε1 és f ′′( x0 ) / f ′( x0 ) előjele megegyezik, akkor a tétel feltételei mellett a szelőmódszer is monoton konvergens sorozatot állít elő, mert a formulájában szereplő osztott differencia mindig helyettesíthető egy intervallum-beli deriválttal, aminek az előjele [a, b] -ben állandó.

16.4.4 Tétel, lokális konvergencia Legyen f ∈ C 2 [ a, b], f ( x∗ ) = 0, f ′( x) ≠ 0, x, x∗ ∈ [a, b] , és az x0 ∈ [ a, b] kezdőpontra teljesüljön x0 − x* <

2 min[ a ,b ] f ′( x) 1 = . max[ a ,b ] f ′′( x) M

(16.12)

Ilyen x0 -ból indítva a Newton-Raphson módszer konvergál x* -hoz. A szelőmódszer konvergál, ha x0 mellett x1 is kielégíti a (16.12) feltételt. Bizonyítás. Az első lépéstől kezdve van kontrakció, ha (16.9) vagy (16.10) alapján ε0

max[ a ,b ] f ′′( x) f ′′(ξ ) ≤ x0 − x∗ <1. 2 f ′(η ) 2 min[ a ,b ] f ′( x)

Az állítás innen átrendezéssel adódik. A szelőmódszernél a második lépéshez még ε1 -re is meg kell követelnünk ugyanezt a feltételt. ■ A fentiek alapján a Newton-Raphson módszernél megbecsüljük az n + 1 -edik hibát. Bevezetve a d k = M ε k jelölést d n +1 = M ε n +1 ≤ M 2ε n2

→ d n +1 ≤ d 02

n

→

ε n +1 ≤ ( M ε 0 ) . 2n

(16.13)

16.4.5 Tétel, szelőm ódszer konvergencia-sebessége A 16.4.4 Tétel feltételei mellett az x0 , x1 kezdőpontokból indítva a szelőmódszer p = (1 + 5) / 2 ≈ 1,62 aszimptotikus sebességgel konvergál x∗ -hoz. Bizonyítás. Most ε n +1 ≤ M ε n ε n −1 érvényes (16.9) alapján, ahol M ugyanaz, mint (16.12)-ben. Ismét a d k = M ε k jelöléssel

6


d n +1 ≤ d n d n −1 , n = 1, 2,… Az indításkor x0 − x* < 1/ M és x1 − x* < 1/ M , ezzel d 0 , d1 < 1 . Igaz tehát, hogy ∃d < 1: d 0 , d1 ≤ d , amellyel d 2 ≤ d 2 , d3 ≤ d 3 , d 4 ≤ d 5 , általában d n ≤ d fn , f 0 = f1 = 1, f n +1 = f n −1 + f n , n = 1, 2,… Itt f n -ek a jólismert Fibonacci-sorozat tagjai, melyeknek explicit előállítása ismert: 1

b1n +1 − b2n +1  , b1 = 5

fn =

1+ 5 1− 5 , b2 = . 2 2

(16.14)

Mivel b2 < 1 , a növekvő hatványai zérushoz fognak tartani. Emiatt létezik egy K szám, hogy minden

n -re d sn+1 ≤ K , sn +1 = −b2n +1 / 5 . Tehát írható

(

d n ≤ K d b1 /

5

)

b1n

( )

=K d

b1n

, d = d b1 /

5

.

Kaptuk, hogy a szelőmódszerhez tartozó hibák majorálhatók egy olyan sorozattal, amelynek 1+ 5 konvergenciarendje b1 = ≈ 1,62 , azaz a módszer szuperlineáris. ■ 2

16.5. Példák 1. Legyen f ∈ C 3 [ a, b] . Parabola interpolációval készítsünk három pontra támaszkodó iterációs módszert f ( x) egy [a, b] -beli lokális minimumának meghatározására! Megoldás. Legyen [a, b] -ben három pont ( xi − 2 , f i − 2 ),( xi −1 , f i −1 ),( xi , f i ) . Newton-interpolációval p2 ( x) = f i − 2 + f [ xi − 2 , xi −1 ]( x − xi − 2 ) + f [ xi − 2 , xi −1 , xi ]( x − xi − 2 )( x − xi −1 ) . A deriváltjának zérushelye: xi +1 =

xi − 2 + xi −1 f [ xi − 2 , xi −1 ] − . 2 2 f [ xi − 2 , xi −1 , xi ]

(16.15)

2. Egyenletes lépésközzel haladva hogyan derítenénk fel egy minimumhelyet? Megoldás. Legyen a lépésköz h és x j = a + jh, j = 0,1,… . Az x j −1 , x j , x j +1 alappont-hármas megfelelő, ha f [ x j −1 , x j ] < 0 és f [ x j , x j +1 ] > 0 . Ekkor a lokális minimumot a következő egyszerűsített formulával becsülhetjük, ha (16.15)-ben x j -t vesszük a középső pontnak: xmin ≈

x j −1 + x j +1 2

−

f [ x j −1 , x j +1 ] 2 f [ x j −1 , x j , x j +1 ]

= xj −

f j +1 − f j −1 h . 2 f j +1 − 2 f j + f j −1

(16.16)

3. A kapott iterációs módszerre fogalmazzunk meg ahhoz hasonló tételt, mint amit a Newtonmódszernél láttunk a lokális konvergenciára! Megoldás. Az egyszerűség kedvéért tekintsük (16.15)-ben azt az esetet, amikor i = 2 . Az osztott differenciák tulajdonságait kihasználva a hibák terjedésére próbálunk egy összefüggést származtatni. Vonjuk le mindkét oldalból a minimumhelyet adó x∗ -ot és legyen ε i = xi − x∗ , ezzel ε3 =

ε 0 + ε1 f [ x0 , x1 ] − . 2 2 f [ x0 , x1 , x2 ]

7 A számlálóban lévő osztott differenciát átalakítjuk, kihasználva, hogy f [ x∗ , x∗ ] = 0 és az alappontok sorrendje az osztott differenciákban tetszőleges: f [ x0 , x1 ] = f [ x0 , x1 ] − f [ x1 , x∗ ] + f [ x1 , x∗ ] − f [ x∗ , x∗ ] = = ε 0 f [ x0 , x1 , x∗ ] + ε1 f [ x1 , x∗ , x∗ ] . Beírva a fenti formulába és közös nevezőre hozva: ε3 =

ε 0 ( f [ x0 , x1 , x2 ] − f [ x0 , x1 , x∗ ]) + ε1 ( f [ x0 , x1 , x2 ] − f [ x1 , x∗ , x∗ ]) 2 f [ x0 , x1 , x2 ]

.

A számláló első két tagja továbbírva ε 0ε 2 f [ x0 , x1 , x2 , x∗ ] . A utolsó két tag átalakítása kicsit hosszabb: ε1 ( f [ x0 , x1 , x2 ] − f [ x0 , x1 , x∗ ] + f [ x0 , x1 , x∗ ] − f [ x1 , x∗ , x∗ ]) = ε1ε 2 f [ x0 , x1 , x2 , x∗ ] + ε 0ε1 f [ x0 , x1 , x∗ , x∗ ] .

Ezekkel ε3 =

( ε 0 + ε1 ) ε 2 f [ x0 , x1 , x2 , x∗ ] + ε 0ε1 f [ x0 , x1 , x∗ , x∗ ] . 2 f [ x0 , x1 , x2 ]

2 f [ x0 , x1 , x2 ]

Legyen δ 2 = max { ε 0 , ε1 , ε 2 } és M=

max f (3) ( x) x∈[ a ,b ]

2 min f ′′( x)

.

(16.17)

x∈[ a ,b ]

Felhasználva, hogy az osztott differenciák kifejezhetők a nekik megfelelő rendű deriváltakkal, kapjuk: 3 2! ε 3 ≤ δ 22 2M = δ 22 M . 2 3!

(16.18)

Így ε 3 biztosan kisebb a három megelőző ε abszolút maximumánál, ha δ 2 M < 1 , vagy másképp

δ 2 < 1/ M . Tehát a kapott módszer biztosan konvergens, ha a három induló pont a minimumhely 1/ M -sugarú környezetében van.

16.6. Gyakorlatok 1. Bizonyítsuk be, hogyha f ∈ C1[a, b] , f ( a) f (b) < 0 és f ′( x) nem vált előjelet [a, b] -ben, akkor ott az f ( x) függvénynek csak egy gyöke van. 2. A fixponttétel alkalmazásával mutassuk meg, hogy a cos x − 4 x + 2 = 0, x ∈ egyenletnek egy zérushelye van és x -et 4x felől kifejezve a fixpont iteráció minden kezdőértékre konvergens! 3. Az előző feladatban a gyök milyen környezetéből konvergál biztosan a Newton-iteráció? 4. Oldjuk meg az f ( x) = 1/ x − a = 0 egyenletet Newton-iterációval! Milyen kezdőértékekre van konvergencia? A kapott formula érdekesége, hogy nincs benne osztás, aminek régebben külön jelentősége volt az osztás műveletével nem rendelkező gépi aritmetikákban. 5. Oldjuk meg az konvergenciát!

f ( x) = x 2 − a = 0, a > 0

egyenletet Newton-iterációval és tisztázzuk a

6. Az előző feladat megoldása alapján készítsünk módszert a1/ k meghatározására, ahol a pozitív valós szám. 7. Mutasssuk meg, hogy a 16.4.3 Tételt módosíthatjuk úgy, hogy az f ( x0 ) f ′′( x0 ) > 0 feltételt elhagyjuk és helyette azt követeljük meg, hogy az első lépés után x1 ∈ [a, b] . 8. Mutassuk meg, hogy az F ( xn ) = xn −

( f ( xn ) )

2

f ( xn ) − f ( xn − f ( xn ))

iteráció konvergenciája másodrendű!

8


9. Ellenőrízzük, hogy a Newton-módszer többszörös gyök esetén csak elsőrendben konvergál. 10. Igazoljuk, hogy r -szeres multiplicitású gyöknél a kvadratikus konvergencia megmarad, ha a Newton-módszer formuláját a következőre módosítjuk: xn +1 = xn − rf ( xn ) / f ′( xn ) . 11. Adott ε pontosság elérése érdekében dolgozzuk ki annak feltételét, hogy mikor állítsuk le a Newton-módszert. 12. Mi történik a szelőmódszernél, ha a 16.4.3 Tétel feltételeitől csak annyi a különbség, hogy ε 0 > 0 , de ε1 < 0 ? 13. Mikor állítsuk le a szelőmódszert, hogy a megoldás előírt pontosságú legyen?

17. Nemlineáris egyenletek megoldása II. Néhány speciális esettel folytatjuk.

17.1. Az intervallumfelezés módszere Tegyük fel, az [a, b] intervallum tartalmaz 1 db gyököt: f ( a) f (b) < 0 és a függvény folytonos [ a, b] ben. Az intervallumfelezés módszere szerint ekkor megfelezzük az intervallumot és a két intervallum közül megtartjuk azt, ahol az előjelváltás megmarad. Így az algoritmus: 1.

∃f ∈ C[a, b], f (a ) f (b) < 0 és adott ∈ előírt pontosság.

2. Indulás: [a0 , b0 ] = [a, b], x1 = (a + b) / 2 .  [a , x ], ha f (an −1 ) f ( xn ) < 0, [an , bn ] =  n −1 n 3. [ xn , bn −1 ], egyébként, xn +1 = (an + bn ) / 2. 4. Megállás: ha f ( xn ) = 0 , vagy bn − an <∈. Ez nem túl gyors, de biztos módszer. Az előjelváltásból nem mindig következik a gyök léte. Gondoljunk az 1/ x függvényre, amikor az algoritmust −1 és 2 között indítjuk.

17.1.1 Tétel Az intervallumfelezéssel kapott xn , n = 1, 2,… sorozat elsőrendben konvergens és εn ≤

b−a 2n

, n = 0,1,…

(17.1)

Bizonyítás. A konvergencia abból következik, hogy mindig a gyököt tartalmazó intervallumot tartjuk meg. A hibára minden lépésben teljesül: ε n +1 ≤ ez pedig elsőrendű konvergenciát jelent. ■

1 εn , 2

9

17.2. A húrmódszer (regula falsi) Itt csak annyi az eltérés az intervallumfelezés módszerétől, hogy nem az intervallum közepét vesszük, hanem az ( an , f (an ) ) és ( bn , f (bn ) ) pontokra illesztett egyenes, más névvel: húr zérushelye a következő közlítés: bn − an . f (bn ) − f (an )

xn +1 = an − f (an )

(17.2)

Megfelelő feltételek mellett bizonyítható, hogy a húrmódszer konvergencája lineáris, így nem gyorsabb, mint az intervallumfelezés. Még az is megeshet, hogy annál lassúbb. Ez történik például olyan esetben, amikor a függvény értékei az x -tengelyhez közel vannak és az egyik végponthoz ( an vagy bn ) nagyon közeli a gyök.

17.3. A Newton-iteráció többváltozós esetben Legyen f :

n

→

n

egy n -változós leképezés, amelynek keressük azt a vektorát, amelyre f ( x) = 0 .

Tételezzük fel a differenciálhatóságot, ekkor az xk ∈

n

körüli sorfejtésből közelítve

f ( xk ) + f ′( xk )( x − xk ) = 0 , ahol most f ′( x) =  ∂fi ( x) / ∂x j  ∈

n× n

(17.3)

mátrix – a rendszer un. Jacobi-mátrixa - , amelyről feltesszük,

hogy invertálható. (17.3)-et x -re megoldva a következő iterációt kapjuk: xk +1 = xk − [ f ′( xk )] f ( xk ) . −1

(17.4)

Ha van megoldás és elég közel vagyunk hozzá, akkor remélhetjük, hogy a többváltozós Newtoniteráció konvergens lesz. A módszer azt követeli meg, hogy minden lépésben elkészítsük a deriváltak mátrixát és megoldjunk vele egy lineáris egyenletrendszert. Mivel ez nagyon munkaigényes lehet, szokás alkalmazni a következő egyszerűsítést: Elkészítjük az f ′( xk ) = LU faktorizációt és utána az egyszerűbb xk +1 = xk − ( LU )

−1

f ( xk )

(17.5)

iterációt alkalmazzuk. Ez 1-dimenzióban annak felel meg, hogy lépésenként a derivált értékét nem változtatjuk. Az ilyen módszereket kvázi-Newton módszereknek nevezzük.

17.4. Polinomok gyökei A polinomok gyökeinek meghatározása talán leggyakrabban a mátrixok sajátértékeinek keresésekor jön elő, de ekkor nem érdemes a hatványösszeg alakot használni, mert a lineáris algabrai algoritmusok numerikusan előnyösebb megoldásokat kínálnak. Valóban magasabbfokú polinomok esetén a hatványösszeg reprezentáció n

p ( x) = ∑ ai x i

(17.6)

i =0

nem előnyös, mert gépi számként ábrázolva az együtthatók n növekedésével egyre bizonytalanabb információt nyújtanak a gyökök pontos értékéről. Példának álljon itt Wilkinson kisérlete, aki az 1,2,...,19,20 gyökökkel rendelkező huszadfokú polinomot (17.6) alakban előállította, majd visszaszámolta a gyököket. Az eredmény annyira más volt, hogy több komplex gyökpárt is kapott. A jelenséget magyarázza a gyökök és együtthatók összefüggése: például a nulladfokú tag a gyökök szorzata: 20!, aminek a pontos ábrázolására messzi nem elegendő 15 decimális jegy. Így a gépi számábrázolás folytán sok fontos információ elvész.

10


A (17.6) alak összefüggésbe hozható az un. Frobenius-féle kisérő mátrix-szal, amellyel már találkoztunk a 7.3 szakaszban: 0 0 … … 0 − a0 / an  1 0 … … 0 − a / a  1 n    1  F = .    0 − an − 2 / an    1 − an −1 / an  

(17.7)

1 n ∑ ai λ i . Ennek a an i = 0 mátrixnak ismerete két szempontból is hasznos. Egyrészt mutatja, hogy a polinom xk gyökei lineáris algebrai módszerekkel kereshetők, amelyek a legstabilabbnak tekinthető módszerek közé tartoznak. Másrészt rögtön lehetőségünk van egy olyan körlemez megadására a komplex síkon, amelyben a polinom összes gyöke benne van: Az utolsó oszlopa mentén kifejtve könnyen igazolható, hogy det ( λ I − F ) =

F

∞

= max (1 − δ i 0 + ai / an ) = R ≥ xk , 0≤i < n

ahol δ ij a Kronecker delta. R nagyobb vagy egyenlő F spektrál sugaránál, ami most a polinom gyökök abszolút értékeinek maximuma. Megadhatunk egy másik kisebb körlemezt is, amelyen kívül van az összes gyök. Vezessük be az x = 1/ y transzformációt és írjuk át a polinomot y szerint. Eredményül egy olyan polinomot kapunk, ahol az együtthatók fordított sorrendben vannak és ennek a polinomnak a gyökei az eredeti gyökök reciprokai. Az új polinomhoz tartozó Frobenius-féle mátrix sornormáját véve kapjuk: 1/ xk ≤ max (1 − δ in + ai / a0 ) = 1/ r , ahol természetesen feltételeztük, hogy a0 ≠ 0 . A két eredményt 0
egybevetve látjuk, hogy a polinom gyökei az r ≤ xk ≤ R, k = 1, 2,… , n

(17.8)

körgyűrű tartományba esenek. A (17.6)-tal adott polinomoknál előnyösen alkalmazható a Newton-módszer, mert a polinom értéke és a deriváltja egy lépésben egyszerűen számolható. Ha például a ξ helyen szeretnénk ezeket kiszámolni, nem kell mást tennünk, mint a polinomot maradékos osztással elosztani ( x − ξ ) 2 -tel: p ( x) = q ( x)( x − ξ ) 2 + α x + β .

(17.9)

Könnyen meggyőződhetünk róla, hogy a helyettesítési érték αξ + β , a derivált pedig α lesz a ξ helyen. A többszörös gyökök kiszűrésére alkalmazhatjuk az Euklidészi algoritmust. Ekkor a két induló polinom p0 ( x) = p( x), p1 ( x) = − p′( x) , az i + 1 -edik polinomot pedig úgy készítjük, hogy pi −1 ( x) -et osztjuk pi ( x) -szel és a maradékot képezzük:

pi −1 ( x) = qi ( x) pi ( x) − ci pi +1 ( x), i = 1, 2,…

(17.10)

A sorozatban a polinomok fokszáma csökkenő, ci > 0 , egyébként tetszőleges. Az algoritmus m ≤ n lépés után befejeződik:

pm −1 ( x) = qm ( x) pm ( x), pm ( x) ≠ 0. Az utolsó polinom a két kezdő polinom legnagyobb közös osztója. Mivel a derivált polinom az 1-nél nagyobb multiplicitású gyököket tartalmazza, így ezek a gyökök megjelennek pm ( x) -ben.

11 Abban az esetben, amikor minden gyök valós és egyszeres, akkor az Eulidészi algoritmus olyan polinomsorozatot készít, amely Sturm-sorozat tulajdonságú. Legyen az a helyen a sorozat előjelváltásainak száma V (a ) , a b helyen pedig V (b) , ekkor megmutatható, hogy az [a, b] intervallumban a gyökök száma V (b) − V (a ) .

17.5. Gyakorlatok 1. Készítsük el a (17.9) osztás algoritmusát! 2. Adjuk meg a 4 x 5 − 3 x 4 + 6 x 3 − 5 x 2 − 8 x + 2 polinom gyökeit tartalmazó körgyűrűt!

12


18. Numerikus integrálás (kvadratúra) I. Az integrálok kiszámításakor nem mindig ismert a primitív függvény, vagy ha igen, némely esetben nagyon bonyolult, nehezen számítható. Ilyenkor a numerikus módszerek a kívánt pontosságú eredmény előállítására egyszerűbb alternatívát kínálnak. A továbbiakban az interpolációból nyerhető kvadratúra-formulákkal fogunk foglalkozni. Láttuk, a függvény az [a, b] intervallumban a következő módon állítható elő:

f = Ln + rn ,

(18.1)

ahol Ln a Lagrange-interpolációs polinom és rn a hibatag. (Feltesszük, az alappontok nagyság szerint rendezettek: xi −1 < xi és x0 = a , xn = b .) A kvadratúra-formulák származ-tatási elve: b

∫ a

b

b

n

a

a

i =0

f = ∫ Ln + ∫ rn = ∑ ai f ( xi ) + Rn ,

(18.2)

ahol az b

ai = ∫ li ( x) dx

(18.3)

a

súlyok a Lagrange alappolinomok integrálásából adódnak. Következmény. Az így nyert formulák legfeljebb n -edrendű polinomig pontosak. Ekvidisztáns alappontok esetén nyerjük a Newton-Cotes formulákat.

18.1. Zárt és nyilt Newton-Cotes kvadratúra formulák 18.1.1 Definició Az alappontok halmaza legyen Ω n = {x0 ,…, xn } . Zárt a kvadratúra-formula, ha a, b ∈ Ω n , h = (b − a ) / n , xk = a + k ⋅ h, k = 0,1,…, n . Nyilt a formula, ha a, b ∉ Ω n , h = (b − a ) /(n + 2) , xk = a + (k + 1) ⋅ h, k = 0,1,…, n , x−1 = a, xn +1 = b . A továbbiakban rátérünk a zárt Newton-Cotes formulák együtthatóinak előállítására. b

b

ω n ( x) dx. ( x − xk )ω n′ ( xk ) a

ak = ∫ lk ( x) dx = ∫ a

Vegyük észre: xk − x j = (k − j )h és vezessünk be új változót: t = ( x − a ) / h , ahonnan x = a + th és x − x j = (t − j )h , s ezzel ( t − k ) hiányzik n

t (t − 1) ……… (t − n)h n hdt = n 0 k ( k − 1)…1( −1)( −2)… ( − n + k ) h

ak = ∫

= (b − a )

n−k

(18.4)

n

(−1) t (t − 1)… (t − n) dt = (b − a) Bkzárt ,n , ∫ nk !(n − k )! 0 t−k

ahol a Bkzárt , n együtthatók az intervallumtól függetlenül egyszer s mindenkorra kiszámíthatók. Hasonló módon nyerhetjük a nyilt Newton-Cotes formulák együtthatóit:

13 ak = (b − a)

(−1) n − k (n + 2)k !(n − k )!

n+2

∫ 0

(t − 1)(t − 2)… (t − n − 1) dt = (b − a) Bkny,n , t − k −1

(18.5)

Az első néhány Newton-Cotes együttható: Zárt

Nyílt

1

1

1

4

1

1

3

3

1

7

32

12

32

Trapéz

1

Simpson

1

1

2

-1

2

11

1

1

7

Érintő formula

11

A táblázatban minden sort osztani kell az együtthatók összegével, mert az együtthatók összegének 1nek kell lenni. Például az 1 4 1 súlyok az 1/6 4/6 1/6 valódi súlyokra utalnak.

18.1.2 Tétel n

1.

∑ Bk ,n = 1,

2. Bk ,n = Bn − k ,n .

k =0

(18.6)

Bizonyítás. Az első állítás az f ( x) ≡ 1 függvény integrálásából adódik, kihasználva, hogy az integrál 0-adfokú polinomra pontos. A második állítást az y = n − t új változóra való áttéréssel nyerjük. ■

18.2. Néhány egyszerű integráló formula 1. Az érintőformula (nyílt Newton-Cotes): n = 0 ,

ny B0,0

2

1 = 1 ⋅ dt , tehát 2 ⋅ 1 ⋅ 1 ∫0

a+b I 0 ( f ) = (b − a ) f  .  2 

(18.7)

Az érintőformula úgy is értelmezhatő, hogy a függvényt [a, b] -ben a középponthoz húzott érintő egyenessel közelítjük, és az egyenes alatti területet vesszük. Ez mutatja, hogy legfeljebb elsőfokú polinomig pontos.

18.2.1 Tétel, érintő formula hibája Legyen c = (a + b) / 2, f ∈ C 2 [ a, b] , ekkor az érintő formulával b

∫

f ( x)dx = (b − a ) f (c) +

a

(b − a ) 3 f ′′(η ), η ∈ [a, b]. 24

(18.8)

Bizonyítás. A c körüli sorfejtésből f ( x) = f (c) + ( x − c) f ′(c) +

( x − c) 2 f ′′(ξ x ). 2

Integráláskor az első tag adja a közelítő formulát, a második tag eredménye zérus, így elegendő a hibatagot vizsgálni: b

R1 ( f ) =

1 ( x − c) 2 f ′′(ξ x )dx. 2 ∫a

14


Az integrálszámítás középértéktétele szerint b

b f ′′(η ) f ′′(η )  ( x − c)3  (b − a)3 2 R1 ( f ) = ( x c ) dx f ′′(η ). − = =   2 ∫a 2  3 a 24

■

A gyakorlatban ezt a formulát nem az egész [a, b] intervallumra alkalmazzuk, hanem azt m részre osztjuk, és az egyes részintervallumokban az érintőformulával integrálunk. Például, ha m = 3 : h = (b − a) / 3 és a három részintervallumra alkalmazzuk a (18.7) szabályt. A részintervallumon nyert eredmények felösszegzésével jutunk az érintőszabályhoz: b

∫ a

f =

b−a m (b − a ) 3 f ( a h / 2 ih ) f ′′(η ), − + + ∑ m i =1 24m 2

(18.9)

1 ( f ′′(η1 ) + f ′′(η2 ) + … + f ′′(ηm )) , mert a Darboux-tulajdonság miatt f ′′ ezt az m átlagértéket is felveszi valahol a teljes intervallumban.

ahol most f ′′(η ) =

2. A trapézformula. Elsőfokú polinom interpolációból nyert zárt formula: n = 1, B0z = B1z = 1/ 2 : I1 ( f ) =

b−a ( f (a ) + f (b)). 2

(18.10)

18.2.2 Tétel, trapéz formula hibája Legyen f ∈ C 2 [a, b] , ekkor b

∫ a

f =

b−a (b − a ) 3 ( f (a) + f (b)) − f ′′(η ), η ∈ [ a, b]. 2 12

(18.11)

Bizonyítás. Az interpoláció hibatagjának integrálja az integrálszámítás középérték tételének felhasználásával: b

R1 ( f ) = ∫ a

b f ′′(ξ x ) f ′′(ξ x ) f ′′(η ) ( x − a)( x − b)dx = − ∫ ( x − a )(b − x) dx = − (b − a)3 . 2 2 12 a

■

≥0

A teljes intervallumot m részre osztva, a részintervallumok eredményét felösszegezve nyerjük a trapéz szabályt: b

∫ a

f =

b−a (b − a ) 3 [ f ( x0 ) + 2 f ( x1 ) + … + 2 f ( xm −1 ) + f ( xm )] − f ′′(η ). 2m 12m 2

(18.12)

18.2.3 Definició Egy kvadratúra formulát akkor mondunk k -adrendűnek, ha k -adfokú az a legkisebb fokszámú polinom, amelyre a formula már nem pontos. Eszerint az érintő- és trapézformula másodrendű. 3. A Simpson formula: másodfokú polinom interpolációból nyert zárt formula, n = 2 , B0z = 1/ 6 , B1z = 4 / 6 , B2z = 1/ 6 és I2 ( f ) =

b−a a+b ( f (a) + 4 f ( ) + f (b)). 6 2

15 a+b )( x − b) szerepel, ennek az integrálja [a, b] 2 ben zérus. Ezt legegyszerűbben úgy tudjuk belátni, hogy [a, b] -t a [−1,1] intervallumba transzformáljuk. Ekkor ω 2 ( x) páratlan függvény, amelynek az integrálja zérus. Emiatt a hibatételt az Hermite-interpolációból származtatjuk, Az interpoláció maradéktagjában ω 2 ( x) = ( x − a)( x −

f ( x) = H 3 ( x) +

f (4) (ξ x ) a+b 2 ( x − a)( x − ) ( x − b) , 4! 2

(18.13)

ahol az (a + b) / 2 középpontban az első deriváltat is interpoláljuk. Az általánosított osztott differenciák táblázatára gondolva tudjuk, hogy az interpoláló polinom a következő alakú: a+b H 3 ( x) = L2 ( x) + C ( x − a )( x − )( x − b) . A második tag együtthatójának értéke nem fontos, mert az 2 integrálja az előbbiek alapján zérus, s ezzel

b

b

a

a

∫ H 3 = ∫ L2 .

18.2.4 Tétel, Simps on-formula hibája Legyen f ∈ C 4 [a, b] . Ekkor létezik η ∈ [a, b] , amelyre 5

b

∫ a

f =

(4) b−a a+b  f (η )  b − a  f ( a ) 4 f ( ) f ( b ) + + −     = 6  2 90  2  

= I2 ( f ) −

f

(18.14)

(4)

(η ) 5 h , 90

ahol h = (b − a ) / 2. Bizonyítás. Kiindulunk az Hermite-interpoláció (18.13) alakjából, ahonnan integrálással kapjuk: b

I ( f ) − I2 ( f ) = ∫ a

f (4) (ξ x ) a+b 2 ( x − a)( x − ) ( x − b) dx. 4! 2

Ahhoz, hogy az integrálszámítás középértéktételét alkalmazhassuk, az f (4) mellett álló tényező nem lehet negatív. Ezt úgy biztosíthatjuk, hogy ( x − b) helyett (b − x) -et írunk, s így b

5

f (4) (η ) a+b 2 f (4) (η )  b − a  I ( f ) − I2 ( f ) = − ( x − a )( x − ) ( b − x ) dx = −   . 4! ∫a 2 90  2 

■

Simpson-szabály. A teljes b − a intervallumot páros számú m részintervallumra osztva és a Simpsonformulát a szomszédos intervallumpárokra alkalmazva kapjuk a Simpson-szabályt, mint összetett formulát. Ekkor három pontonként fogjuk össze a formulákat és az összetett formula: b

∫ a

f =

 2h f (4) (ηk ) 5  + + − ( f ( x ) 4 f ( x ) f ( x )) h = ∑  k −1 k k +1  90 k =1,3,…  6 

   h5 h =  f ( x0 ) + 2 ∑ f ( xk ) + 4 ∑ f ( xk ) + f ( xm )  − ∑ f (4) (ηk ). 3 90 k páros k ptlan k ptlan  belső pont   A hibatag még tovább írható:

(18.15)

16

Hegedűs: Numerikus Analízis II. (4) h5 (b − a )5  ∑ f (ηk )  (b − a )5 (4) (4) − f (ηk ) = − = − f (η ),   ∑  90 k ptlan m/2 180m 4  180m 4 

(18.16)

mivel a Darboux-tulajdonság miatt van egy η , amelyre a negyedik derivált az átlagértéket felveszi.

18.3. Példák 1

1) Az

dx

∫ 2+ x

integrált az érintőszabállyal közelítjük. Hány osztópontot kell választanunk, hogy az

−1

integrált 10−2 -nél kisebb hibával kapjuk? Megoldás.

Azt

kell

biztosítani,

(b − a ) 3 M 2 ≤ 10−2 24m 2

hogy

teljesüljön,

b−a=2

ahol

és

M 2 = 2 max (2 + x) −3 = 2 . A számokat helyettesítve: 200 / 3 ≤ m 2 , így m = 9 megfelel. x∈[ −1,1]

2

2)

Határozzuk

meg

az

A0 , A1 , A2

paramétereket

úgy,

hogy

a

∫

x f ( x)dx ≈

0

≈ I 2 ( f ) = A0 f (0) + A1 f (1) + A2 f (2) kvadratúra legfeljebb másodfokú polinomokra pontos legyen! Megoldás. Két megoldás is létezik. Az egyik, hogy kiszámítjuk a kijelölt integrálokat a Lagrange2

alappolinomokkal: Ai = ∫ li ( x) xdx , ahol

x -et súlyfüggvénynek tekintjük. A másik módszer szerint

0

felírjuk azt a lineáris egyenletrendszert, ami a pontossági követeléseket tartalmazza. A következő egyenletrendszer első sora azt fejezi ki, hogy a kvadratúra az 1 polinomra pontos, a második sor szerint az x polinomra pontos, a harmadik sor szerint pedig az x 2 polinomra pontos:  2 x1/ 2 dx = 4 2 / 3  ∫  1 1 1   A0   0   0 1 2   A  = 2 x ⋅ x1/ 2 dx = 8 2 / 5  .     1   ∫0   0 1 4   A2   2 2 1/ 2  x x dx = 16 2 / 7   ∫0  Ennek a megoldása: A0 =

8 2 32 2 12 2 , A1 = , A2 = . 105 35 35

17

19. Numerikus integrálás, Gauss-kvadratúrák II. Az eddigi, interpolációból származtatott kvadratúra-formulák legalább annyiad fokú polinomra pontosak, ahányad fokú polinomból származtattuk őket. A Gauss-kvadratúrák abból az észrevételből származnak, hogy az alappontok speciális megválasztásával a kavadratúra-formula rendje növelhető. Ismét szükségünk lesz az ortogonális polinomokra.

19.1. Tétel, ortogonális polinom gyökei Legyen { pk ( x)} egy ortogonális polinom rendszer. Ekkor bármely n -re pn +1 ( x) gyökei valósak, egyszeresek és az [a, b] intervallumban vannak, ahol [a, b] a skalárszorzat integrálási tartománya. Bizonyítás. Legyenek x0 , x1 ,…, xk pn +1 ( x) páratlan multiplicitású gyökei [a, b] -ben, azaz ott pn +1 ( x) előjelet vált. Ha k = n , akkor a tétel állítása rendben van. Ha nem, akkor indirekt úton feltesszük: k < n és megmutatjuk, hogy az állítás ellentmondásra vezet. Ehhez tekintsük a q( x) = ( x − x0 )( x − x1 )… ( x − xk ) k + 1 -edfokú polinomot. Mivel k + 1 < n + 1 , az ortogonalitás miatt ( pn +1 , q ) = 0 . De ezzel ellentmondásra jutunk, mert pn +1 ( x)q ( x) nem vált előjelet [a, b] -ben, mivel

pn +1 minden előjelválátását q ( x) megszünteti és így

∫ pn+1qα ≠ 0

volna. Vegyük észre, a

gondolatmenet akkor is jó, ha egyetlen páratlan multiplicitású gyök sincs, mert ekkor q ( x) 0 -adfokú. ■ Az n + 1 -pontos Gauss-kvadratúrát úgy kapjuk, hogy a pn +1 ( x) ortogonális polinom gyök-helyein készítjük az interpolációból származtatott kvadratúra-formulát. A séma a következő: b

∫ a

b

b

n

a

i =0

b

f α = ∫ Lnα + ∫ rnα = ∑ ai f ( xi ) + Rn , ai = ∫ li ( x)α ( x)dx. a

(19.1)

a

19.2. Tétel, Gauss-kvadratúra pontossága Legyenek a pn +1 ( x) ortogonális polinom gyökei x0 , x1 ,…, xn , ai = ∫ liα , ahol li az i -edik Lagrange alappolinom a fenti alappontokon. Ekkor a Gauss-kvadratúra n

Gn ( f ) = ∑ ai f ( xi ) i =0

pontos minden legfeljebb 2n + 1 -edfokú polinomra: f ∈ P2 n +1 → ∫ f α = Gn ( f ) . Bizonyítás. Az interpolációból való származatatás miatt Gn ( f ) biztosan pontos a legfeljebb n -edfokú polinomokra. Tegyük fel, f ∈ P2 n +1 , f = pn +1 ⋅ q + r , q, r ∈ Pn , így n

n

i =0

i =0

Gn ( f ) = ∑ ai f ( xi ) = ∑ ai [ pn +1 ( xi ) ⋅ q ( xi ) + r ( xi )] = = 0 minden i − re

n

= ∑ ai r ( xi ) =Gn (r ) = ∫ rα =

(mert n-edfokig pontos)

i =0

= ∫ ( pn +1 ⋅ q + r )α = =∫ f α.

(mert q ∈Pn , így ortogonális pn +1 -re)

18


19.2.1 Következmé ny Az ai együtthatók pozitívak. Bizonyítás.

Tudjuk, li ( x j ) = li2 ( x j ) = δ ij , ahol

δ ij a Kronecker-delta, li2 ( x) ≥ 0 és li2 ( x) ∈ P2 n ,

következésképp a Gauss-kvadratúra pontos : n

0 < ∫ li2α = Gn (li2 ) = ∑ ai li2 ( x j ) = ai . j =0

Az f ( x) = 1 függvény integrálásával most a következőt kapjuk az együtthatók összegére: n

∑ ai = ∫ α = µ0 , i =0

µi = ∫ xiα ,

ahol µ0 a nulladik momentum. Vegyük észre, ez egyenlő b − a -val, ha a súlyfüggvény 1.

19.2.2 Tétel, Gaus s-kvadratúra hibaformulája n

Legyen f ∈ C 2 n + 2 [a, b] és Gn ( f ) = ∑ ak f ( xk ) , ahol az alappontok pn +1 ( x) gyökei. Akkor k =0

I ( f ) − Gn ( f ) =

f (2 n + 2) (η ) ( pn +1 , pn +1 ) , (2n + 2)!

(19.2)

itt pn +1 ( x) 1-főegyütthatós ortogonális polinom. Bizonyítás. Hermite-Fejér interpolációból (amikor az interpolációban a függvényértékek és az első deriváltak vesznek részt) kapjuk a következő hiba-előállítást: f ( x) = H 2 n +1 ( x) +

f (2 n + 2) (ξ x ) ( x − x0 ) 2 ( x − x1 ) 2 … ( x − xn ) 2 . (2n + 2)! 2 = pn +1 ( x )

Innen az integrálás középérték-tételének alkalmazásával b

I ( f ) − Gn ( f ) = ∫ a

f (2 n + 2) (ξ x ) 2 pn +1 ( x)α ( x) dx (2n + 2)! ≥0

nyerjük az állítást, mert H 2 n +1 ( x) -re a Gauss-kvadratúra pontos. ■ Megadunk néhány 1-főegyütthatós ortogonális polinomot: Név

[ a, b]

α ( x)

µ0

α n +1

βn

p0

p1

p2

Legendre

[−1,1]

1

2

0

n 2 /(4n 2 − 1)

1

x

x 2 − 1/ 3

Csebisev

[−1,1]

π

0

1

x

x 2 − 1/ 2

Laguerre

[0, ∞]

e− x

1

2n + 1

n2

1

x −1

x2 − 4 x + 2

Hermite

[−∞, ∞]

e− x

π

0

n/2

1

x

x 2 − 1/ 2

1 1 − x2

2

1/4, de

β1 = 1/ 2

19

19.3. Példák 1

1. Három-pontos Gauss-Csebisev kvadratúrával közelítjük a

2 ∫ (1 − x )

−1/ 2

e − x dx integrált. Becsüljük

−1

meg a hibát! Megoldás. A három-pontos kvadratúránál n = 2 , ezt alkalmazzuk a (19.2) formulánál: M 6 = e és mivel 1-főegyütthatós polinomoknak kell szerepelni, emiatt p3 ( x) = T3 ( x) / 4 . Így a hiba kisebb, mint e ⋅ (T3 , T3 ) /(16 ⋅ 6!) = eπ /(32 ⋅ 720) , mert (T3 , T3 ) = π / 2 , (lásd a 7.4 feladatot). 2. Készítsük el a két-pontos Gauss-Hermite kvadratúra súlyait! Ellenőrízzük, hogy az így kapott kvadratúra legfeljebb harmadfokú polinomokra pontos! Megoldás. A két-pontos kvadratúránál a másodfokú Hermite-polinom gyökei − x0 = x1 = 2−1/ 2 . Így ∞ ( x − x1 ) xµ a0 = ∫ exp(− x 2 ) dx = 1 0 = µ0 / 2 , mert az elsőfokú tag integrálja zérus, lévén az −∞ ( x − x ) 2 x1 0 1 integrandus páratlan függvény. Hasonlóan adódik, hogy a1 = a0 . A kapott kvadaratúra pontos az 1 függvényre, mert az eredmény µ0 . Az x és x3 függvényre is pontos, mert a két tag a gyökhelyeken a páratlanság miatt kiejti egymást. Így már csak azt kell igazolunk, hogy a pontosság a másodfokú x 2 polinomra is teljesül. Ennek az integrálja a (9.8) formula alapján nem más mint ( p1 , p1 ) = β1 ( p0 , p0 ) , ami az előző oldalon látható táblázat szerint egyenlő µ0 / 2 -vel. A Gauss-Hermite kvadratúra µ 1 1 eredménye pedig 0 ( + ) , tehát igaz az egyezés. 2 2 2 3. Határozzuk meg a következő integrál értékét: 1

∫

−1

(2 x 2 + x)dx 1 − x2

.

Megoldás. A számlálóban szereplő polinomot előállítjuk az első három Csebisev polinom lineáris kombinációjaként: 2x 2 + x = c0T0 + c1T1 + c2T2 . Ezt felhasználva az integrálunk így is írható: (T0 , c0T0 + c1T1 + c2T2 ) = c0 (T0 , T0 ) = c0π . Az első három Csebisev polinom együtthatóival a következő lineáris egyenletrendszert írhatjuk fel (9.5) alapján: 1 0 −1  c0   0   1 0   c  = 1  ,   1    2   c2   2  amiből c0 = 1 , tehát az integrál értéke π .

20


20. Közönséges differenciálegyenletek 20.1. Alapfogalmak Az F ( x, y , y ′, y ′′,… , y ( m ) ) = 0

(20.1)

alakú egyenletet közönséges differenciálegyenletnek nevezzük, ahol keresendő az y = y ( x) függvény, mely a fenti kifejezésben deriváltjaival együtt szerepel. Ha y m -szer differenciálható [0,1] -ben és kielégíti (20.1)-et, akkor y egy megoldás. Ismeretes, a differenciálegyenleteknek sok megoldásuk van. A megoldást úgy tudjuk egyértelművé tenni, hogy a peremen még feltételeket teszünk a függvény vagy deriváltjainak értékére. Ha ezen feltételek mindegyike a kezdőpontban (általánosabban fogalmazva: csak egy pontban) van megadva, akkor kezdetiérték feladatról beszélünk, ha pedig több ponthoz kapcsolódnak a feltételek, akkor peremérték feladatunk van. A differenciálegyenlet lineáris, ha y és deriváltjainak lineáris kombinációja szerepel (20.1)-ben és a differenciálegyenlet explicit alakú, ha a legmagasabb derivált explicit módon kifejezhető: y ( m ) = f ( x, y , y ′,…) .

(20.2)

Világos, minden lineáris differenciálegyenlet ebbe a kategóriába tartozik, és a gyakorlatban előforduló nemlineáris differenciálegyenletek zöme is ilyen. A továbbiakban explicit elsőrendű differenciálegyenletek megoldásának numerikus közelítéseivel fogunk foglalkozni, amikor kezdetiérték feladatról van szó. Ekkor a kezdetiérték feladat

y (0) = y0 , y ′ = f ( x, y ), x ∈ [0,1],

(20.3)

függvényt, amely felveszi x0 ahol f : 2 → , x0 és y0 adottak, és keressük azt az y :[0,1] → ban az y0 értéket és az egyenletet kielégíti. Az egyszerűség kedvéért szoritkozunk a [0,1] intervallumra, hiszen tudjuk, az [a, b] intervallum lineáris transzformációval ide átvihető. A megoldás létezésére és egyértelműségére vonatkozik a következő tétel, melyet bizonyítás nélkül idézünk.

20.1.1 Tétel, a kez detiérték feladat egyértelműsége Ha f folytonos egy ( x, y ) ∈ [0,1] × [c, d ] téglán és f a második változója szerint eleget tesz a Lipschitz-feltételnek, azaz létezik L : f ( x, y1 ) − f ( x, y2 ) ≤ L y1 − y2 , ∀x ∈ [0,1], y1 , y2 ∈ [c, d ] ,

(20.4)

akkor a kezdetiérték feladatnak egyértelmű megoldása van. A numerikus eljárás során a megoldást az xn = nh osztópontokban közelítjük, ahol h = 1/ N és a numerikus közelítést xn -ben yn -nel jelöljük. Természetesen azt szeretnénk, hogy yn minél közelebb legyen a pontos értékhez, y ( xn ) -hez.

20.1.2 Az Euler-mó dszer Ez a legegyszerűbb módszer, amit a differencia-hányadosból származtatunk a következő közelítő érték előállítására: yn +1 − yn = f ( xn , yn ) → yn +1 = yn + hf ( xn , yn ). h

(20.5)

21 Példa: y (0) = 1, y′ = y . Ennek a megoldása y = e x és (20.5)-ből yn = (1 + h) n . Ha n → ∞, akkor az analízisből tudjuk, hogy yn = (1 + h) n = (1 + xn / n) n → e xn .

20.1.3 Definíció, ko nvergencia Egy numerikus módszer lokálisan konvergens az x ∈ [0,1] pontban, ha h = x / n, x = xn mellett lim yn = y ( x)

n →∞

teljesül. Ha a módszer konvergens ∀x ∈ [0,1] -re, akkor azt mondjuk, a módszer konvergens.

20.1.4 Definíció, ko nvergencia-sebesség Egy konvergens numerikus módszer sebessége p -edrendű, (1 ≤ p ) , ha h = x / n, x = xn mellett ∃M úgy, hogy az

y ( x) − yn ≤ Mh p

(20.6)

hibabecslés teljesül, ahol M független h -tól és n -től.

20.1.5 Definíció, lo kális hiba v. képlethiba Ez annak a képletnek a hibája, amellyel a következő függvényértéket közelítjük. Ilyenkor a képletbe mindenütt a pontos értéket írjuk. Például az Euler-módszernél a

gi −1 = y ( xi ) − y ( xi −1 ) − hf ( xi −1 , y ( xi −1 )), i = 1, 2,…, N

(20.7)

mennyiségek a lokális hibák vagy képlethibák.

20.1.6 Definíció, ko nzisztencia Ha ∃p ≥ 1 és M konstans úgy, hogy

gi ≤ Mh p +1 , i = 1, 2,…, N

(20.8)

akkor a módszer p -edrendben konzisztens. A definició alapján világos, nagyobb p -re pontosabb megoldás várható.

20.1.7 Definició, gl obális hiba Jelölje a pontos és numerikus megoldás eltérését

ei = y ( xi ) − yi , i = 1, 2,…, N , h = 1/ N , xi = ih.

(20.9)

Az ei mennyiségeket globális hibának nevezzük.

20.1.8 Definició, st abilitás A numerikus módszer stabil, ha van olyan K konstans, amellyel i  ei ≤ K  e0 + ∑ g j  j =1 

  , i = 1, 2.… , N  

teljesül, vagyis a globális hiba felülről becsülhető a lokális hibák abszolút összegével.

(20.10)

22


20.1.9 Tétel, nume rikus módszer konvergenciája Ha egy numerikus módszer stabil és p -edrendben konzisztens, akkor p -edrendben konvergens. Bizonyítás. Ha x = 0 , akkor a tétel igaz. Ha x ∈ (0,1] , legyen h = x / n, xn = x ∀n -re. Először a stabilitás, majd a konzisztencia definicióját felhasználva n  ei ≤ K  e0 + ∑ g j  j =1 

n   ≤ K ∑ ch p +1 ≤  j =1 

≤ Kcnh p +1 ≤ Kc(nh)h p ≤ ( Kc)h p , ahol e0 = y ( x0 ) − y0 = 0. Ezzel p -edrendű konvergencia-sebességre jutottunk.

■

20.2. Az Euler-módszer vizsgálata Látjuk, az Euler-módszerre a konzisztenciát és stabilitást kéne megmutatni ahhoz, hogy a konvergenciát igazoljuk.

20.2.1 Tétel, konzis ztencia Az Euler-módszer elsőrendben konzisztens, vagyis gi = O(h 2 ) teljesül, ha a megoldás kétszer folytonosan differenciálható. Bizonyítás. Tekintsük a lokális hibát az i -edik pontban és y ( xi +1 ) -et fejtsük sorba az xi hely körül: gi = y ( xi +1 ) − y ( xi ) − hf ( xi , y ( xi )) = = y ( xi ) + hy ′( xi ) +

h2 y ′′(ξi ) − y ( xi ) − hy ′( xi ), 2!

emiatt gi =

h2 h2 y ′′(ξ i ) ≤ y ′′ ∞ , 2! 2!

tehát p + 1 = 2 , p = 1 , elsőrendű a konzisztencia.

■

20.2.2 Tétel Az Euler-módszer stabil, ha f eleget tesz a második változója szerint a Lipschitz-feltételnek. Bizonyítás. Vegyük az i -edik pontban a lokális hiba képletét, a módszer képletét és vonjuk ki őket egymásból: y ( xi +1 ) = y ( xi ) + hf ( xi , y ( xi )) + gi yi +1 = yi + hf ( xi , yi )

/(+ ) /(−)

ei +1 = ei + h[ f ( xi , y ( xi )) − f ( xi , yi )] + gi Mivel f eleget tesz a Lipschitz-feltételnek: ei +1 ≤ ei + h f ( xi , y ( xi )) − f ( xi , yi ) + gi ≤ ei (1 + hL ) + gi . Fejtsük vissza a rekurziót e0 -ig:

23 ei +1 ≤ ei (1 + hL ) + gi ≤ ( ei −1 (1 + hL ) + gi −1 ) (1 + hL ) + gi ≤ ≤ (1 + hL ) ei −1 + (1 + hL ) gi −1 + g i ≤ (1 + hL ) 2

i +1

i

e0 + ∑ (1 + hL ) k =0

ebből megkapjuk a kívánt becslést, ha felhasználjuk az (1 + hL ) ≤ e jhL = e j

xjL

i −k

gk ,

≤ e L relációt:

i i   ei +1 ≤ e L e0 + ∑ e L g k = e L  e0 + ∑ g k  . k =0 k =0  

■

A továbbiakban rátérünk néhány fontosabb módszercsalád rövid ismertetésére.

20.3. Taylor-polinomos módszerek Az Euler-módszer nem egyéb, minthogy vesszük a függvény elsőfokú Taylor-polinomját xn körül és annak segítségével lépünk a következő pontba. Ebből az ötletből kiindulva készíthetünk magasabbrendű módszereket is. A következő módszert m -edrendű Taylor-polinomos módszernek nevezzük: yn +1 = yn + y ′( xn ) + … + y ( m ) ( xn )h m / m !, n = 0,1,… , N .

(20.11)

Fontos kérdés, egyáltalán kiszámíthatók-e ezek a deriváltak. Ha f elegendően sokszor differenciálható, ennek elvi akadálya nincs. Sokszor azonban súlyos gyakorlati nehézség, hogy a deriváltak nagyon hosszú, nehezen kezelhető formulákat eredményeznek. A konstrukcióból látható, m edrendben konzisztens módszerre vezet a fenti eljárás.

20.4. Runge-Kutta módszerek Ezek a Taylor-polinomos módszerek fenti nehézségét küszöbölik ki: nem kell magasrendű deriváltakat számolni, a magasabb konzisztencia-rend rekurzív függvényhívásokkal is elérhető. Az általános alak: k1 = f ( xn , yn ), k2 = f ( xn + ha2 , yn + hb21k1 ), (20.12) j −1

k j = f ( xn + ha j , yn + h∑ b jl kl ), j = 1, 2,… , s, l =1

s

yn +1 = yn + h∑ c j k j . j =1

Az előre kiszámolt ai , bij , ci paraméterek határozzák meg a konkrét módszert. Az s szám a rekurzió mélysége, más szóval: az egy lépés megtételéhez szükséges függvényhívások száma. A következő, ún. módosított Euler-módszer Runge-től származik: k1 = f ( xn , yn ), k2 = f ( xn + h / 2, yn + (h / 2)k1 ),

(20.13)

yn +1 = yn + hk 2 . Deriválással megmutatható, hogy másodrendben konzisztens. Hasonlóképp másodrendű a következő Runge-Kutta módszer, amelyet az egyszerűsége miatt egy sorban írunk fel: h yn +1 = yn + [ f ( xn , yn ) + f ( xn +1 , yn + hf ( xn , yn ))]. 2

(20.14)

24


Az igen népszerű negyedrendű Runge-Kutta módszer negyedrendben konzisztens és stabil, vagyis negyedrendű konvergens módszer: k1 = f ( xn , yn ), k2 = f ( xn + h / 2, yn + hk1 / 2), k3 = f ( xn + h / 2, yn + hk2 / 2),

(20.15)

k4 = f ( xn + h, yn + hk3 ), h yn +1 = yn + (k1 + 2k2 + 2k3 + k4 ). 6 A Runge-Kutta módszerek lokális hibája: s

g n = y ( xn ) − y ( xn −1 ) − h∑ c j k j ( xn −1 , y ( xn −1 ), h), n = 1,… , N ,

(20.16)

j =1

ahol k j ( xn −1 , y ( xn −1 ), h) azt jelenti, hogy a k j sorozatot a pontos y ( xn −1 ) értékkel indítva számoljuk végig. Egy adott rend eléréséhez a lokális hibát sorbafejtjük és a paramétereket úgy választjuk, hogy a tagok minél magasabb rendig eltűnjenek.

20.5. Lineáris többlépéses módszerek Ezek is az Euler-módszer általánosításának tekinthetők. Az s -lépéses módszer általános alakja: s

s

k =0

k =0

∑ α k yi + k = h ∑ β k fi + k , fi +k = f ( xi+ k , yi + k )

(20.17)

ahol α k , β k , k = 0,1,…, s adottak, α s = 1 valamint α 0 + β 0 ≠ 0 teljesül. Például az Euler-módszerre

s = 1, α 0 = −1, α1 = 1, β 0 = 1, β1 = 0. A módszer indításához az y0 értékén túl kellenek még az y1 , y2 ,…, ys −1 értékek. Ha β s = 0 , akkor a módszer explicit és könnyen számolható. Ha β s ≠ 0 , akkor a módszer implicit és a következő yi + s érték az adódó nemlineáris egyenletből számítandó. Vegyük észre, az explicit módszernél minden lépésben egy új f ( x, y ) típusú függvény-kiértékelés kell, (mivel a többi már korábbról megvan), ugyanakkor az implicit módszernél még ügyes esetben is legalább 2-3 kiértékelésre szükség van. A mondottakat két egyszerű példán szemléltejük. Középpontszabály. Explicit, 2-lépéses módszer, a centrális differenciahányadosból származtatható:

yn +1 = yn −1 + 2hf n , n = 1, 2,…, N

(20.18)

tehát α 0 = −1, α1 = 0, α 2 = 1, β 0 = 0, β1 = −2, β 2 = 0 . Bebizonyítható, hogy a módszer másodrendben konzisztens és stabil, ha f a második változójában eleget tesz a Lipschitz-feltételnek. A középpontszabály indításához legalább másodrendűen pontos y1 értéket kell előállítani, amit megtehetünk valamely Runge-Kutta vagy Taylor-polinomos módszerrel, különben a másodrendű konvergencia nem lesz igaz. Trapézszabály. Ez implicit módszer. Úgy származtatható, hogy a differenciálegyenletet átírjuk integrál alakba és a jobb oldalon keletkező integrált a trapézmódszerrel közelítjük: yn +1 = yn +

h ( f n + f n+1 ) , n = 0,1,…, N − 1. 2

(20.19)

így s = 1, α 0 = −1, α1 = 1, β 0 = β1 = 1/ 2. Itt a következő pont előállításához meg kell oldanunk y -ra az y = yn +

h ( f n + f ( xn + h, y ) ) = F ( y ) 2

25 fixpont egyenletet. F ( y ) az általános feltételünk alapján eleget tesz a Lipschitz-feltételnek hL / 2 állandóval, így F ( y ) kontrakció, ha h < 2 / L . Célszerű formája az iterációnak: yn0+1 = yn + hf n , ynk++11 = yn +

(

kezdőérték az Euler-módszerből

)

h f n + f ( xn +1 , ynk+1 , k ≥ 0. 2

(20.20)

Leállás: ha ynk++11 − ynk+1 < ε (1 − q ) , ahol q a konvergencia-tényező és ε a kívánt hibakorlát. A trapézmódszer másodrendben konzisztens és stabil, ha f ∈ C 2 ([0,1] × ) és h < 1/ L , tehát ekkor másodrendű konvergenciára számíthatunk. További többlépéses módszereket a zárt vagy a jobbról nyilt kvadratúra-formulák segítségével lehet származtatni.

20.6. Aszimptotikus stabilitás A gyakorlati számítások szempontjából nem elegendő, ha egy módszer konzisztens és stabil. A kezdeti hiba ( - ha y0 is hibával terhelt, -) továbbterjedése szempontjából fontos egy további stabilitási tulajdonság. Egy differenciálegyenlet-megoldó numerikus módszer aszimptotikusan stabil, ha az

y (0) = 1, y ′( x) = q y ( x), q < 0

(20.21)

tesztfeladatra alkalmazva a kapott numerikus közelítések sorozatára minden h > 0 lépésköz esetén fennáll yn → 0 , ha n → ∞ . Ezt a stabilitási fogalmat a szakirodalom gyakran A0 -stabilitásnak nevezi. A tesztfeladat megoldása x növekedésével zérushoz tart: y ( x) = e qx → 0, x → ∞ , mivel q negatív. Így a definició azt követeli a módszertől, hogy a lépésköztől függetlenül a numerikus megoldás is tartsa meg ezt a lecsengő jelleget. Kimutatható, az Euler-módszer nem aszimptotikusan stabil, de a trapézmódszer az.

16. Nemlineáris egyenletek megoldása I

Recommend Documents