Numerikus módszerek - jegyzet. Kupán Pál

Numerikus módszerek - jegyzet

Kupán Pál

Tartalomjegyzék 1. fejezet. Számábrázolás, hibaszámítás.

5

1.1. Számítógépes számábrázolás

5

1.2. A hibaszámítás alapfogalmai

9

2. fejezet. Számsorok

17

2.1. Pozitív tagú sorok kiszámítása

18

2.2. Váltakozó el®jel¶ sorok

19

2.3. A sorok konvergenciájának a javítása

21

3. fejezet. Egyenletek numerikus megoldása

31

3.1. A gyökök elkülönítése

31

3.2. A felez® módszer

32

3.3. A húr módszer (regula falsi)

34

3.4. Az érint® (Newton) módszer

36

3.5. A szel® módszer

38

3.6. A Steensen-féle módszer

39

3.7. A fokozatos közelítések módszere

40

3.8. Vektor és mátrix normák

42

3.9. Algebrai egyenletek numerikus megoldása

46

4. fejezet. Lineáris egyenletrendszerek numerikus megoldása

53

4.1. Sajátos alakú egyenletrendszerek

53

4.2. Direkt módszerek

54

4.3. Lineáris egyenletrendszerek iteratív megoldása

63

4.4. Hibabecslés, kondícionálás

70

4.5. Túlhatározott lineáris egyenletrendszerek

71

5. fejezet. Sajátérték, sajátvektor numerikus kiszámítása 5.1. Krylov módszer

75 76

3

TARTALOMJEGYZÉK 5.2. Hatvány módszer 6. fejezet. Függvény approximáció, interpoláció

4 77 78

6.1. Analitikus függvények numerikus kiszámítása Taylor sorral

78

6.2. Interpoláció

79

6.3. Kétváltozós lineáris (bilineáris) interpoláció

97

6.4. Függvény approximáció a legkisebb négyzetek módszerével 100 7. fejezet. Bézier görbék, Bézier felületek

104

7.1. Bézier görbe szerkesztése "divide et impera" algoritmussal 104 7.2. Négyzetes Bézier görbék

106

7.3. Harmadfokú Bézier görbék

108

7.4. Bézier felületek

111

8. fejezet. Dierenciálegyenletek

115

8.1. Els®rend¶ dierenciálegyenletek

115

8.2. Els®rend¶, dierenciál egyenletrendszerek

121

8.3. Magasabb-rend¶ dierenciálegyenletek

124

Irodalomjegyzék

128

1. FEJEZET

Számábrázolás, hibaszámítás. 1.1. Számítógépes számábrázolás Tradicionálisan az emberi gondolkodás a számokat tízes (decimális)

1205 = 1000 + 200 + 0 + 5 = Általánosabban, ha r jelöli a számrendszer

számrendszerben kezeli-értelmezi. Például

3

2

1

0

1·10 +2·10 +0·10 +5·10 . alapját akkor egy a szám felírása

a következ®:

a = am · rm + am−1 · rm−1 + ..., ahol

ai ∈ {0, 1, 2, ..., r − 1}.

A számítógépek bináris, oktális, illetve hexadecimális számrendszert használnak. Számrendszer

Alap=r

Szimbólumok=ai

decimális

10 2 8 16

0, 1, 2, 3, 4, 5, 6, 7, 8, 9 0, 1 0, 1, 2, 3, 4, 5, 6, 7 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, A (= 10) , B (= 11) , C (= 12) , D, E, F

bináris oktális hexadecimális

10010110101(2) = 1 · 20 + 0 · 21 + 1 · 22 + 0 · 23 + 1 · 24 + 1 · 25 + 0 · 26 + 1 · 27 + 0 · 28 + 0 · 29 + 1 · 210 = 1205(10) 2265(8) = 5 · 80 + 6 · 81 + 2 · 82 + 2 · 83 = 1205(10) 4B5(16) = 5 · 160 + B · 161 + 4 · 162 = 5 · 160 + 11 · 161 + 4 · 162 = 1205(10) . 1. Példa.

Decimális számrendszerb®l egy osztással kapjuk meg (mod

r).

helyiérték¶ számjegyet adja, stb.

r számrendszerbe való áttérést maradékos

Az els® maradék a szám legkissebb Például decimális számrendszerb®l

binárisba való átszámításkor a keresett szám számjegyei a kett®vel való osztási maradékok, melyeket fordított sorrendbe olvasunk.

5

Pl.

1.1. SZÁMÍTÓGÉPES SZÁMÁBRÁZOLÁS

1205(10) = 10010110101(2)

mert

osztó

maradék

1205

1

602

0

301

1

150

0

75

1

37

1

18

0

9

1

4

0

2

0

1

1↑

.

6

Az osztást addig

végezük amíg a hányados nulla.

Hasonlóan

4B5(16)

mert

1205(10) = 2265(8)

osztó

maradék

1205

5

75

11=B

4

4

mert

osztó

maradék

1205

5

150

6

18

2

2

2

,

1205(10) =

↑

.

↑

Az oktális, illetve hexadecimális számrendszerek el®nye a tömörség, vagyis egy szám ábrázolásához kevesebb szimbolumra van szükség. Ugyanakkor, a kettes számrendszerb®l nagyon könny¶ az áttérés az oktális

(23 = 8),

illetve hexadecimális

(24 = 16)

számrendszerbe. Egy

bináris szám oktális alakját úgy kapjuk meg, hogy a bináris szám számjegyeit hármasával csoportosítjuk (jobbról), majd a megfelel® oktális

1.1. SZÁMÍTÓGÉPES SZÁMÁBRÁZOLÁS *(2)

*(8)

000

0

001

1

010

2

értéket 011

3

100

4

101

5

110

6

111

7

rendeljük hozzá. Például

7

10 |{z} 010 |{z} 110 |{z} 101

(2)

= 2265(8) .

A bináris-hexadecimális átalakítás hasonlóan történik, csak a bináris számjegyeket négyesével csoportosítjuk, például

4B5(16) .

100 |{z} 1011 0101 |{z}

(2)

=

Egy másik különbség az emberi, illetve a számítógép számábrázolása között a valós számok ábrázolásánál jelenik meg, nevezetesen az emberek xpontosan értelmezik a számokat, a számítógépben pedig lebeg®pontosan. A xpontos módszert a gépek csak egész számok ábrázolására használják. Ha egy

a szám nem fér el a xpontos ábrázolási tartományban akkor át

kell térni lebeg®pontos számábrázolásra. Egy decimális szám lebeg®pontos felírása a következ®képpen lehetséges:

a = (−1)s · 0.f · re

(1.1.1)

s =el®jel (sign), f =a szám törtrésze (mantissza) normalizált alakban, r =a számrendszer alapja, e =a szám kitev®je (exponent). 1 4 Például a −1234.5 = (−1) ·0.12345·10 szám egy 8 bites regiszteren a következ®képpen ábrázolható: 1 0 4 1 2 3 4 5 . Tehát a 8 bitb®l egy bit a szám el®jele (az els® 1 vagyis negatív), két bit a kitev® (04), ®t bit pedig a mantissza (12345). Mivel a kitev®t két biten ábrázoljuk az általa leírt tartomány 0 : 99. További bitre volna ahol

szükség a kitev® el®jelére. Ezt egyszer¶ eltolással megoldható, vagyis a kitev®höz

(04-hez)

hozzáadunk

50-et,

a kitev® tartomány felét. így az

el®bbi szám a következ®képpen ábrázolható:

1

5 4

1

2

3

4

5 .

Tehát a kitev®t úgy számítjuk ki, hogy a regiszterben szerepl® értékb®l kivonunk

50-t.

1.1. SZÁMÍTÓGÉPES SZÁMÁBRÁZOLÁS 2. Példa. A 0

5 5

5

55−50

4

3

10 = 54321 számot jelöli, míg 1048−50 = −0.0024680 számot. 3. Példa. A

357.02468

2 1

8

1 nyolc bites ábrázolás a

4 8

2

4

6

8

0 a

+0.54321·

−0.24680·

szám 8-bites ábrázolásához a következ®

lépéseket végezzük:

0.35702468 · 10−3 ; −3 jegy¶re vágjuk: 0.35702 · 10 ; 4 7 3 5 7 0 2 .

- normál alakba írjuk: - a mantisszát öt - ábrázoljuk:

0

A számítógépekben a valós számok ábrázolása hasonlóan történik a leírtakkal, csak nem tízes- hanem kettes számrendszerben, és nem 8 biten hanem (általában) 32 vagy 64. Kettes számrendszerben a szám lebeg®pontos ábrázolása a következ®:

a = (−1)s · (1 + f ) · 2e = (−1)s · 1.f · 2e .

(1.1.2) ahol az

s, f, e paraméterek jelent®sége ugyanaz mint a (1.1.1) képletben.

4. Példa.

3.55(10) = 11.10001100110011... = (−1)0 ·1.110001100110011...·

21 Ahhoz, hogy különböz® processzorokon végzett számítások eredményei ne térjenek el egymástól számottev®en, a lebeg®pontos ábrázolást 1985ben szabványosították (IEEE Standards 754). A 64 bites regiszter esetében (double precision) egy bit a szám el®jelét tárolja, 11 a kitev®t és 52 a mantisszát: s

e e e e e e e e e e e f f f f f f f f f f f f ... 0 1 10 Mivel a kitev® tartománya 0 és 2 + 2 + ...2 = 211 − 1 között van,

a kitev® el®jele szabályázoható ha eltoljuk a tartomány felével vagyis

1023-al.

A binárisan ábrázolt mantisszát mindig

1-re

normalizáljuk.

Innen következik, hogy még egy bitet spórolhatunk meg ha ezt az

1-t

nem ábrázoljuk (implicitbites ábrázolás).

11.34(10) = 1011.010101110...(2) = 1.011010101110... · 23 szám el®jele + (s = 0) , a kitev®je 3+1023 = 1026(10) = 10000000010(2) , és a mantisszája 011010101110..., tehát a szám ábrázolása: 0 1 0 0 0 0 0 0 0 0 1 0 0 1 1 0 1 0 1 0 1 1 1 0 5. Példa. A

...

1.2. A HIBASZÁMÍTÁS ALAPFOGALMAI

9

6. Példa. Az

1 0 0 0 0 0 0 0 1 0 1 1 1 1 0 1 1 0 1 1 1 0 ábrázolás esetén az el®jel negatív (s = 1), a kitev® e = 10000000101(2) = 210 + 22 + 20 = 1029(10) , a mantissza 1.f = 1.111011011100... tehát a 1029−1023 szám =−1.111011011100... · 2 = −123.45 számot jelöli. 1

Gépi epszilonnak nevezzük az

0

...

a = 1 és a rákövetkez® szám különbségét.

A két szám ábrázolása: 0

0 1 1 1 1 1 1 1 1 1 1

0

0

0

0

0

0

...

0

0

0

0

0

0

0

0

0

0

...

0

0

0

1

és a rögtön rákövetkez® szám: 0

0 1 1 1 1 1 1 1 1 1 1

ahonnan

eps = (1 + 0.0 . . . 01)·22

10 −1−1023

10 −1−1023

−(1 + 0)·22

A kitev® (eltolt) két széls®értéke:

= 0.0...01·20 = 1·2−52 ' 2.2·10−16 (10)

−1023 és 1024, különös jelent®séggel

bírnak, ezért ezeket csak sajátos esetben használjuk. Tehát a számok ábrázolásához a kitev® használható tartománya:

−1022 ≤ e ≤ 1023. Innen következik, hogy a legkisebb (pozitív) ábrázolható számot (jel.

e = −1022 és f = 0 értékekre kapjuk, vagyis realmin=2−1022 (2) ≈ 2.2251 · 10−308 (10) . A legnagyobb számot e = +1023 és f = 0.11 . . . 1 1023 értékekre kapjuk, vagyis: realmax=(2 − eps) 2 ≈ 1.7977 · 10308 (10) . realmin) az

Túlcsordulásról beszélünk ha a számítások során egy eredmény nagyobb realmax-nál, illetve alulcsordulásról ha az eredmény kisebb realmin-nél.

1.2. A hibaszámítás alapfogalmai A gyakorlatban, a m¶veleteket általában a számok közelít® értékével végezzük, és az eredmény is közelít® érték¶ lesz. Lényeges követelmény, hogy mindig becslést tudjunk adni a közelítés pontosságáról. Hibaforrás lehet már a bemen® adat; ez bekövetkezhet a híbás mérés, vagy a pontatlan m¶szerek miatt.


1.2.1. A közelít® érték és hibája A gyakorlatban

Tekintsük az

10

A pontos értéket.

A helyett ennek egy közelít® érték ével dolgozunk, jelöljük

a-val. 7. Példa. A

π közelít® értékének az a = 3.14 értéket lehet tekinteni.

A pontos es a közelít® érték abszolút különbségét az

hibáj ának nevezzük: |A − a| = ∆.

a szám abszolút

Mivel a pontos érték legtöbbsz®r ismeretlen ezért az abszolút hiba is ismeretlen. Helyette egy

α

legkisebb

hibakorlátot

fogunk használni:

∆ = |A − a| ≤ α. A fenti képlet felírható a következ®képpen:

a − α ≤ A ≤ a + α. 8. Példa. Ha egy tetsz®leges mérésnél egy bizonyos szavatolunk úgy

0.1

gramm a hibakorlát

0.1

grammért

(A = a ± 0.1).

Az abszolút hiba es a hibakorlát a mérés és számolás pontosságának a jellemzésére nem elegend®. Például ha mérések alapján

L1 = 145cm±

0.1cm, L2 = 6.2cm±0.1cm, habár a két hibakorlát azonos, az els® mérés min®sége jobb, mint a másodiké. Az

A

érték

relatív (viszonylagos) hiba ∇=

értelmezése a kovetkez®:

∆ |A − a| = . |a| |a|

A gyakorlatban, hasonlóan az abszolút hibához, nem a relatív hibát vizsgáljuk hanem annak a legkisebb korlátját.

hibakorlátnak

és jelöljük

Ezt nevezzük

relatív

δ -val: ∆ = ∇ ≤ δ. |a|

Ismerve az abszolút hibakorlátot mivel

δ=

∆≤α

α . |a|

9. Példa. Az el®bbi példát felhasználva

−4

10 , δ2 =

0.1 6.2

−2

= 1. 612 9 × 10 .

felírhatjuk:

δ1 =

0.1 145

= 6. 896 6 ×

1.2. A HIBASZÁMÍTÁS ALAPFOGALMAI A továbbiakban az abszolút, illetve relatív hibát

11

∆,

illetve

δ -vel

azonosítjuk.

1.2.2. A helyes jegyek száma

A számok relatív hibája könnyen

kiszámítható ha ismert e számok tizedes alakja illetve a helyes jegyek száma.

Els® értékes számjegy nek

nevezzük az els® nem nulla számjegyet

amit balról jobbra haladva találunk. 10. Példa. A értékes számjegye

3.14 = 3 × 100 + 1 × 10−1 + 4 × 10−2 számnak az els® −3 a 3, a 0.0027 = 2 × 10 + 7 × 10−4 -nek pedig a 2. a = am × 10m + am−1 × 10m−1 + ... számnak az els® am . Innen következik, hogy:

Általánosan egy értékes számjegye

am 10m ≤ a ≤ am 10m + 10m = (1 + am )10m . A tizedestört alakjaban felírt

a szám valamely meghatározott jegyének

helyén álló mindegyik egységnek megvan a maga

m−1

10 -el egyenl®, a második 10 m−n+1 egyenl® 10 .

helyen álló egység álló egység

m

helyiértéke : ..., az

az els®

n-edik helyen

3.141 számban a harmadik jegy helyiértéke 10−2 , míg 0.3141 számban ugyanannak a jegynek 10−3 helyi értéke van. 11. Példa. A

a

Hogy egyszer¶bbé teggyünk egy tizedes számot használhatjuk a

csonkitási szabaly t, azaz, elhanyagolhatunk bizonyos számu számjegyet a végér®l. Az abszolút hiba, amelyet ebben az esetben elkövetünk, nem fogja meghaladni a meghagyott jegyek utolsójának helyi értékét.

2.10527 számnak az els® három jegyére szorítkozunk, −2 akkor a 2.10 számot kapjuk, az elkövetett hiba pedig kisebb mint 10 . 12. Példa. Ha a

Ha nagyobb pontosságot szeretnénk elérni használhatjuk a

szabály t.

kerekítési

Ez abban áll, hogy a meghagyott számjegyek utolsóját változatlanul

hagyjuk vagy növeljük egy egységgel attól függ®en, hogy az els® elhanyagolt számjegy kisebb (<), vagy nagyobb vagy egyenl® (≥)

5-el.

Ebben az

esetben, az elkövetett abszolút hiba nem haladja meg a meghagyott jegyek utolsójának helyi értékenek a felét.


12

2.10527 számban, a kerekítesi szabályt alkalmazva, jegyet tartjuk meg, akkor kapjuk a 2.11 számot, az −2 pedig kisebb mint 1/2 × 10 .

13. Példa. Ha a az els® három elkövetett hiba

Altalánosan, ha az

A számot az a érték n helyes számjeggyel közelíti

meg, akkor az abszolút hibára a következ® egyenl®tlenségeket kapjuk:

∆ = |A − a| ≤ 10m−n+1 csonkítási szabállyal ill. 1 ∆ = |A − a| ≤ 10m−n+1 kerekítési szabállyal. 2 Ezen feltételek mellett kifejezhet® természetesen a relatív hiba is:

10m−n+1 10m−n+1 1 ≤ =⇒ δ ≤ m |a| am 10 am 10n−1 1 kerekítés esetben. δ ≤ 2am 10n−1

δ =

Nagy vagy igen kicsi számokat ajánlatos

10-nek

csonkítás ill.

a hatványai szerint

írni.

12345678 számot három helyes számjeggyel a következ®képpen írjuk: 123 × 10 . Hasonlóan a 0.000001234 számot három helyes jeggyel −7 =12.3 × 10 . 14. Példa. A

5

A problemát fordítva is meg lehet fogalmazni, vagyis egy adott pontosságra hány helyes számjegyet kell számításba venni.

1.2.3. Hibaterjedés

1.2.3.1. Az alapm¶veletek hibái közelít® értékei.

a1 , a2 az A1 , A2 pontos számok eltérés hibakorlátjait: α1 , α2

Legyen

Ismerjük továbbá az

illetve δ1 , δ2 . Szeretnénk meghatározni a hibakorlátokat abban az esetben mikor az

•

A1 , A2

értékekkel a négy alapm¶veletet végezzük.

Az összeg hibakorlátja

Ebben az estben az

a = a1 + a2

összeg

α

A = A1 + A2

mennyiség megközelítéséül szolgaló

abszolút hibakorlátot határozzuk meg.

|A − a| = |A1 − a1 + A2 − a2 | ≤ |A1 − a1 | + |A2 − a2 | ≤ α1 + α2

1.2. A HIBASZÁMÍTÁS ALAPFOGALMAI vagyis az

α

hibakorlátot egyenl®nek lehet venni az

α1 , α2

13 hibakorlátok

összegével:

α = α1 + α2 .

(1.2.1)

•

A különbség hibakorlátja

Hasonlóan lehet eljárni mint az el®bbi pontnál:

A = A1 − A2

,

a =

a1 − a2 =⇒ |A − a| = |A1 − a1 − A2 + a2 | ≤ |A1 − a1 | + |A2 − a2 | ≤ α1 + α2 vagyis, mint az el®bb, az venni az

α1 , α2

α

különbség hibakorlátját egyenl®nek lehet

hibakorlátok ®sszegével:

α = α1 + α2 .

(1.2.2)

Külön gyelmet kell szentelni ha

A1 es A2 közeli értékek ugyanis, ebben

az esetben sokat veszíthetünk a pontosságból.

15. Példa. Számítsuk ki két kocka

2.140 m. A térfogatok három pontos jeggyel a következ®k lesznek: V1 = 9.81, V2 = 9.80, ahonnan a különbségük V1 − V2 = 0.01 = 10−2 . Ugyanakkor felhasználva a Taylor képletet: ha oldalaik

2.141 m

V1 , V2 térfogatainak különbségét

ill.

f (x + ∆x) − f (x) = f 0 (x)∆x f (x) = x3 , x = 2.140, ∆x = 0.001 azt kapjuk, hogy: V1 − V2 = 3 ∗ 2.1402 ∗ 0.001 = 1.37 × 10−2 ami pontosabb eredmény az el®bbinél. ahol

•

A szorzat hibakorlátja

|A1 A2 − a1 a2 | = |A1 A2 − A1 a2 + A1 a2 − a1 a2 | = |A1 (A2 − a2 ) + (A1 − a1 ) a2 | ≤ ≤ |A1 | |A2 − a2 | + |a2 | |A1 − a1 | ≤ |A1 | ∆2 + |a2 | ∆1 . Az ismeretlen

|A1 |

tagot majoráljuk a következ® taggal:

|A1 | = |A1 − a1 + a1 | ≤ |A1 − a1 | + |a1 | ≤ ∆1 + |a1 |


14

∆ ≤ (∆1 + |a1 |) ∆2 + |a2 | ∆1 = ∆1 ∆2 + |a1 | ∆2 + |a2 | ∆1 . Elosztva egyenl®tlenség mindkét oldalát |a1 a2 | taggal kapjuk, hogy:

tehát az

∆ ∆1 ∆2 ∆1 ∆2 ≤ + + |a1 a2 | |a1 a2 | |a1 | |a2 | vagyis:

δ = δ1 δ2 + δ1 + δ2 . Mivelhogy a relatív hiba kicsi szokott lenni ezért a δ1 δ2 szorzat elhanyagolható és a szorzat relatív hibájának vehetjük a

A1 A2 −

relatív hibák ®sszegét:

δ = δ1 + δ2 .

(1.2.3)

•

δ1 , δ2

A hányados hibakorlátja

a1 2 a1 = A1 aA2 −A a2 2 a2 |A1 −a1 ||a2 |+|a1 ||a2 −A2 | |A2 ||a2 |

=

|A1 a2 −a1 a2 +a1 a2 −A2 a1 | |A2 a2 |

=

|(A1 −a1 )a2 +a1 (a2 −A2 )| |A2 a2 |

≤

|A2 | ismeretlent minoráljuk a következ® képpen: |a2 | = |a2 − A 2 + A2 | ≤ |a2 − A2 |+|A2 | =⇒ |A2 | ≥ |a2 |−|a2 − A2 | ahonnan A1 a1 ∆1 |a2 |+|a1 |∆2 ∆ = A − a2 ≤ |a2 |(|a2 |−∆2 ) . 2 a Hogy a relatív hibakorlátot megkapjuk, elosztjuk ∆-t 1 -vel es azt a2

Az

kapjuk, hogy:

∆1 + |a∆22| ∆ |a | + |a1 | ∆2 |a2 | ∆ |a1 | ≤ 1 2 = . a1 |a2 | (|a2 | − ∆2 ) |a1 | 1 − |a∆22| a2 Figyelembe véve, hogy

∆2 elhanyagolható |a2 |

1-hez

képest azt kapjuk,

hogy: (1.2.4)

δ = δ1 + δ2 .

1.2.3.2. Függvény hibakorlátja

Ebben az esetben szeretnénk megvizsgálni

hogyan befolyasolja a véltozó hibakorlátja a függvény hibakorlátját.

f : Rn −→ R egy deriválható függvény, illetve X 0 = (X10 , X20 , ..., Xn0 ) ∈ Rn egy pontos érték. Kérdés mennyi f (X 0 ) =? 0 0 Az X ideális értéket megközelítjük egy x = (x01 , x02 , ..., x0n ) ∈ Rn ponttal, illetve az f (X 0 ) értéket az f (x0 ) értékkel. A függvény 0 0 0 megközelítésében elkövetett abszolút hibát jelöljük ∆f (x ) = |f (X ) − f (x )|. Legyen


15

A középérték tételb®l:

f X 0 − f x0 = df (ξ) X 0 − x0 , ξ ∈ X 0 , x0 következik, hogy:

∆f x

0

n n X X ∂f ∂f 0 0 = (ξ) Xi − xi ≤ (ξ) Xi0 − x0i . ∂xi ∂xi i=1

i=1

Tételezzük továbbá fel, hogy amelyre

x0 ∈ Ω

es

f

X 0 pontnak létezik egy olyan Ω környezete

parciális deriváltjai korlátosak, vagyis:

∂f sup (ξ) = Mi ∈ R. x∈Ω ∂xi Következik, hogy:

n X f X 0 − f x0 = ∆f x0 ≤ Mi Xi0 − x0i . i=1 Ha

egy el®re megadott pontosság akkor kikötjuk, hogy:

0 Xi − x0i ≤

, ∀i nMi

ahonnan az el®bbi egyenl®tlenség alapján kapjuk:

f X 0 − f x0 ≤ . f (x) = x1 + x2 , akkor Mi = 1 és visszakapjuk az összegre 2 vonatkozó hibakorlátot, x = (x1 , x2 ) ∈ R . Hasonlóan, visszakapjuk a hibabecslést a kivonásra, szorzásra, osztásra ha f (x) rendre f (x) = x1 − x2 , f (x) = x1 x2 , f (x) = xx12 . Ha

16. Példa. Legyen

f : R2 −→ R, f (x) = x21 + 2x22 + x1 x2 ,

ahol

2

x = (x1 , x2 ) ∈ R . ∆f x0 = M1 ∆x1 + M2 ∆x2 ∂f ∂f = sup |2x + x | , M = sup ahol M1 = sup ∂x2 = sup |4x2 + x1 | . 1 2 2 ∂x1 0 0 Konkrétan ha x = (0, 0) , X = (0.1, −0.025), Ω = [−0.1, 0.1] × [−0.025, 0.025], akkor M1 = supx∈Ω |2x1 + x2 | = 0.225, M2 = supx∈Ω |4x2 + x1 | =


0.2

ahonnan

∆f x0 = 0.225 · 0.1 + 0.2 · 0.025 = 0.027 5.

16

2. FEJEZET

Számsorok 17. Definíció. Számsornak nevezzük az alábbi végtelen összeget:

∞ X

(2.0.5)

an = a1 + a2 + ... + an + ...

n=1 A sor természetét -konvergencia, divergencia- az

(Sn )n

részletösszeg

sorozat segítségével tanulmányozzuk

Sn = a1 + a2 + ... + an .

(2.0.6)

Ha az (Sn )n sorozatnak van véges S = limn→∞ Sn határértéke, akkor a P∞ n=1 an sort konvergensnek nevezzük és a sor összege egyenl® S -el; ha viszont ez a határérték végtelen, vagy a sorozatnak nincs határértéke, akkor a sort divergensnek nevezzük.

Könnyen belátható hogy ha az

(an )n

sorozat nem konvergál zéróhoz

akkor a sor divergens, tehát szükséges (de nem elégséges) feltétele P∞ annak hogy a n=1 an sor konvergens legyen az hogy:

lim an = 0.

(2.0.7)

n→∞

18. Példa. A

sor általános

p p 0, 01 + 2 0, 01 + ... + n 0, 01 + ... √ n tagja an = 0, 01 → 1 6= 0, tehát a sor

divergens.

19. Példa. Az úgynevezett harmonikus sor esetén:

∞ X 1 1 1 = 1 + + + ..., n 2 3 n=1

17

2.1. POZITÍV TAGÚ SOROK KISZÁMÍTÁSA annak ellenére hogy

an → 0, a sor divergens.

18

Az általánosított harmonikus

sor:

∞ X 1 nα n=1

(2.0.8)

konvergens ha

α>1

és divergens ha

α ≤ 1.

A konvergens sorok összegének a megközelítéséhez a sor els®

n

tag

összegét használjuk:

∞ X

(2.0.9)

an ' a1 + a2 + ... + an ,

n=1 és ekkor az elkövetett hiba a következ® lesz:

Rn = |an+1 + an+2 + ...|

(2.0.10)

Természetesen, minél több tagot veszünk gyelembe a (2.0.9) közelítéshez, annál pontosabb eredményt kapunk, de az összegezend® tagok számát úgy fogjuk meghatározni hogy az abszolút hiba egy bizonyos

> 0

ponotossági küszöbön belül maradjon

Rn < .

(2.0.11)

2.1. Pozitív tagú sorok kiszámítása Tekintsük a

P∞

n=1

an , an > 0,

pozitív tagú sort amelyre teljesül az

alábbi D'Alembert-féle konvergencia kritérium: (2.1.1)

an+1 ≤ q < 1, ∀n ≥ N. an

Tehát:

aN +1 ≤ qaN , aN +2 ≤ qaN +1 ≤ q 2 aN , ... és az eredeti sort majorálhatjuk a következ®képpen:

∞ X

an = a1 + ... + aN + aN +1 + ... ≤ a1 + ... + aN + qaN + q 2 aN + ...

n=1

= a1 + ...aN −1 + aN (1 + q + q 2 + ...) = a1 + ...aN −1 + aN lim( n

= a1 + ...aN −1 + aN (

1 ), 1−q

mert

q ∈ (0, 1) .

1 − qn )= 1−q

2.2. VÁLTAKOZÓ ELJEL SOROK

19

Tehát, a sor összegének a közelítésére használjuk a

∞ X

(2.1.2)

an ' a1 + a2 + ... + aN −1 ,

n=1 és a hiba

RN = aN (

(2.1.3) Az

N

1 ). 1−q

küszöbszám kiszámításához megoldjuk a (2.0.11) kikötésb®l ered®

egyenl®tlenséget:

aN (

(2.1.4)

1 ) ≤ . 1−q

Mivel numerikusan az (2.1.2) összeget egy DO-WHILE ciklusban számítjuk ki, a (2.1.4) egyenl®tlenség konkrét megoldására nincs szükség, ugyanis az index számot addig növeljük (és összegezünk egyúttal) ameddig az (2.1.4) egyenl®tlenség igaz lesz. 20. Példa. A

an+1 = an

P∞

n=1

1 2n+1 (n+1)2 1 2n n 2

=

an =

P∞

1 n=1 2n n2 sor esetében:

1 2n n2 1 ≤ = q < 1. = 2 2 2 2n+1 (n + 1) 2 1 + n1

2 −2 akkor a (2.1.4) szerint N 2 ≤ 10 =⇒ 2 N P∞ 1 1 1 1 és a sor közelít® értéke (2.1.2) n=1 2n n2 ' 2 + 22 22 + 23 32 =

Tehát, ha például

= 10−2

N = 4, 0.576 39.

2.2. Váltakozó el®jel¶ sorok A váltakozó el®jel¶ sorok az alábbi alakban irhatók:

(2.2.1)

∞ X

(−1)n+1 an = a1 − a2 + a3 − ... , an > 0

n=1 és esetükben alkalmazhatjuk a Leibniz-féle kritériumot, azaz ha

(an )n

monoton csökken® és zéróhoz konvergál, akkor a (2.2.1) sor konvergens. A sor összegének a közelítéséhez újból az els®

(2.2.2)

∞ X n=1

N

tag összegét vesszük

(−1)n+1 an ' a1 − a2 + a3 − ... + (−1)N +1 aN ,

2.2. VÁLTAKOZÓ ELJEL SOROK

20

és a hiba:

RN = (−1)N +2 aN +1 + (−1)N +3 aN +2 + ... .

(2.2.3)

Feltételezve hogy teljesül a Leibniz-féle kritérium, a (2.2.3) hiba majorálható az els® elhanyagolt taggal (moduluszban)

RN ≤ aN +1 .

(2.2.4) Valóban, ha

N

páros azt kapjuk hogy

RN = | aN +1 − aN +2 + aN +3 − aN +4 + ...| = = aN +1 − aN +2 + aN +3 − aN +4 + ... ≤ aN +1 , ha pedig

N

páratlan

RN = | −aN +1 + aN +2 − aN +3 + aN +4 + ...| = = aN +1 − aN +2 + aN +3 − aN +4 + ... ≤ aN +1 . A fenti képletekben felhasználtuk hogy

(an )n monoton csökken®: aN +1 ≥ aN +2 ≥

aN +3 ≥ .... Az összegezend® tagok számát

N -et,

újból a

> 0

pontosságtól

tesszük függ®vé, vagyis kikötjük hogy

RN ≤ aN +1 ≤ ,

(2.2.5)

bebiztosítva ezáltal hogy az összegzésb®l ered® hiba nem lépi túl a megengedett határt.

21. Példa. Tekintsük a Leibniz sort

∞

X 1 1 1 1 − + − ... = (−1)n+1 . 2 3 n n=1 1 sorozat telejesti a Leibniz feltételeket, tehát a sor konvergens. n −2 A sor összegének = 10 pontossággal való kiszámításához felhasználjuk

Az

an =

a (2.2.5) kikötést

1 ≤ , N +1

2.3. A SOROK KONVERGENCIÁJÁNAK A JAVÍTÁSA ahonnan

N = 99.

∞ X

21

Tehát

(−1)n+1

n=1

1 1 1 1 ' S99 = 1 − + − ... + = 0, 6982 n 2 3 99

= 10−2 pontossággal. Ismerve, hogy a sor pontos összege=ln 2, kiszámítható −3 a hiba: RN = |ln 2 − 0, 6982| ' 5 ∗ 10 < . A fenti összegzést egy FOR ciklusban oldjuk meg; ugyanakkor, ha egy DO WHILE ciklust használunk akkor a (2.2.5) egyenl®tlenséget nem kell el®zetesen megoldani mert az összegzést addig folytatjuk míg ez teljesül (tehát egyesével lépegetve csak leellen®rizzük).

P∞ 22. Definíció. Egy n=1 P∞ ha n=1 |an | sor konvergens.

an sort abszolút konvergensnek nevezünk

A pozitív tagú sorok esetében az abszolút konvergenciának semmi jelent®sége nincs, viszont a váltakozó sorok esetében ha egy sor abszolút konvergens akkor az konvergens is lesz. A fordított állítás nem igaz. P∞ P∞ n+1 1 n+1 1 Például a Leibniz-féle sor konvergens, de n=1 (−1) n=1 (−1) n n P∞ 1 n=1 n divergens.

2.3. A sorok konvergenciájának a javítása A sorok összegének a kiszámításánal a pontos közelítés mellett fontos a gyorsaság is.

A fenti példán láttuk hogy a Leibniz sor

= 10−2

N = 99 tagot kellett összegezni. A pontosság −9 −9 növelésével a N is arányosan n®. Például ha = 10 akkor N = 10 − 1. Ahhoz hogy ugyanazt az eredményt kevessebb tag összegezésével pontosságú kiszámításához

érjük el szükséges az eredeti sor átalakítása. 23. Definíció. Azt mondjuk hogy a P∞ mint a n=1 an sor ha

P∞

n=1 bn sor gyorsabban konvergál

bn = 0. n→∞ an lim

(2.3.1)

2.3.1. Az Euler transzformáció S=

∞ X n=1

A váltakozó el®jel¶ sorok esetében

(−1)n+1 an = a1 − a2 + a3 − ..., ai > 0

=

2.3. A SOROK KONVERGENCIÁJÁNAK A JAVÍTÁSA az

Sn

22

részletösszeg helyett használjuk ezek középértékét:

Tn =

(2.3.2)

Sn−1 + Sn 2

és így a következ® sorhoz jutunk:

T =

(2.3.3)

a1 a2 − a1 a3 − a2 − + + ... 2 2 2

E sor összege megegyezik az eredeti sor összegével:

Tn − Sn = (−1)n

an → 0, 2

ha

n → ∞,

tehát

T = S,

viszont gyorsabban konvergál:

lim

an −an−1 2 (−1)n+1 an

(−1)n+1

n

1 an−1 = lim 1 − = 0. 2 n an

24. Példa. A Leibniz sorra alkalmazzuk az (2.3.3) átalakítást

∞ X

(−1)

n=1

n+1

1 1 = − n 2 =

1 2

−1 + 2

1 3

− 2

1 2

1 n+1 n

+ ... + (−1)

1 − n−1 + ... 2

1 1 1 1 + − + ... + (−1)n + ... 2 4 12 2n (n − 1)

Figyelembe véve a (2.2.5) hiba majorálási kritériumot:

1 ≤ 2N (N + 1) azt kapjuk hogy a sor

= 10−2 pontossággal való kiszámításához N = 7

tagra van szükség és

∞ X

(−1)n+1

n=1

1 1 1 1 1 = ln 2 ' + − + ... − . n 2 4 12 84

Az átalakítás többször is alkalmazható.

2.3.2. A Kummer transzformáció pozitív tagú sorból, válusszunk egy

P∞

n=1

ismert és létezik a (2.3.4)

lim n

an = λ ∈ R, un

un

sort

P∞

an , a n > 0 melynek u összege

Kiindulva

n=1

2.3. A SOROK KONVERGENCIÁJÁNAK A JAVÍTÁSA

23

határérték. Ekkor

∞ X

(2.3.5)

an = λ

n=1 Tehát, a

P∞

∞ X

un +

n=1

n=1

an

∞ X

(an − λun ) = λu +

n=1

∞ X

(an − λun ) .

n=1

sor tanulmányozása visszavezethet® az alábbi sor

elemzésére:

∞ X

(2.3.6)

(an − λun ) ,

n=1 ami viszont gyorsabban konvergál az eredeti sornál:

lim n

an − λun un = 1 − λ lim = 0. n an an

25. Példa. Alkalmazzuk az átalakítást az alábbi konvergens sorra

∞ X

an =

n=1 Segédsornak a

1-el

P∞

n=1

un =

∞ X 1 . 2 n n=1

P∞

1 n=1 n(n+1) sort vesszük aminek az összege

egyenl® és

λ = lim n

n (n + 1) an = lim = 1. n un n2

Tehát, a (2.3.5) képlet szerint

∞ ∞ X X 1 1 1 an = 1 + − =1+ 2 2 n n (n + 1) n (n + 1) n=1 n=1 n=1

∞ X

ami nyilván gyorsabban konvergál. A gyorsító módszerek többször is alkalmazható egymásután.A módszer többszer is alkalmazható egymásután.

2.3.3. Függvénysorok, hatványsorok

Az olyan sorokat, amelyeknek

tagjai egy (vagy több) változó függvényei, függvénysoroknak nevezzük. Jelölésük:

∞ X

(2.3.7)

fn (x) .

n=1

1 + x + x2 + ... + xn + ... = P sin nx + ... + sinn2nx + ... = ∞ n=1 n2 .

26. Példa. a) b)

sin x +

sin 2x 22

P∞

n=0

xn


24

A legfontossabb függvénysorok a hatványsorok:

∞ X

(2.3.8)

an (x − x0 )n ,

n=1 vagy

x := x − x0

transzformációt használva

∞ X

(2.3.9)

an x n .

n=1 A hatványsorok el®nye hogy a részletösszeg sorozat egy polinom függvény: P Sn = nk=1 ak xk .

(−R, R) intervallumot amelyre a (2.3.9) hatványsor abszolút konvergens konvergencia intervallumnak nevezzük, az RRet pedig konvergencia sugárnak. Ha R = 0 akkor a hatványsor csak x = 0-ban konvergens, ha pedig R = ∞ akkor az egész R tengelyen. Ha x eleme a konvergencia intervallumnak akkor a függvénysor egy f Azt a legnagyobb

függvényt fog el®állítani

f (x) =

∞ X

an x n ,

n=1

> 0 számhoz létezik egy N (, x) ∈ N küszöbszám (ami x-t®l) amelyre ha n > N , akkor

vagyis minden függ

és

|f (x) − Sn (x)| < .

(2.3.10) Ha az

N

küszöbszám nem függ

x-t®l

akkor a sorrol azt mondjuk hogy

egyenletesen konvergál.

Ha 0 < R < ∞ akkor a (2.3.9) hatványsor egyenletesen konvergál bármilyen [−r, r] intervallumon ahol 0 < r < R. 27. Tétel.

Az

R

konvergencia sugárnak a kiszámításához alkalmazhatjuk a

D'Alembert (hányados), illetve Cauchy (gyök) kritériumot: ha (2.3.11)

an+1 , l = lim n an

vagy (2.3.12)

l = lim n

p n

|an |,

2.3. A SOROK KONVERGENCIÁJÁNAK A JAVÍTÁSA határérték létezik, akkor

Ha

f

28. Tétel. P n (x) = ∞ n=1 an x ∞

25

R = 1l .

an xn hatványsor konvergencia sugara R, akkor , f : (−R, R) → R összegfüggvény P∞

n=1

(i) C (−R, R) osztályú és

an =

(2.3.13)

f (n) (0) , ∀n ≥ 0, n!

(ii) bármilyen [a, b] ⊂ (−R, R) intervallumon ˆ

b

f (x) dx =

(2.3.14)

a

∞ X n=1

ˆ

b

xn dx,

an a

vagyis a sor tagonként integrálható. 29. Példa. a)

1 − x + x2 − ... = lim n

1 − (−x)n 1 = , ∀x (−1, 1) . 1 − (−x) 1+x

Tagonként integrálva azt kapjuk hogy

1 1 c + x − x2 + x3 − ... = ln (x + 1) 2 3 x = 0 -ra következik c = 0. Tehát x − 21 x2 + 31 x3 − ... = ln (x + 1) x → 1-re 1 − 21 + 31 − ... = ln 2. és

és

b)

1 , ∀x (−1, 1) . 1 + x2 1 3 1 5 Tagonként integrálva kapjuk hogy x − x + x − ... = arctgx, és 3 5 x → 1 =⇒ 1− 13 + 15 −... = π4 , ahonnan π értéke kiszámítható bármilyen 1 − x2 + x4 − ... =

pontossággal.

2.3.4. Taylor sor

Az el®bbi fejezetben láttuk miként rendelhet®

hozzá egy konvergens sorhoz egy függvény. A kérdés fordítva is felmerül, vagyis ha adott egy

f

függvény el®állítható-

e hatványsorral?

I ⊂ R egy nyílt intervallum, x0 I és f : I→ R egy végtelenszer ∞ deriválható függvény: f C (I) . Legyen

2.3. A SOROK KONVERGENCIÁJÁNAK A JAVÍTÁSA Minden

f

26

függvényhez hozzárendelhetünk az úgy nevezett Taylor

sorát:

1 1 1 f (x) ≈ f (x0 )+ f 0 (x0 ) (x − x0 )+ f 00 (x0 ) (x − x0 )2 +...+ f (n) (x0 ) (x − x0 )n +... 1! 2! n! de ez nem jelent automatikus egyenl®séget a függvény és a Taylor sora között.

30. Definíció. Egy

x0 ∈ I -re, f

f

függvényt analitikusnak nevezzünk ha bármely

el®állítható hatványsorként

x0

környezetében (lokálisan):

∞ X 1 (n) f (x) = f (x0 ) (x − x0 )n , x ∈ Vx0 . n! n=0

(2.3.15)

Minden

R-en értelmezett elemi függvény analitikus, de nem minden

végtelenszer deriválható függvény analitikus.

31. Példa. Az

f : R → R, ( f (x) =

függvény végtelenszer deriválható: de nem analitikus. értékre

0

1 1

ex

0,

, x>0 x≤0

f ∈ C ∞ (R)

Különben (2.3.13) szerint

f (n) (0) = 0, n ≥ 0 f (x) = 0, minden x

és

környezetében ami hamis.

Ha f : I→ R egy végtelenszer deriválható függvény az I ⊂ R nyílt intervallumon: f C ∞ (I), és ∃M > 0 úgy, hogy ∀x ∈ 32. Tétel.

I, ∀n ≥ 0, (n) f (x) ≤ M,

(2.3.16)

akkor az f függvény analitikus. Numerikusan az nevezetesen az

n-ed

f

függvényt véges tagú összegként számítjuk ki,

fokú

Tn

Taylor polinom segítségével:

(2.3.17)

1 1 1 Tn (x) = f (x0 )+ f 0 (x0 ) (x − x0 )+ f 00 (x0 ) (x − x0 )2 +...+ f (n) (x0 ) (x − x0 )n . 1! 2! n!

2.3. A SOROK KONVERGENCIÁJÁNAK A JAVÍTÁSA 33. Tétel.

akkor:

27

Ha f ∈ C (n+1) (I) , x0 ∈ I és M = supx∈I f (n+1) (x) ,

f (x) = Tn (x) + Rn (x)

(2.3.18)

ahol az Rn hibatagra igaz az alábbi egyenl®tlenség: |R n (x)| ≤ M

(2.3.19)

1 |x − x0 |n+1 . (n + 1)!

Ha f C n+1 (I) és x0 I, akkor az Rn maradéktag felírható az alábbi Lagrange-féle alakban: (2.3.20)

Ha

(x − x0 )n+1 (n+1) f (ξ) , ξ (x0 , x) (vagy ξ (x, x0 ) ). R n (x) = (n + 1)!

x0 = 0,

Tn (x) = f (0) +

(2.3.21) Ha a

akkor a Taylor sorfejtés MacLaurin nevet viseli:

h = x − x0

x2 xn x 0 f (0) + f 00 (0) + ... + f (n) (0) . 1! 2! n!

jelölést használjuk, akkor a (2.3.18) Taylor képlet

az alábbi alakban adható meg: (2.3.22)

h2 hn hn+1 (n+1) h f (x0 + θh) f (x0 + h) = f (x0 )+ f 0 (x0 )+ f 00 (x0 )+...+ f (n) (x0 )+ 1! 2! n! (n + 1)! ahol

θ (0, 1) .

34. Példa. a) Az függvényekre

M =1

f, g : R → R, f (x) = sin (x) , g (x) = cos (x)

és a (2.3.21) MacLaurin sorfejtésük

∞

X x2n+1 x x3 − + ... = (−1)n , sin (x) = 1! 3! (2n + 1)! n=0 ∞

X x2 x4 x2n cos (x) = 1 − + − ... = (−1)n . 2! 4! (2n)! n=0 R = ∞. f : (−1, 1) → R, f (x) = ex

A sorok konvergencia sugara b) Az

függvényre

M =3

MacLaurin sorfejtése

∞

X xn x x2 e =1+ + + ... = . 1! 2! n! n=0 x

és a (2.3.21)


R = ∞. f : (−1, 1) → R, f (x) = (1 + x)α , α ∈ R − N

28

A sor konvergencia sugara c) Az

M=

függvényre

és a (2.3.21) MacLaurin sorfejtése

α α (α − 1) 2 x+ x + ... 1! 2! R = 1.

(1 + x)α = 1 + A sor konvergencia sugara

2.3.5. Hibabecslés

O−nagy

ordó segítségével

A függvények

aszimptotikus viselkedésének a tanulmányozására használjuk az úgy nevezett nagy-ordó

O,

kis-ordó

o

(Landau-féle) jelöléseket.

Va az aR egy környezete és f, g : Va → R két f = O (g) (f egyenl® nagy Ordó g -vel), ha ∃c > 0

35. Definíció. Ha függvény, akkor konstans ú.h.

|f (x)| ≤ c |g (x)| , x ∈ Va .

(2.3.23)

f (x) = 3x2 −x+5 akkor f = O (x2 ) mert |3x2 − x + 5| < 4x2 , ∀x ≥ 2. Hasonlóan sin (x) = O (1) , ∀x ∈ R. 36. Példa. Ha

A

O használata egyszer¶síti a függvények f ' f

közelítéséb®l létrejöv®

hiba tanulmányozását, ugyanis a hibafüggvényt egy egyszer¶bb (általában polinom) függvény viselkedésével helyettesítjük. 37. Definíció. Az meg az

f : Va → R

f : Va → R

függvény

O (xn )

renddel közelíti

függvényt ha

f (x) − f (x) ≤ c |xn | , x ∈ Va

(2.3.24) vagyis

hiba = f (x) − f (x) = O (xn ) . A közelítés, illetve az abból ered® hibarend az alábbi egyszer¶ kifejezéssel írható le:

f (x) = f (x) + O (xn ) .

(2.3.25)

2 x x3 + x2! + 1! 3! + x közelítésb®l ered® abszolút hiba: e − 1 + x

Például, a MacLaurin sorfejtésb®l

ex ' 1 + x +

x2 2

ex = 1 +

...

és az

+

x2 2

,

2.3. A SOROK KONVERGENCIÁJÁNAK A JAVÍTÁSA zéró környezetében kisebb mint

ex = 1 + x +

|x3 | (egy konstanstól eltekintve).

29 Tehát

x2 + O x3 , x ∈ V0 . 2

Hasonlóan

x x2 xn + + ... + + O xn+1 , x ∈ V0 , 1! 2! n! 3 x x x2n+1 sin (x) = − + ... + (−1)n + O x2n+3 , x ∈ V0 . 1! 3! (2n + 1)! ex = 1 +

Általánosan, a (2.3.18),(2.3.19) képletekb®l

f (x) = Tn (x) + Rn (x) = Tn (x) + O xn+1 ahol

Tn

az

f

függvény

n-ed

rend¶ Taylor polinomja.

A nagyordóra a következ® összeadási, illetve szorzási szabályok érvényesek:

• O (xp ) + O (xq ) = O (xr ) ahol r = min {p, q} ; • O (xp ) · O (xq ) = O (xs ) ahol s = p + q. xV0 : ex = 1 + 2 3 2 4 x + x2 + x6 + O (x4 ), cos (x) = 1 − x2 + x24 + O (x6) =⇒ 2 3 2 4 ex + cos (x) = 1 + x + x2 + x6 + 1 − x2 + x24 + O xmin{4,6} = 38. Példa. A MacLaurin sorfejtéseket használva

3

4

3

4

x 2 + x + x6 + x24 + O (x4 ) = 2 + x+ x6 + O (x4 ) mert + O (x4 ) = O(x4 ) ; 24 2 3 2 4 2 3 ex ·cos (x) = 1 + x + x2 + x6 · 1 − x2 + x24 + 1 + x + x2 + x6 O (x6 ) + x4 5 4 1 5 1 6 x2 x − 24 x + 48 x + + 1 − 2 + 24 O (x4 ) + O (x4+6 ) = (1 + x − 31 x3 − 24 1 1 3 7 min(4,6) 4 x )+O x = 1 + x − 3 x + O (x ) . 144 Összefoglalva, ha

r = min {p, q},

f (x) = f (x) + O (xp ) , g (x) = g (x) + O (xq ) ,

és

akkor:

• f (x) + g (x) = f (x) + g (x) + O (xr ) ; • f (x) g (x) = f (x) g (x) + O (xr ) ; (x) (x) = fg(x) + O (xr ) , (g (x) , g (x) 6= 0) . • fg(x) A sorozatok sajátos függvények, ezért az esetükben is használhatók a

O-ra

említett tulajdonságok.

∞ (xn )∞ n=1 , (yn )n=1 sorozatok esetében, (xn ) = egyenl® nagy O rend¶ (yn )), ha ∃c > 0 és N

39. Definíció. Adott

O ((yn ))

(az

(xn )

sorozat


30

számok úgy, hogy

|xn | ≤ c |yn | , ∀n ≥ N. n−1 és n2

xn = xn = O ((yn )) .

40. Példa. Ha tehát

yn =

1 akkor n

|xn | ≤ |yn | , ∀n ≥ 1,

(xn )n , limn xn = x∗ , konvergens konvergenciarendje p ∈ R, p ≥ 1, ha

41. Definíció. Az mondjuk, hogy

lim

(2.3.26)

n

sorozatról azt

|xn+1 − x∗ | = c, |xn − x∗ |p

vagy

en+1 = O (epn ) ,

(2.3.27) ahol

en = |xn − x∗ |

az

xn

hibája.

A konvergenciát lineárisnak nevezzük ha ha

p > 1,

és négyzetesnek ha

p = 1,

szuperlineárisnak

p = 2.

Az algoritmusok komplexitásának az elemzéséhez úgyszintén használatos

O jelölés. Igy például két n×n-es mátrix szorzásának a komplexitása O (n3 ) . Az O (np ) , p ≥ 1 algoritmusokat polinomiális típusúnak nevezzük (p = 1 lineáris, p = 2 négyzetes, p = 3 köbös). Gyakran fordul el® n még a logaritmikus O (log n) , loglineáris O (n log n) = O (log n ) = O (n!) , exponenciális O (cn ), c > 1, illetve faktoriális O (n!) típusú az

algoritmusok. A függvények aszimptotikus vizsgálatát más mér®számmal is elvégezhet®, például a kis-ordó feltétellel. Ha

f, g : Va → R

Az értelmezésb®l

f = O (g) .

(x) = 0. f = o (g) ha limx→a fg(x) látszik, hogy f = o (g) feltétel er®sebb

akkor

mint a

3. FEJEZET

Egyenletek numerikus megoldása Az

f (x) = 0, f : D ⊂ R −→ R

egyenleteket algebrainak, illetve

transzcendensnek nevezzük attól függ®en hogy az

f

polinom e (esetleg

azzá alakítható) vagy sem.

f (x) = 4x3 + 16x − 3 = 0 egyenlet sin x + 1 − x = 0 transzcendens egyenlet. Pl.

algebrai, míg az

f (x) =

Az egyenletek numerikus megoldása két lépésben történik: az els® lépés a gyökök elkülönítése (szeparálása), majd következik, a meghatározott intervallumokban, a gyökök kiszámítása.

3.1. A gyökök elkülönítése A gyökök elkülönítése egy olyan eljárás amellyel a

D

értelmezési

tartományt diszjunkt intervallumokra bontjuk úgy, hogy minden intervallum az egyenlet egyetlenegy gyökét tartalmazza. A módszerek között meg lehet említeni a grakus, illetve az analitikus eljárást.

A grakus módszer

Az eljárás a következ®: vázoljuk az

függvény görbéjét,

f (x) = 0 egyenlet közelít® gyökeit. Gyakran elonyosebb felbontani az f (x) = 0 egyenletet f1 (x) − f2 (x) = 0 alakra, es ezutan keresni az f1 illetve f2 majd az

Ox

y = f (x)

tengellyel való metszéspontok szolgálják az

fuggvenyek metszespontjat. Pl.

f (x) = sin x + 1 − x = 0.

El®nyösebb (és könnyebb) felbontani

sin x = x − 1. ∗ található: x ∈ (0, π).

az egyenletet a következ®képpen: két függvény metszeténél

Az egyenlet gyöke a

ABRA

Analitikus modszer A modszer a Rolle tetel ket folyomanyan alapszik.

31

3.2. A FELEZ MÓDSZER

Folyomány1.

32

a es b az f ∈ C 1 (D) fuggvenynek ket egymasutani 0 0 stacionarius pontja: f (a) = f (b) = 0, akkor f -nek legfeljebb egy gyöke van az (a, b) intervallumban. Valoban, ha ∃c1 , c2 ∈ (a, b) ugy hogy f (c1 ) = f (c2 ) = 0 akkor 0 Rolle tetele szerint ∃x ∈ (c1 , c2 ) u.h. f (x) = 0, de akkor a es b nem Ha

egymasutani.

Folyomány2. kovetkezo allitas: az

(a, b)

f ∈ C (D) folytonos fuggvenyre teljesul a f (a) f (b) < 0, akkor f -nek legalabb egy gyöke van Ha az

intervallumban.

Ha a ket folyomány igaz, akkor

(a, b)

f -nek

egyetlen egy gyöke lesz az

intervallumban.

Tehat: - meghatarozzuk az

f

fuggveny stacionarius pontjait:

x1 , x2 , ..., xn

es ezekkel -kepezzuk a Rolle-fele sorozatot:

f (x1 ) , f (x2 ) , ..., f (xn ) .

Ha a Rolle-fele sorozatban elojel valtas van, akkor az illeto intervallumban egyetlen egy gyöke van az egyenletnek.

2x3 − 9x2 + 12x − 4.5 = 0 egyenlet gyökeit. f (x) = 2x3 − 9x2 + 12x − 4.5 =⇒ f 0 (x) = 6x2 − 18x + 12 . f 0 (x) = 0 =⇒ x1 = 1 es x2 = 2 . A Rolle sorozatot tablazatba Pl. Kulonitsuk el az

irjuk:

x f (x)

-∞

1 2 +∞ -∞ 0.5 0.5 +∞ − + − +

Tehát az egyenlet gyökeit elkulonítettük a

(−∞, 1) , (1, 2) , (2, +∞)

intervallumokban.

3.2. A felez® módszer f (x) = 0 egyenletnek elkül®nítettük a gyökét az [a, b] intervallumban és f (a) f (b) ≤ 0, f ∈ C [a, b] . Az [a, b] intervallumot felbontjuk két egyforma hosszúságú részintervallumra: a+b a+b a, 2 , 2 , b . Ezek közül csak az egyik tartalmazza az x∗ gyököt Tételezzük fel, hogy az

(hacsak nem éppen

a+b a gyök, de akkor a feladat meg van oldva). Ezt 2

az intervallumot újból felezzük, stb.

3.2. A FELEZ MÓDSZER

33

ABRA Az algoritmus ket sorozatnak

(an )n∈N , (bn )n∈N

a felépítéséb®l áll.

a0 = a, b0 = b. a+b a+b Ekkor vagy f (a) f ≤ 0, vagy f (a) f 2 ≥ 0 ; ennek megfelel®en 2 a+b a1 = a és b1 = 2 , vagy a1 = a+b és b1 = b. Az n-edik lépésnél 2 ∗ meghatározzuk az [an , bn ] intervallumot mely tartalmazza az x gyököt A sorozatok tagjait a kovetkez®képpen értelmezzük:

es amit újból felezünk.

an +bn , ha pedig 2

bn+1 = bn+1 = bn .

n f (an ) f an +b ≤ 0 =⇒ an+1 = an 2 an +bn n ≥ 0 =⇒ an+1 = an +b f (an ) f 2 2

Ha

es és

Az (an )n sorozat monoton növekv®, a (bn )n sorozat pedig monoton csökken®. A sorozatok határértéke megegyezik: 42. Tétel.

lim an = lim bn = x∗ , ahol f (x∗ ) = 0. Bizonyítás. A sorozatok felépítéséb®l következik a monotonitás.

Mivel

an ≤ x ∗ ≤ b n

következik hogy a sorozatok korlátosak.

a sorozatok konvergensek.

Határértéket számítva a

képletb®l, a fogó szabályt alkalmazva, azt kapjuk,

lim an = x∗ .

Tehát

bn − an = b−a 2n hogy: lim bn =

A gyakorlatban természetesen csak véges lépést fogunk végrehajtani, viszont minden lépésben csak egy intervallumot kapunk (kivéve ha a gyök az intervallum közepére esik, ekkor viszont az algoritmus véget ér). A keresett gyök közelíthet® az

n-ik

x ∗ ≈ cn =

(3.2.1)

intervallum közepével

an + b n . 2

Mivel minden lépésben az intervallum felez®dik, az hossza

b−a . 2n hiba en kisebb

[an , bn ]

bn − an = vagyis,

n

lépés után az abszolút

(3.2.2) 43. Tétel.

en = |cn − x∗ | <

mint

b−a . 2n

A módszer konvergenciarendje lineáris en+1 = O (en ) .

intervallum

3.3. A HÚR MÓDSZER (REGULA FALSI) Bizonyítás.

en

és

34

en+1 -vel jelölve a hibákat az n, illetve (n + 1)-ik

intervallumban és ismerve, hogy az intervallumok felez®dnek következik, hogy

1 en+1 = en 2

vagyis

en+1 = O (en ) . Ha adott egy következik, hogy

> 0 pontosság, akkor a b−a ≤ egyenl®tlenségb®l 2n / ln (2) lépésre van szükség a pontosság n ≥ ln b−a

eléréséhez.

f (x) = sin x + 1 − x = 0. A [0, π] intervallumon az f π és f (0) = 1 > 0, f (π) = 1−π < 0, függvény el®jelt vált, közepe c1 = 2 π π π ∗ f 2 = 2 − 2 > 0, tehát x ∈ 2 , π mert itt vált el®jelt a függvény. A π π 3π 3π 3π ∗ , π intervallum k®zepe c2 = és f , < 0 tehát x ∈ . Ha a 2 4 4 2 4 π ∗ gyököt x ≈ c2 -vel közelítjük akkor a hiba kisebb mint e2 < 2 . Ahhoz, 2 π −4 hogy = 10 pontosságot elérjünk ln / ln (2) = 14.939 ≤ n = 15 10−4 44. Példa.

felezést kell végrehajtani.

f függvény nem vált el®jelt a f (x) = x2 = 0, x ∈ [−1, 3] egyenlet nem

A módszer nem alkalmazható ha az gyök két oldalán, például az

oldható meg felezési módszerrel.

3.3. A húr módszer (regula falsi) A felez® módszerhez hasonlóan, feltételezzük, hogy vagyis a gyök el van kül®nítve az

(a, b)

f (a) · f (b) < 0,

intervallumba. A húr módszer

alapgondolata hasonló a felez® módszeréhez, a különbség az, hogy a

(a, b) intervallum közepével közelítjük, hanem az A = (a, f (a)), B = (b, f (b)) pontokat összeköt® húr az Ox tengellyel való gyököt nem az

metszetével. ÁBRA A húr egyenlete

y =f (a) =

f (b)=f (a) (x=a), b−a

3.3. A HÚR MÓDSZER (REGULA FALSI)

35

ahonnan

f (a) c = a= f (b)=f (a) .

(3.3.1)

b−a

Az eljárást folytatjuk az

[a, b] := [a, c], vagy az [a, b] := [c, b] intervallummal

attól függ®en, hogy az els® vagy a második intervallumban találhatóe a gyök.

A felez® módszert®l eltér®en a húr módsdzernél az

intervallum nem feltétlenül konvergál nullához. látható, ahol

f (x) > 0

és

[a, b]

Ez az alábbi ábrán

f (b) > 0.

ÁBRA Ebben az esetben a (3.3.2)

B = (b, f (b))

pont rögzített, a gyök pedig a

xn+1 = xn −

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

sorozattal közelíthet® (x0

,

= a).

Numerikus eljárások esetében általánosan elfogadott az alábbi kilépési kritérium: az iterációt befejezzük ha az utolsó két tag abszolút:

|xn+1 − xn | < ,

(3.3.3) illetve relatív eltérése:

|xn+1 − xn | < , |xn |

(3.3.4) kisebb egy adott

>0

értéknél.

Sajátosan, mivel az egyenletek esetén (3.3.5)

f (xn ) ≈ 0,

használható az

|f (xn )| < ,

feltétel. A fenti kilépési feltételek nem garantálják, hogy az pontos gyökt®l kisebb mint

xn

eltérése a

.

A húr módszer konvergenciarendje lineáris.

f (x) = sin (x) + 1=x, x ∈ [0, π]. A húr metszete az Ox-vel c1 = 1 és f (0) = 1 > 0, f (π) = 1=π < 0, f (1) = sin(1) > 0, ∗ tehát x ∈ [1, π] mert itt vált el®jelt a függvény. Az [1, π] intervallumon a húr metszete az Ox-el c2 = 1.6. 45. Példa.

3.4. AZ ÉRINT (NEWTON) MÓDSZER

36

3.4. Az érint® (Newton) módszer f (x) = 0 egyenletnek, f ∈ C 1 ([a, b]) , elkül®nítettük a gyökét az [a, b] intervallumban. A módszer segítségével egy olyan (xn )n konvergens sorozatot állítunk ∗ el® amely konvergál az adott egyenlet x megoldásához. A módszer Tételezzük fel, hogy az

mértani jelentése és a sorozat el®állításának lépései az alábbi ábrán láthatók. ABRA

M0 (x0 , f (x0 )) pontból érint®t húzzunk a függveny görbéjéhez, ami az Ox tengelyt az x1 -ben metszi. Az érint®k szerkesztését folytatjuk az M1 (x1 , f (x1 )) ponttal, illetve az n-ik lépésben az Mn (xn , f (xn )) Legyen

x 0 = b.

Az

ponttal. Ebben a pontban az érint® egyenlete a következ®:

y − f (xn ) = f 0 (xn ) (x − xn ) . Az érint® metszete az

Ox

tengellyel adja az

(xn+1 , 0)

pontot:=⇒

−f (xn ) = f 0 (xn ) (xn+1 − xn ) ahonnan

xn+1 = xn −

(3.4.1)

f (xn ) . f 0 (xn )

3.4.1. A kiindulópont megválasztása és az érint® módszer konvergenciája Az iteratív módszereknél nagyon fontos a kezdeti pont helyes kiválasztva, ugyanis gyakran függ ett®l a módszer konvergenciája. A fenti ábrán is érzékelhet®: ha az valasztottuk volna, akkor A legegyszer¶bb ha

x1

x0

[a, b] intervallum masik végpontját

kívül esett volna ennek az intervallumon.

az

[a, b]

intervallum egyik végpontja, de ez

nem szükségszer¶. Feltételezzük, hogy

f ∈ C 2 ([a, b]) és a f 0 , f 00 deriváltak el®jeltartóak

az adott intervallumban. Ebben az esetben a következ® esetek lehetségesek: ABRAK

(3.4.2)

f (x0 ) f 00 (x0 ) > 0

3.4. AZ ÉRINT (NEWTON) MÓDSZER

37

Az érint® módszerrel (3.4.1) kapott sorozat xn+1 = xn − a fentiek szerint van meghatározva, konvergens.

46. Tétel. f (xn ) x0 f 0 (xn )

, ahol

Bizonyítás. Feltételezzük, hogy

eset tárgyalása hasonlóan történik).

f 0 (x) > 0, f 00 (x) > 0 (a többi Igazoljuk, hogy az (xn )n sorozat

monoton és korlátos, tehát konvergens. A (3.4.2)-ból következik, hogy

f (x0 ) > 0. Mivel f 0 (x) > 0 következik, hogy f növekv®, és negatív az [a, x∗ ) illetve pozitív az (x∗ , b] intervallumon, tehát x0 ∈ (x∗ , b] , ahol x∗ az f gyöke f (x∗ ) = 0. Feltételezve, hogy x∗ < xn ≤ b a Taylor képletb®l következik, hogy

1 1 0 = f (x∗ ) = f (xn )+ f 0 (xn ) (x∗ − xn )+ f 00 (ξn ) (x∗ − xn )2 , ξn ∈ (x∗ , xn ) 1! 2! 00 és mivel f (ξn ) > 0 következik, hogy xn+1 = xn − tehát

x∗

az

(xn )n

f (xn ) > x∗ , f 0 (xn )

sorozat alsó korlátja. Mivel

xn+1 − xn = − tehát az

f (xn ) > 0, ∀xn ∈ (x∗ , b]

f (xn ) < 0, f 0 (xn )

(xn )n sorozat monoton csökken®. λ-val jel®lve az (xn )n sorozat

határértékét a (3.4.1)-b®l

λ=λ− ahonnan

f (λ) = 0,

47. Tétel.

vagyis

f (λ) f 0 (λ)

λ = x∗ .

Az érint® módszer konvergenciarendje négyzetes en+1 = O e2n ,

ahol en = |xn − x∗ | . Bizonyítás. Taylor sorba fejtve az

f

függvényt következik, hogy

1 1 0 = f (x∗ ) = f (xn )+ f 0 (xn ) (x∗ − xn )+ f 00 (ξn ) (x∗ − xn )2 , ξn ∈ (x∗ , xn ) 1! 2! ∗ ahonnan gyelembe véve, hogy f (x ) = 0 következik, hogy xn −

f (xn ) f 00 (ξn ) ∗ ∗ − x = (x − xn )2 , ξn ∈ (x∗ , xn ) . f 0 (xn ) 2f 0 (xn )

3.5. A SZEL MÓDSZER

38

Felhasználva a Newton iterációt (3.4.1) következik, hogy

f 00 (ξn ) ∗ (x − xn )2 2f 0 (xn ) 00 f (ξn ) −−−−→ f 00 (ξn ) |xn+1 − x∗ | = n → ∞ 0 ∗ = c, (x∗ − xn )2 2f 0 (xn ) 2f (x ) xn+1 − x∗ =

⇒

vagyis

en+1 = O (e2n ) .

f (x) = sin x+1−x = 0, x ∈ (0, π) ⇒ f 0 (x) = cos x−1, f 00 (x) = − sin x < 0 ha x ∈ (0, π) , tehát x0 = π. Az (3.4.1) iterációból 0) 1) ⇒ x1 = x0 − ff0(x = π − 1−π ≈ 1.940 2. ≈ 2. 070 8, x2 = x1 − ff0(x (x0 ) −2 (x1 ) 48. Példa.

n-edik lépésben az iterálást abba lehet hagyni ha a függvény értéke f (xn ) < , vagy két egymásutáni tag abszolút |xn − xn−1 | (relatív |xn −xn−1 | ) eltérése elég kicsi (< ). |xn | Az

49. Példa. A Newton módszer jól alkalmazható egész kitev®j¶ gyökvonásra

√ p

a.

Például

√

a

kiszámítása ekvivalens az

x2 − a = 0

egyenlet (pozitív) megoldásával amit Newton módszerrel határozunk meg. A sorozat rekurzíós képlete a következ® lesz: (3.4.3)

ahonnan

1.414 2,

x2n − a 1 a xn+1 = xn − = xn + 2xn 2 xn √ √ √ a = 2-re ⇒ 2 ≈ x0 = 2, 2 ≈ 23 , 2 ≈

17 , 12

√

2≈

577 408

=

...

A Newton módszer hátránya, hogy minden iterációban szükséges az

f

függvény deriváltja. A következ® módszerek ezt a hátrányt küszöbölik

ki a konvergenciarend róvására.

3.5. A szel® módszer A (3.4.1) képletben az

f 0 (xn )

f 0 (xn ) ≈

tagot felcserélve a közelít® értékkel

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

a következ® iterációs képletet kapjuk: (3.5.1)

xn+1 = xn −

f (xn ) , f (xn )−f (xn−1 ) xn −xn−1

n ≥ 1, x0 , x1 =kezdeti

értékek

3.6. A STEFFENSEN-FÉLE MÓDSZER

39

A mértani jelentése a szel® módszernek a következ®: kiindulva kezdeti értékekb®l az

Ox

szel® metszete az

(xn−1 , f (xn−1 )) , (xn , f (xn ))

x0 , x1

pontokon áthaladó

tengellyel adja a közelít® sorozat

xn+1

tagját.

ABRA Az

(xn )n

sorozat els® tagját

x0 -t

ugyanúgy választjuk meg mint

x∗

x1 ∈ (x , x0 ) . Az x0 , x1 pontoknak megfelel a függvény görbéjén az M0 (x0 , f (x0 )) , és M1 (x1 , f (x1 )) . Az M0 M1 szel® metszi az Ox tengelyt x2 -ben. Az eljárást folytatjuk M1 M2 szel®vel. Az n-edik lépésnél Mn−1 Mn ∩ Ox = xn+1 . Az Mn−1 Mn szel® egyenlete:

az érint® módszernel (3.4.2),

x1

pedig egy pont

x0

es

között:

∗

f (x) − f (xn ) =

f (xn ) − f (xn−1 ) (x − xn ) xn − xn−1

amelybe behelyettesítjük a metszéspont koordinátáit:

−f (xn ) = ahonnan kifejezve

(xn+1 , 0) =⇒

f (xn ) − f (xn−1 ) (xn+1 − xn ) xn − xn−1

xn+1 -t

kapjuk a (3.5.1) képletet. √ 1+ 5 (< A szel® módszer konvergenciarendje p = 2

2).

3.6. A Steensen-féle módszer f 0 (xn ) derivált közelítéséhez felhasználjuk, hogy a gyök közelében f (xn ) ≈ f (x∗ ) = 0, tehát a derivált Az

f (xn + h) − f (xn ) f (xn + h) − f (xn ) ≈ h→0 h h

f 0 (xn ) = lim

átírható a következ® alakban:

f 0 (xn ) ≈

f (xn + f (xn )) − f (xn ) . f (xn )

Innen az iteratív eljárás a következ® lesz: (3.6.1)

xn+1 = xn −

f (xn ) . f (xn +f (xn ))−f (xn ) f (xn )

A Steensen módszer konvergenciarendje (a gyök környezetében) négyzetes p = 2. 50. Tétel.

3.7. A FOKOZATOS KÖZELÍTÉSEK MÓDSZERE

40

3.7. A fokozatos közelítések módszere A fokozatos közelítések (sukcesszív approximációk) módszere a numerikus eljárások egyik alapvet® és gyakran használt módszere. A megoldandó egyenletet

f (x) = 0

átírjuk a következ® alakra:

x = ϕ (x) .

(3.7.1)

ϕ (x) = x + f (x), ϕ (x) = x − f (x), ϕ (x) = x + ωf (x) .

Erre több lehet®ség is van, például vagy általánosabban Az

(xn )n

sorozatot ahol az approximációk sorozatának nevezzük.

Ha az (xn )n sorozat konvergens es ϕ folytonos akkor a határértéke megegyezik az egyenlet gyökével. 51. Tétel.

x∗ = lim xn+1 = lim ϕ (xn ) = ϕ (lim xn ) = ϕ (x∗ ) ∗ vagyis x gyöke az (3.7.1) egyenletnek ami ekvivalens f (x) = 0 egyenlettel. Bizonyítás.

A

ϕ-t®l

függ®en az eljárás lehet konvergens vagy divergens.

Az

alábbi ábrákon szemléltettünk egy pár esetet: ABRAK

y = x illetve y = ϕ (x) görbék metszeténél taláható. Kiindulva az (xn , xn+1 ) = (xn , ϕ (xn )) pontból a (xn+1 , xn+2 ) pontot a következ® keppen szerkesztjük meg. Az (xn , xn+1 ) ponton keresztül húzzuk az y = xn+1 egyenest az y = x egyenessel való metszésig, majd berajzolva az x = xn+1 egyenest az y = ϕ (x) görbével való metszésig megkapjuk a kivánt (xn+1 , xn+2 ) pontot. A gyök az

Legyen ϕ : [a, b] → [a, b] egy deriválható függvény. Ha teljesül a következ® feltétel: 52. Tétel.

|ϕ0 (x)| ≤ q < 1, ∀x ∈ (a, b)

akkor a) az (xn )n sorozat: xn+1 = ϕ (xn ) konvergens tetsz®leges x0 ra, b) a sorozat x∗ határértéke az ϕ (x) = x egyenlet egyetlen gyöke az [a, b] intervallumban.

3.7. A FOKOZATOS KÖZELÍTÉSEK MÓDSZERE

41

Bizonyítás. a)

|xn+1 − xn | =def |ϕ (xn ) − ϕ (xn−1 )| =Lagrange |ϕ0 (ξn ) (xn − xn−1 )| ≤ ≤T etel q |xn − xn−1 | ≤ ... ≤ q n |x1 − x0 | , ahol

ξn ∈ (xn , xn+1 ) =⇒

|xn+p − xn | = |xn+p − xn+p−1 + xn+p−1 − ... − xn | ≤ |xn+p − xn+p−1 | + ... + |xn+1 − xn | ≤ 1 − qp |x1 − x0 | . ≤ q n+p−1 |x1 − x0 | + ... + q n |x1 − x0 | = q n+p−1 + ... + q n |x1 − x0 | = q n 1−q Tehát: (3.7.2)

|xn+p − xn | ≤

qn |x1 − x0 | , n, p ∈ N. 1−q

q ∈ (0, 1) a jobboldal tetsz®legesen lekicsinyíthet® vagyis a (xn )n Cauchy-féle sorozat. Az R-en minden Cauchy sorozat konvergens, tehát (xn )n konvergens.

Mivel

(xn )n sorozat határértékét x∗ -al: limxn = x∗ =⇒ x∗ = lim xn = lim ϕ (xn−1 ) = ϕ (lim xn−1 ) = ϕ (x∗ ) tehát x∗ az x = ϕ (x) egyenlet gyöke. Kimutatjuk hogy x∗ az egyeduli gyöke. ∗∗ Feltételezzük, hogy az egyenletnek még létezik egy gyöke: x = ϕ (x∗∗ ) , x∗∗ 6= x∗ =⇒ |x∗ − x∗∗ | = |ϕ (x∗ ) − ϕ (x∗∗ )| = |ϕ0 (ξ ∗ ) (x∗ − x∗∗ )| ≤ q |x∗ − x∗∗ | , ahol ξ ∗ ∈ (x∗ , x∗∗ ) =⇒ (1 − q) |x∗ − x∗∗ | ≤ 0 és mivel (1 − q) > 0 =⇒ |x∗ − x∗∗ | = 0 ⇔ x∗ = x∗∗ hamis. A (3.7.2) képletb®l ha p → ∞ kapjuk a következ® hibabecslást: b) Jelöljük az

|xn − x∗ | ≤

(3.7.3)

qn |x1 − x0 | . 1−q

Hasonlóan kimutatható, hogy: (3.7.4)

|xn − x∗ | ≤

53. Példa. Az

q |xn − xn−1 | , n ≥ 1. 1−q

x3 + 12x − 1 = 0

egyenletet aminek a gyöke a 3

[0, 1] intervallumba esik, átírhatjuk többféleképpen: x = 1−x , x = 12 √ 3 1 1−x3 0 1 − 12x, vagy x = x2 +12 . Az els® esetben ϕ (x) = 12 és |ϕ (x)| = √ 3 x2 1 ≤ , x [0, 1] ; a második esetben pedig a ϕ (x) = 1 − 12x függvényre 4 4 0 nem érvényes |ϕ (x)| < 1. Ha az utolsó esetet vesszük gyelembe akkor

3.8. VEKTOR ÉS MÁTRIX NORMÁK

42

1 ϕ (x) = x2 +12 , ϕ : [0, 1] → [0, 1] es |ϕ0 (x)| < 1. Meghatarozzuk q -t: 2x 2 q = supx∈[0,1] |ϕ0 (x)| = supx∈[0,1] (x2 +12) 2 = 169 . Kezd®értéknek legyen 144 1 , x2 = ϕ (x1 ) = 1729 . Ha el®re meg van x0 = 0 =⇒ x1 = ϕ (x0 ) = 12 −4 ∗ adva a hibakorlat ε = 10 akkor (3.7.3)-ben kikötjük |xn − x | ≤ n q |x1 − x0 | < ε =⇒ n = 2 vagyis x2 megközelíti a gyököt a kívánt 1−q pontossággal. A

| ϕ0 (x)| < 1

függvény meghatározására nem létezik általános

eljárás, de egyes esetekben ennek el®állítására bizonyos sémát követhetünk. Például ha

f (x) = 0

0 < m ≤ f 0 (x) ≤ M

egyenlet ekvivalens az

M ≤ f 0 (x) ≤ m < 0) akkor az x = x − M1 f (x) egyenlettel, tehát ha (vagy

ϕ (x) = x − 0 ≤ ϕ0 (x) = 1 −

f 0 (x) M

≤1−

1 f (x) , M

m M

< 1. Az ismertett módszereknél az iteraló ϕ függvény a következ®: Newton f (x) f (x) , Steensen módszer ϕ (x) = x − f (x+f (x))−f (x) , módszer ϕ (x) = x − 0 f (x)

akkor

f (x)

szel® módszer

ϕ (x, y) =

yf (x)−xf (y) stb. f (x)−f (y)

Mivel a fokozatos kozelítések módszerét nem csak az

R térben alkalmazzuk,

a továbbiakban ismertetjük a fenti tétel általanosítását normált terekben.

3.8. Vektor és mátrix normák A norma a modulusz (abszolút érték) fogalom általánosítása. 54. Definíció. Az

n

nevezzük (R

||·|| : Rn → R

függvényt vektor normának

térben) ha eleget tesz az alábbi tulajdonságoknak bármilyen

n

x, y ∈ R és α ∈ R: (i) ||x|| ≥ 0 és ||x|| = 0 ⇔ x = 0Rn ; (ii) ||αx|| = |α| ||x|| ; (iii) ||x + y|| ≤ ||x|| + ||y|| . Egy vektort egységvektornak nevezzünk ha

||x|| = 1.

Az egyik leggyakrabban használt az úgynevezett 1

p−norma:

||x||p = (|x1 |p + ... + |xn |p ) p , p ≥ 1, x = (x1 , ..., xn )t ∈ Rn


p = 1, p = 2

ezek között is a

és

p=∞

43

a legismertebbek:

||x||1 = |x1 | + ... + |xn | , 1 ||x||2 = x21 + ... + x2n 2 ,

(3.8.1) (3.8.2)

||x||∞ = max (|x1 | , ..., |xn |) .   1   3 55. Példa. Ha x =  −5  ∈ R , akkor ||x||1 = 8, ||x||2 = 2 √ 30, ||x||∞ = 5.

(3.8.3)

56. Tétel.

egyenl®tlenség:

Ha p1 + 1q = 1, p > 1, q ∈ R akkor igaz az alábbi (Hölder) t x y ≤ ||x|| ||y|| . p q

(3.8.4)

Fontos következménye az el®bbi tételnek az úgynevezett CauchySchwartz egyenl®tlenség:

t x y ≤ ||x|| ||y|| 2 2 vagy komponensekre lebontva: (3.8.5)

|x1 y1 | + ... + |xn yn | ≤

q

x21

+ ... +

x2n

q

y12 + ... + yn2 .

A Rn -n értelmezett normák ekvivalensek, vagyis ha ||·||p és ||·||q normák az Rn térben akkor léteznek c1 , c2 pozitív konstansok úgy, hogy: 57. Tétel.

c1 ||x||p ≤ ||x||q ≤ c2 ||x||p , ∀x ∈ Rn . Például, ha

x ∈ Rn

akkor:

√ n ||x||2 , √ ≤ ||x||2 ≤ n ||x||∞ ,

||x||2 ≤ ||x||1 ≤ ||x||∞

||x||∞ ≤ ||x||1 ≤ n ||x||∞ . Ha

x ∈ Rn

az

X ∈ Rn

vektor approximációja, akkor a közelítés

abszolút, illetve relatív hibája:

abs = ||X − x|| , rel =

||X − x|| . ||x||


44

A mátrix norma deníciója hasonló a vektor normáéhoz.

||·|| : Rm×n → R függvényt mátrix normának m×n nevezzük (Mmn (R) ≈ R térben) ha eleget tesz az alábbi tulajdonságoknak m×n bármilyen A, B ∈ R és α ∈ R: (i) ||A|| ≥ 0 és ||A|| = 0 ⇔ A = 0Rm×n ; (ii) ||αA|| = |α| ||A|| ; (iii) ||A + B|| ≤ ||A|| + ||B|| . 58. Definíció. Az

Ha egy mátrix norma teljesíti a

||AB|| ≤ ||A|| ||B|| , A, B ∈ Rn×n .

(3.8.6)

tulajdonságot, akkor (a normát) szubmultiplikatívnek nevezzük. A leggyakrabban használt mátrix normák a Frobenius:

v uX n p u m X t ||A||F = a2ij = T r (At A),

(3.8.7)

i=1 j=1 illetve egy

x ∈ Rn ,

(3.8.8)

vektor által származtatott (indukált)

||A||p = sup x6=0

p−norma:

||Ax||p . ||x||p

A (3.8.8)-ból következik, hogy:

(3.8.9)

||A||p = sup A x6=0

x ||x||p

! = max ||Ay||p . ||y||p =1 p

p-normák kapcsolódnak a vektor p-normákhoz, következik, hogy úgyszintén a p = 1, 2, ∞ esetek a legfontosabbak. m×n Ha A ∈ R akkor m X (3.8.10) ||A||1 = max |aij | , Mivel a mátrix

1≤j≤n

(3.8.11)

||A||∞ =

max

1≤i≤m

i=1 n X

|aij | .

j=1

p = 2-re a ||A||22 norma egyenl® a p (λ) = det (At A − λI) karakterisztikus t polinom legnagyobb gyökével (vagyis az A A legnagyobb sajátértékével). Ha csak az ||A||2 norma nagyságrendje fontos, akkor használhatjuk az


45

alábbi egyenl®tlenségeket:

√ n ||A||2 , √ max |aij | ≤ ||A||2 ≤ mn max |aij | , ||A||2 ≤ ||A||F ≤

i,j

(3.8.12)

i,j

√ 1 √ ||A||∞ ≤ ||A||2 ≤ m ||A||∞ , n √ 1 √ ||A||1 ≤ ||A||2 ≤ n ||A||1 , m q ||A||2 ≤ ||A||1 ||A||∞ . 

 2 −1 √   59. Példa. A =  1.5 3  =⇒ ||A||F = 4 + 1 + 2.25 + 9 + 4 = 0 −2 4.5, ||A||1 = max {3.5, 6} = 6, ||A||∞ = max {3, 4.5, 2} = 4.5, 3 ≤ √ ||A||2 ≤ 3 6 (||A||2 ' 3.8388) . Az ismertetett normák mindegyike szubmultiplikatív tulajdonságú, viszont a

||A|| := maxi,j {|aij |}

norma nem.

3.8.1. Banach-féle xpont tétel normált terekben 60. Definíció. A a

||·||

normára nézve,

Φ : Rn → Rn leképezést kontrakciónak nevezzük ha ∃q ∈ (0, 1) ú.h.

||Φ (X) − Φ (Y )|| ≤ q ||X − Y || , ∀X, Y ∈ Rn .

Ha Φ : Rn → Rn leképezés kontrakció, akkor Φ-nek egyetlenegy X ∗ ∈ Rn xpontja van: Φ (X ∗ ) = X ∗ . Az X (k+1) = Φ X (k) iterációk sorozata konvergál bármilyen X (0) ∈ Rn kezd®értékre. A hibára a következ® becslések igazak: 61. Tétel.

(3.8.13) (3.8.14)

qk 1−q (k) q X − X ∗ ≤ 1−q (k) X − X ∗ ≤

(1) X − X (0) , (k) X − X (k−1) , k ≥ 1.

3.9. ALGEBRAI EGYENLETEK NUMERIKUS MEGOLDÁSA

46

3.9. Algebrai egyenletek numerikus megoldása 3.9.1. Polinomok alakjai

Bármely

n-ed

fokú

pn ∈ Pn

polinom

felírható kanonikus alakban:

pn (x) = an xn + an−1 xn−1 + ... + a1 x + a0 =

(3.9.1)

n X

ai x i ,

i=0 ami azt jelenti hogy a

(1, x, ..., xn−1 , xn )

polinom-rendszer egy bázist

alkot, vagyis a tagok lineárisan függetlenek. Különböz® okokból indokolt más alakban is felírni a polinomokat. Egy ilyen alak az úgynevezett Bézier-féle polinom alak. 62. Definíció. A Bernstein-féle alap-polinomokat a következ®képpen értelmezük:

bni

(3.9.2)

Az

[0, 1]

n i (x) = x (1 − x)n−i , i = 0, ..., n; x ∈ [0, 1] . i

intervallum nem jelent megszoritást mivel az

x= változócserével az 63. Tétel. (3.9.3)

y−a b−a

[a, b] 3 y intervallum átalakítható [0, 1] intervallumá.

A bni polinomokra a következ® rekurrens képlet igaz:

bni (x) = (1 − x) bin−1 (x) + xbn−1 i−1 (x) , i = 1, ..., n − 1 n−1 bn0 (x) = (1 − x) b0n−1 (x) , bnn (x) = xbn−1 (x) .

Bizonyítás. Az értelmezést felhasználva kapjuk, hogy:

bn0

n 0 n−1 0 n (x) = x (1 − x) = (1 − x) x (1 − x)n−1 = (1 − x) b0n−1 (x) . 0 0

bnn (x) = xbn−1 n−1 (x) . Ha i = 1, ..., n − 1 felhasználva a kombinaciók rekurzív képletet kapjuk: n i n−1 n−1 n−i n bi (x) = x (1 − x) = + xi (1 − x)n−i = i i i−1 n−1 i n − 1 i−1 n−i−1 = (1 − x) x (1 − x) +x x (1 − x)n−i = i i−1 Hasonlóan

= (1 − x) bn−1 (x) + xbn−1 i i−1 (x)


47

64. Tétel.

A bn0 , bn1 , ..., bnn polinomok lineárisan függetlenek.

Bizonyítás. Legyen

c0 , c1 , ..., cn ∈ R

ú.,h.:

c0 bn0 (x) + c1 bn1 (x) + ... + cn bnn (x) = 0 , x ∈ [0, 1] .

(3.9.4)

n = 0 akkor c0 b00 (x) = 0 hogy c0 = 0. Ha

és mivel (3.9.2)-bol

b00 (x) = 1

azt kapjuk

n ≥ 1. Ha x = 0 (3.9.4)-ból azt kapjuk, hogy c0 = 0 , 1 n hasonlóan x = 1-re cn = 0. =⇒ c1 x (1 − x)n−1 +c2 n2 x2 (1 − x)n−2 ...+ 1 n−1 n cn−1 n−1 x (1 − x) = 0 , x ∈ [0, 1] . Újból x = 0-ra c1 = 0, míg x = 1-re cn−1 = 0 stb.=⇒ ci = 0, i = 0, ..., n tehát bni lineárisan függetlenek. n n n A fenti tételb®l következik, hogy a (b0 , b1 , ..., bn ) rendszer egy bázist alkot a Pn térben, tehát bármilyen p ∈ Pn polinomot egyértelm¶en fel n lehet írni a bi polinomok lineáris kombinációjáként: Legyen

p=

(3.9.5)

c0 bn0

+

c1 bn1

+

...cn bnn

=

n X

ci bni .

i=0 A fenti alakot Bézier-féle alaknak nevezzük,

ci pedig a p polinom Bézier

együtthatói.

n = 3-ra a harmadfokú alap-polinomok: (1 − x)3 , 3x (1 − x)2 , 3x2 (1 − x) , x3 , egy harmadfokú polinom Bézier-féle alakja: (3.9.6)

p (x) = c0 (1 − x)3 + c1 3x (1 − x)2 + c2 3x2 (1 − x) + c3 x3 .

65. Példa. Legyen

c3 = p (1) = 1.

p ∈ P3 , p (x) = 4x3 − 3x. =⇒ c0 = p (0) = 0,

Innen kapjuk, hogy

p (x) = c1 3x (1 − x)2 + c2 3x2 (1 − x) + x3 c1 , c2 tagokat c1 = −1, c2 = −2.

majd a

tartalmazó egyenletrendszert megoldva kapjuk


48

A Bézier-féle alakban lév® polinom deriváltját a következ® képlet adja meg:

p0 (x) = n

n−1 X

(ci+1 − ci ) bin−1 (x) .

i=0 Bizonyítás.

n X

!0 ci bni

(x)

i=0

=

n X

ci (bni

i=0

h n i0 X n (x)) = ci xi (1 − x)n−i = i i=0 0

h n i X n ci ixi−1 (1 − x)n−i − (n − i) xi (1 − x)n−i−1 = i i=0 n X i=1

n−1

X n! n! ixi−1 (1 − x)n−i − ci (n − i) xi (1 − x)n−i−1 = ci i! (n − i)! i! (n − i)! i=0

n X

n−1 X (n − 1)! (n − 1)! n−i i−1 ci n x (1 − x) − ci n xi (1 − x)n−i−1 = (i − 1)! (n − i)! i! (n − i − 1)! i=1 i=0 ! n n−1 X X n − 1 n − 1 i−1 n ci x (1 − x)n−i − ci xi (1 − x)n−i−1 =i−1=j i i − 1 i=1 i=0 ! n−1 n−1 n−1 X X X n−1 n−1 n cj+1 bj (x) − ci bi (x) = n (ci+1 − ci ) bin−1 (x) . j=0

i=0

i=0

Igényekt®l függ®en más polinom alakról is beszélhetünk, például Lagrange, Hermite, Csebysev, Legendre stb. polinomokról.

3.9.2. A Horner, illetve de Casteljau elrendezés

Adott polinomra

fontos meghatározni minél pontosabban a polinom értékét bizonyos pontokban. Erre szolgál a Horner séma ha a polinom kanonikus alakban van megadva.

p (x) = a0 + a1 x + ... + an xn egy polinom, illetve egy t ∈ R. p (t) = a0 + a1 t + ... + an−1 tn−1 + an tn értéket a következ® rekurzív Legyen

A

képlettel számítjuk ki: (3.9.7)

p (t) = a0 + t (a1 + ...t (an−1 + tan )) .


49

Gyakorlatban a fenti Horner elrendezés a következ®képpen valósul meg:

q n = an ,

(3.9.8)

qn−1 = an−1 + tqn , ... q0 = a0 + tq1 , p (t) = q0 .

tehát

A Horner elrendezést táblázatba is lehet írni:

an an−1 an−2 t an an t + an−1 (an + tan−1 ) t + an−2 Ha

p (t) = 0

akkor

t

... ...

a1 a0 . a1 + ... p (t)

gyöke a polinom függvénynek.

Felhasználva a polinomok maradékkal való osztási képletét felírhatjuk:

p (x) = (x − t) q (x) + r, 66. Példa. Legyen

ahol

r

az osztási maradék

p (x) = 4x3 − 3x

és

=⇒ p (t) = r.

t = 1/4. =⇒

4 0 −3 0 1/4 4 1 −11/4 −11/16 Tehát

p (1/4) = −11/16.

A de Casteljau elrendezés ugyanazt a célt szolgálja mint a Horner elrendezés, csak a Bézier alakban felírt polinomok esetében használjuk. Legyen

p (x) =

c0 bn0

(x) +

c1 bn1

(x) +

...cn bnn

(x) =

n X

ci bni (x)

i=0 egy Bézier alakban felírt polinom. Akkor felhasználva a (3.9.3) rekurzív képleteket következik, hogy:

p (t) = c0 bn0 (t) + c1 bn1 (t) + ... + cn−1 bnn−1 (t) + cn bnn (t) = = c0 (1 − t) bn−1 (t) + c1 (1 − t) b1n−1 (t) + tbn−1 (t) + ...+ 0 0 n−1 n−1 +cn−1 (1 − t) bn−1 (t) + tbn−2 n−1 (t) + cn tbn−1 (t) = (t) + ((1 − t) c1 + tc2 ) b1n−1 (t) + ... + ((1 − t) cn−1 + tcn ) bn−1 ((1 − t) c0 + tc1 ) bn−1 0 n−1 (t) .


50

Jelöljük az új együtthatókat:

(1)

ci (3.9.9)

= (1 − t) ci + tci+1 , i = 0, ..., n − 1 =⇒

:

n−1 X

p (t) =

(1)

ci bn−1 (x) . i

i=0 Megismételve az eljárást:

(2)

ci

(1)

(1)

= (1 − t) ci + tci+1 , i = 0, ..., n − 2 =⇒

:

p (t) =

n−2 X

(2)

ci bin−2 (x) ,

i=0 és az utolsó iterációban:

(n)

c0

(n−1)

= (1 − t) c0

:

p (t) =

(n−1)

+ tc1

=⇒

(n) c0 .

Tehát, kezdve az el®zetes együtthatókkal:

(0)

ci := ci , i = 0, ..., n kiszámítjuk a lineáris kombinációkat: (3.9.10)

(k)

ci

(k−1)

:= (1 − t) ci

az utolsó érték adja meg a

p (t)

(k−1)

+ tci+1 , i = 0, ..., n − k;

értékét.

Gyakorlatban az adatokat a következ® táblázat szerint rendezzük:

(0)

c0 (0) (1) c1 c0 (0) (1) c2 c1

c0

... (0) cn

... (2) cn−2

... (1) cn−1

(2)

...

(n)

c0

p polinom eredeti együtthatóit, a második oszlop pedig a (3.9.10) képlet szerint vannak kiszámítva k = 1-re. A (n) táblázat jobb-alsó sarkában lév® c0 érték adja a p (t) értéket. Az els® oszlop tartalmazza a


51

p az el®bbi peldából p ∈ P3 , és a Bézier együtthatók: c0 = 0, c1 = −1, c2 = −2, c3 = 1. A polinom értéke t = 1/4-ben a következ®: Legyen

0 −1 − 41 −2 − 54 − 12 1 − 45 − 54 − 11 16 Tehát

p (1/4) = −11/16.

3.9.3. Algebrai egyenlet megoldása Newton-Horner módszerrel A Newton-Horner módszer egy algebrai egyenlet valós gyökeinek a meghatarozására szolgál. Ennek érdekében alkalmazzuk a Newton-féle módszert kihasználva a polinomok tulajdonságait. Tekintsük a következ® algebrai egyenletet:

pn (x) = an xn + an−1 xn−1 + ... + a1 x + a0 = 0. Felhasználjuk a Newton módszerb®l ismert iteráló sorozatot:

xi+1 = xi −

pn (xi ) , i≥0 p0n (xi )

ahol kezdeti értéknek fel lehet használni például az

x0 =

− aa01 . A

x0 = 0,

vagy

pn (xi ) értéknek a meghatározására a Horner sémát fogjuk

használni:

pn (xi ) = a0 + xi (a1 + ...xi (an−1 + xi an )) . A polinomokra vonatkozó maradékos osztás tétele szerint:

pn (x) = (x − xi ) pn−1 (x) + R,

(3.9.11) ahol

pn−1 ∈ Pn−1 és R = pn (xi ) az (x − xi )-val osztási maradék.

Jelöljük

pn−1 (x) := bn−1 xn−1 + bn−2 xn−2 + ... + b1 x + b0 . Mint ismeretes a

bj , j = 0, ..., n − 1

együtthatók kiszámíthatóak a

(3.9.8) Horner sémából :

bn−1 = an bn−j = an−j+1 + xi bn−j+1 , j = 2, ..., n.


52

Deriválva a (3.9.11) képletet kapjuk hogy:

p0n (x) = (x − xi ) p0n−1 (x) + pn−1 (x) ahonnan:

p0n (xi ) = pn−1 (xi ) .

(3.9.12)

xi

A a

bj

an an−1 an−2 ... a1 a0 an an xi + an−1 (. . .) xi + a1 (. . .) xi + a0 bn−1 bn−2 b0

pn−1 (xi ) érték kiszámításához újból felhasználjuk a Horner sémát

együtthatókra:

pn−1 (xi ) = b0 + xi (b1 + ... + xi (bn−2 + xi bn−1 )) . Az eljárást abbahagyjuk ha a kapott gyök elég pontos. Ha meghatároztuk az egyik gyököt

(α)

elosztjuk a

pn

polinomot

(x − α)-val

es a kapott

n−1-ed fokú polinomra újból alkalmazzuk a Newton-Horner módszert. Legvégül másodfokú polinomot kapunk aminek meghatározzuk a gyökeit. 67. Példa. Határozzuk meg az alábbi egyenlet valós gyökeit (x0

0): x4 − 10x3 + 35x2 − 50x + 24 = 0. x1 = x0 −

24 p4 (0) =0+ = 0.48, |x1 − x0 | = 0.48 0 50 p4 (0)

1 −10 35 −50 24 0 1 −10 35 −50 24 0 1 −10 35 −50 x2 = x1 −

p4 (0.48) 7.011 = 0.48 − = 0.78, |x2 − x1 | = 0.3 ... 0 p4 (0.48) −22.869

1 −10 35 −50 24 0.48 1 −9.52 30.43 −35.39 7.011 . 0.48 1 −9.04 26.09 −22.869

=

4. FEJEZET

Lineáris egyenletrendszerek numerikus megoldása Cramer det

n! 100! = 10157.9

3 Pontos módszerek: ismeretlenek szama kisebb mint10 , Gauss, faktorizacios t modszerek(LU,LL ,QR);

103
(4.0.13) Ax

= b

A = (aij )i,j=1,n , x = (x1 x2 ... xn )t , b = (b1 b2 ... bn )t

4.1. Sajátos alakú egyenletrendszerek Az

(4.1.1)

  a11 x1 + a12 x2 + ... + a1n xn = b1     a22 x2 + ... + a2n xn = b2 ..  .     a x = b nn n

n

sajátos lineáris egyenletrendszert fels®háromszög¶nek nevezzük ha az egyenletrendszer

A

mátrixában a f®átló alatt csak nullák vannak:



a11 a12 ... a1n   0 a22 ... a2n A=  ..  . 0 0 ... ann 53

   .  

4.2. DIREKT MÓDSZEREK Ha

aii 6= 0, i = 1, n,

54

akkor a (4.1.1) megoldását visszahelyettesítéssel

kapjuk meg:

bn , ann P bk − ni=k+1 aki xi = , k = n − 1, 1. akk

xn = (4.1.2)

xk

Az

(4.1.3)

  a11 x1 = b1     a21 x1 + a22 x2 = b2 ..  .     a x + a x + ... + a x = b n1 1 n2 2 nn n n

sajátos lineáris egyenletrendszert alsóháromszög¶nek nevezzük ha az egyenletrendszer

A

mátrixában a f®átló fölött csak nullák vannak:



a11 0 ... 0   a21 a22 ... 0 A=  ..  . an1 an2 ... ann Ha

aii 6= 0, i = 1, n,

  .  

akkor a (4.1.3) megoldása:

x1 = (4.1.4)



xk =

b1 , a11 bk −

Pk−1 i=1

akk

aki xi

, k = 2, n.

4.2. Direkt módszerek 4.2.1. Gauss módszer

A Gauss módszer a lineáris egyenletrendszerek

numerikus megoldásának egyik legrégibb sémája. Lineáris transzformációkat alkalmazva a sorok között, a (4.0.13) egyenletrendszert visszavezetjük a (4.1.1) háromszög alakú egyenletrendszerre.

4.2. DIREKT MÓDSZEREK Az egyszer¶ség kedvéért az eljárást egy

55

4×4-es Ax = b egyenletrendszeren

mutatjuk be:

 (4.2.1)

Ha

a11 6= 0,

   

a11 a21 a31 a41

a12 a22 a32 a42

a13 a23 a33 a43

a14 a24 a34 a44



x1 x2 x3 x4

   





    =  

b1 b2 b3 b4

   . 

akkor az els® sort rendre

−

(4.2.2)

a31 a41 a21 =: e l21 , − =: e l31 , − =: e l41 a11 a11 a11

értékekkel szorozzuk (a szabadtagokat is beleértve), majd (rendre) hozzáadjuk a második, harmadik, illetve negyedik sorhoz:

 (4.2.3)

   

a11 0 0 0

a12 (2) a22 (2) a32 (2) a42

a13 (2) a23 (2) a33 (2) a43

a14 (2) a24 (2) a34 (2) a44

    

 b1   b(2)    2   =  (2)    b3  (2) b4 ai1 = b1 − a11 + bi , i ≥ 2. x1 x2 x3 x4





(2) (2) a ahol aij = a1j − i1 + aij , i, j ≥ 2, és bi a11 (2) Ha a22 6= 0, akkor a második lépésben a második sort szorozzuk be a

(2)

(2)

−

(4.2.4)

a32

(2) a22

=: e l32 , −

a42

(2) a22

=: e l42

értékekkel, majd hozzáadjuk a harmadik, illetve negyedik sorhoz:

 (4.2.5)

   

a11 a12 (2) 0 a22 0 0 0 0

a13 (2) a23 (3) a33 (3) a43

a14 (2) a24 (3) a34 (3) a44

    

x1 x2 x3 x4





    =  

b1 (2) b2 (3) b3 (3) b4

   . 

4.2. DIREKT MÓDSZEREK

56

Az eljárást mindaddig folytatjuk míg egy fels®háromszög lineáris egyenletrendszert kapunk:

    

(4.2.6)

a11 a12 a13 (2) (2) 0 a22 a23 (3) 0 0 a33 0 0 0

a14 (2) a24 (3) a34 (4) a44

    

x1 x2 x3 x4





    =  

b1 (2) b2 (3) b3 (4) b4

    

amit a (4.1.2) képletnek megfelel®en visszahelyettesítéssel oldunk meg.

68. Definíció. Az

Ha valamely

(2) a22

(i)

aii

(1)

(2)

(3)

(4)

a11 , a22 , a33 , a44

tagokat f®elemeknek nevezzük.

nulla - például a (4.2.3) egyenletrendszerben ha

= 0

- a séma nem alkalmazható, ugyanis az e l32 , e l42 szorzokat (2) (2) nem képezhetjük. Tételezzük fel, hogy a22 oszlopában (a22 alatt) (2) vannak nemnulla tagok. Ebben az esetben felcserélhetjük az a22 sorát a nemnulla tag sorával mert az egyenletrendszer megoldása nem módosul. Ezt az eljárást (részleges) f®elemkiválasztásnak nevezzük. Az új f®elem kiválasztásánál az els® nemnulla tag helyett az abszolút értékben o n ajánlatos (2) (2) legnagyobb tagot választani (max a32 , a42 a mi esetünkben), úgyanis így elkerülhet® az értékes számjegyek vesztése az osztás miatt. (2) (2) Ha a22 = 0 oszlopában (a22 alatt) minden elem nulla, akkor az

A

mátrix szinguláris, ugyanis

det A =

a11 a12 0 0 0 0 0 0

a13 (2) a23 (2) a33 (2) a43

a14 (2) a24 (2) a34 (2) a44

vagyis a lineáris egyenletrendszer határozatlan.

Ha a f®elemet a

0 a(2) a(2) 23 24 (2) (2) = a11 0 a33 a34 (2) 0 a(2) a44 43

b-t®l

= 0,

függ®en összeférhetetlen vagy

 (2) (2) (2) a22 a23 a24  (2) (2) (2)   a32 a33 a34  (2) (2) (2) a42 a43 a44 

részmátrixból választjuk ki,

akkor az eljárást teljes f®elemkiválasztásnak nevezzük. Ebben az esetben egy sor és egy oszlop cserét kell végrehajtani. Az oszlop cserét gyelembe kell venni az algoritmus végén a megoldások sorrendjének a leolvasásánál.

4.2. DIREKT MÓDSZEREK Az általános

57

n×n-es egyenletrendszer k+1-ik lépésben a kiküszöbölési

séma a következ®:

(k)

(4.2.7)

(k+1) aij

(4.2.8)

(k+1) bi

=

(k) aij

+

(k) akj

−

aik

(k)

=

+

(k) bk

, i, j ≥ k + 1,

akk (k)

(k) bi

!

−

aik

(k)

! , i ≥ k + 1.

akk

n×n-es lineáris egyenletrendszer Gauss-kiküszöbölési algoritmus 2 3 1 2 1 2 3 2 m¶veletigénye n − n − n = n + O (n ) . 3 2 6 3 Egy

Ha a kiküszöbölést nem csak az átló alatt hanem fölötte is elvégezzük akkor egy diagonális lineáris egyenletrendszerhez jutunk.

Ez a séma

Gauss-Jordan néven ismert. Habár a diagonális lineáris egyenletrendszert egyszer¶bb megoldani mint egy háromszög¶ egyenletrendszert (n összehasonlítva

n2 + O (n)

m¶velettel), a Gauss-Jordan kiküszöbölési m¶veletigénye

kétszer nagyobb a Gauss algoritmusénál, ezért összeségében a Gauss algoritmus gazdaságosabb mint a Gauss-Jordan. A Gauss algoritmus alkalmas determináns kiszámítására is ugyanis (4.2.9) ahol

(2)

det A = det A(1) = ... = det A(n) = a11 a22 ...a(n) nn

A(k) a k iterációban, Gauss módszerrel generált mátrix (A(1) = A).

69. Példa. Az

 (4.2.10)

    10 3 −1 x1 5       7 2 0   x2  =  5  −5 2 3 x3 −1

egyenletrendszer els® oszlopának a kiküszöbölésére az els® sorát beszorozzuk

5 e l31 = 10 -el, majd hozzáadjuk a második, illetve      10 3 −1 x1 5    15  −1 7  harmadik sorhoz =⇒  0 =  10 . A következ® x 10 10   2  25 35 15 0 10 10 x3 10 iterációban a második sort szorozzuk be e l = 35 -el a    23  és hozzáadjuk  10 3 −1 x1 5    15  −1 7  harmadik sorhoz =⇒  0 =  10  , majd a x 10 10   2  0 0 27 54 x3

7 e , l21 = − 10

illetve


58

fels®háromszög egyenletrendszert visszahelyettesítéssel megoldjuk:

x3 =

2, x2 = −1, x1 = 1. 

 10 3 −1  (3) 7  Az A mátrix determinánsát a (4.2.9) képlet szerint az A =  0 −1 10 10  0 0 27 (3) 27 = mátrix f®átlóján lév® tagok szorzata adja: det A = det A = 10 −1 10 −27.

A mátrixnak az inverzének gyelembe vesszük, hogy A ·

A Gauss módszert alkalmazhatjuk egy adott a meghatározására.

A

−1

= In

Ennek érdekében

tehát



1 0   0 1 A · (X1 |X2 | ... |Xn ) =    0 0

(4.2.11)

ahol

Xk

az

−1

A

inverz mátrixnak a

k -ik

oszlopa

 0  0   ..  .  1  x1k  2  xk Xk =   ..  . xnk

   .  

Oszlopokra bontva a (4.2.11)-b®l következik, hogy:



  1 0     0   1   A · X1 =   ..  , A · X2 =  ..  .   . 0 0





 0      0   , ..., A · Xn =  .  ,  .     .  1

Xk , k = 1, n oszlopainak a meghatározására n-szer alkalmazzuk a Gauss módszert egyforma A mátrixxal és különböz® vagyis az inverz mátrixnak az

szabadtag vektorral.


59

 10 3 −1   Az A =  7 2 0  mátrix inverzét A−1 = (X |Y | Z) −5 2 3 

70. Példa.

3 lineáris egyenletrendszer megoldásából kapjuk:



       10 3 −1 x1 1 x1 − 29          7 2 0   x2  =  0  =⇒ X =  x2  =  79 −5 2 3 x3 0 − 89 x3         11 0 y1 10 3 −1 y1 27   25        7 2 0   y2  =  1  =⇒ Y =  y2  =  − 27 35 0 y3 y3 −5 2 3 27         2 10 3 −1 z1 0 z1 − 27         7  7 2 0   z2  =  0  =⇒ Z =  z2  =  27 1 −5 2 3 z3 1 z3 27   11 2 − 29 27 − 27   −1 25 7 . tehát A =  79 − 27 27  35 1 − 89 27 27

4.2.2. LU faktorizáció

LU faktorizáción egy

  ,   ,   ,

A = (ai,j )i,j=1,n

mátrix

A = LU

(4.2.12)

szorzatra való felbontását értjük, ahol

L egy alsó- U

pedig egy fels®háromszög

mátrix:



l11 0 ... 0   l21 l22 ... 0 L=  ..  . ln1 ln2 ... lnn



u11 u12 ... u1n     0 u22 ... u2n , U =  .   .   . 0 0 ... unn

A két mátrix meghatározására a (4.2.12)-ban csak

U

mátrix f®átlóját

2

n 1-el

n2 + n

   .  

elem szükséges, ugyanakkor

egyenlettel rendelkezünk.

Ha az

L

vagy az

tesszük egyenl®vé a két mátrix egyértelm¶en

meghatározható. Az els® eset lii

uii = 1, i = 1, n



= 1, i = 1, n Doolittle, míg a második

Crout faktorizáció néven ismert.


60

A faktorizációs eljárásban a kiszámítási sorrendet úgy állítjuk fel, hogy az éppen sorra kerül® ismeretlen kiszámításánál minden adat rendelkezésünkre álljon. A Doolittle faktorizáció esetében alkalmazhatjuk

k -ik lépésben kiszámítjuk   ... 0 u11 u12 ... u1n     0 u22 ... u2n  ... 0  , U =  .    .    .  ... 1 0 0 ... unn

például a sor-szerinti rendezést, vagyis a



1 0   l21 1 L=  ..  . ln1 ln2

(4.2.13)

mátrixok

a



k -ik sorát lk1 , lk2 , ..., lk,k−1 , ukk , uk.k+1 , ..., uk,n (ebben a sorrendben)

lij =

(4.2.14)

(4.2.15)

aij −

Pi−1

k=1 lik ukj

ujj

uij = aij −

i−1 X

, i > j,

lik ukj , i ≤ j.

k=1



71. Példa. Az

 10 3 −1   A= 7 2 0  −5 2 3

mátrix esetében az

mátrixok tagjait a következ® sorrendben számítjuk ki:

 és a következ® mátrixokat kapjuk:

 L = 

1 7 10 5 − 10

L

és

U

u11 , u12 , u13 ; l21 , u22 , u23 ; l31 , l32 , u33 0 0  1 0 , U = −35 1



 10 3 −1  1 7  .  0 − 10 10  0 0 27 Ha a felbontás megtörtént, akkor az átírhatjuk

LU x = b

Ax = b

egyenletrendszert

alakra, ami ekvivalens az alábbi két háromszög

lineáris egyenletrendszerrel: (4.2.16)

Ly = b, U x = y.

E két lineáris egyenletrendszer (alsó-fels®) megoldásához használjuk a már ismert (4.1.4),(4.1.2) sémákat.

4.2. DIREKT MÓDSZEREK Az

LU

61

felbontás alkalmas a determináns kiszámítására is

det A = det (LU ) = det L · det U = det U =

(4.2.17)

n Y

uii .

i=1

A = LU

72. Példa. Ismerve az szerepl®

A

felbontását az (4.2.10) példában

mátrixnak, a (4.2.16)-ból kapjuk, hogy



1

 

7 10 5 − 10

    5 y1 0 0     1 0   y2  =  5  , −1 y3 −35 1

y = (5, 32 , 54)T , majd      10 3 −1 x1 5    3  1 7  = 2  x  0 − 10 10   2  0 0 27 54 x3

aminek megoldása

ahonnan

x = (1, −1, 2)T .

Ugyanakkor a (4.2.17) képletb®l kiszámítható

A determinánsa: det A =

det U = −27.

4.2.3. A Gauss és az LU módszer kapcsolata

Az

LU

faktorizációs

módszer levezethet® a Gauss módszerb®l, ugyanis ha a Gauss módszerben (k)

ismert szorzók

aik e lik = − (k) akk

segítségével megszerkesztjük az

L(k)

elemi

alsóháromszög mátrixot:



(4.2.18)

L(k)



1

 .  ..    0 =   0  .  ..  0

..

.

1 e lk+1,k

..

.

. . .

e lnk

1 0 1

          

akkor (4.2.19) ahol

A(k+1) = L(k) A(k) ,

A(k) a k iterációban, Gauss módszerrel generált mátrix (A(1) = A).


62

Tehát

U = A(n) = L(n−1) L(n−2) ...L(1) A

(4.2.20) ahonnan

A = L(1)

(4.2.21) Az

−1

... L(n−2)

−1

L(n−1)

−1

U.

L(k) mátrix egyszer¶ szerkezetének köszönhet®en könnyen kiszámítható

az inverze:



L(k)

(4.2.22)

−1

1  .  ..    0 =   0  .  ..  0

 ..

.

1 −e lk+1,k

..

.

. . .

1 0 1

−e lnk

          

és igazolható, hogy (4.2.23)



(1) −1

L

(n−2) −1

... L

L

(n−1) −1

     =    

1 −e l21 1 −e l31 −e l32 e e −l41 −l42

 ..

.

. . .

−e ln1 −e ln2 · · · vagyis (4.2.21)-b®l

1 −e ln,n−1 1

     =L    

A = LU .



10  (1) 73. Példa. Az el®bbi példában A = A =  7 −5 7 5 Gauss a szorzók rendre e l21 = − 10 ,e l31 = 10 , majd e l32

 3 −1  2 0  2 3 = 35

a

voltak

4.3. LINEÁRIS EGYENLETRENDSZEREK ITERATÍV MEGOLDÁSA

63

=⇒ 

   1 0 0 10 3 −1  7   1 7  L(1) =  − 10 1 0  =⇒ A(2) = L(1) A(1) =  0 − 10 10  5 5 7 0 1 0 10 2 2     1 0 0 10 3 −1    7  1 L(2) =  0 1 0  =⇒ A(3) = L(2) A(2) =  0 − 10 =U 10  0 0 27 0 35 1      1 0 0 1 0 0 1 0 0 −1 (2) −1  7    7  L = L(1) L =  10 1 0   0 1 0  =  10 1 0 . 5 5 − 10 0 1 0 −35 1 − 10 −35 1

4.3. Lineáris egyenletrendszerek iteratív megoldása Az iteratív eljárás a numerikus számítások egyik alapvet® módszerének számít. Az

(Xn )n

megoldások sorozatát fokozatosan közelítjük a:

X (k+1) = Φ X (k)

(4.3.1)

iterációkat használva, vagy

X (k+1) = Φ X (k) , X (k−1) , ...

ha az eljárás

többlépéses. A direkt módszerekkel ellentétben az iteratív módszerekben az iterációk száma nem a mátrix méretét®l függ, hanem egy adott

> 0 pontosságtól.

?Míg a direkt módszerekben minden iterációval a hiba kumulálodott, addig az iteratív eljárásban minden iterációval közelebb kerülünk a keresett megoldáshoz.?

4.3.1. Jacobi-féle módszer

Ahhoz, hogy az

AX = b

lineáris

egyenletrendszert (4.3.1) alakra hozzuk kifejezzük az egyenletrendszer minden

i-edik

sorában az

xi

ismeretlent:

  a11 x1 + a12 x2 + ... + a1n xn = b1     a21 x1 + a22 x2 + ... + a2n xn = b2 ⇒ . .  .     a x + a x + ... + a x = b n1 1 n2 2 nn n n


64

 12  x1 = − aa11 x2 · · · − aa1n xn + ab111  11    x2 = − a21 x1 . . . − aa2n xn + ab222 a22 22 . .  .     x = − an1 x − an2 x . . . n + abnn n ann 1 ann 2 vagy mátrix alakban:

X = BX + C,

(4.3.2) ahol (4.3.3)



  . . . − aa1n 0 − aa12 x1 11 11    a21  x2   − a22 0 . . . − aa2n 22   X=  ..  , B =  ..  .   . n2 n1 xn − aann ... 0 − aann Φ-vel

    , C =     

  .  

X (k+1) = Φ X (k) , k ≥ 0

X (0)

egy kezdeti megoldás.

A Banach-féle xpont tételnek megfelel®en ha

X

bn ann



iteratív eljáráshoz jutunk, vagyis

(4.3.5)

(k)

b1 a11 b2 a22 . . .

Φ (X) = BX + C,

X = Φ (X)

ahol



jelölve a (4.3.2) képlet jobboldalát

(4.3.4) az



k

megoldások sorozata konvergál tetsz®leges

X

Φ =kontrakció (0)

az

kezdeti megoldásra.

Felhasználva a (4.3.4) képletet következik, hogy

||Φ (X − Y )|| = ||BX + C − BY − C|| = ||B (X − Y )|| ≤ ||B|| ||X − Y || , tehát ha

||B|| < 1

akkor

Φ

kontrakció és az eljárás konvergál.

Ha ||B|| < 1 akkor a (4.3.5) fokozatos közelítések sorozata konvergens. 74. Tétel.

A tételben szerepl® elégséges feltétel átírható az eredeti Inf típusú normát használva következik, hogy:

||B||∞

n X aij < 1, i = 1, n = max aii i,j=1,j6=i

A mátrixra.


65

ami ekvivalens az alábbi feltétellel: (4.3.6)

n X

|aij | < |aii | , i = 1, n.

j=1,j6=i


A mátrixok átlósan dominánsnak nevezünk ha

eleget tesz az (4.3.6) kikötésnek.



76. Példa. Az

 −7 2 −2   A =  3 5 −1  0 2 −3

mátrix átlósan domináns.

Ha az Ax = b egyenletrendszer A mátrixa átlósan domináns, akkor az (4.3.5)-ban megadott iterációk konvergálnak a lineáris egyenletrendszer megoldásához. 77. Következmény.

   −7x1 + 2x2 − 2x3 = −1 78. Példa. 3x1 + 5x2 − x3 = 14 , AX = b ⇔ X =   2x2 − 3x3 = 7     1 0 27 −2 7 7  14   1  , C = BX + C ahol B =  −3 0  5  . Ha X (0) = 5 5  −7 0 23 0 3   0    1  egy tetsz®leges kezdeti megoldás, akkor X (1) = BX (0) + C = 1     1 0.14 7      3  '  3  és az abszólut eltérés a két megoldás között − 35 −1.67     31 1.48 21 (1)  50   X − X (0) = 8 ' 2.67. Hasonlóan X (2) = BX (1) +C =  ' 2.38 ,  3 21   ∞ − 31 −0.33     45 0.92 49  194    (3) (2) X = BX +C =  105  '  1.85  és a hibák X (2) − X (1) ∞ = − 47 −0.75 63 4 82 (3) (2) ' 1.33, X − X = 147 ' 0.55. Mivel az A mátrix átlósan 3 ∞

4.3. LINEÁRIS EGYENLETRENDSZEREK ITERATÍV MEGOLDÁSA domináns az

limk→∞ X (k)

A

X (k) 

k

66

megoldások sorozata konvergál a pontos megoldáshoz:

 1   =  2 . −1

(k + 1)-ik iterációban az X

(k+1)

=

(k+1) x1

(k+1) x2

...

(k+1) xn

t

megoldás komponenseit az alábbi képlettel számítjuk ki:

(k+1) (4.3.7)xi

(4.3.8)

1 (k) (k) (k) = bi − ai1 x1 − ai2 x2 − . . . / . . . − ain xn aii ! X 1 (k) = bi − aij xj , i = 1, n. aii j=1,j6=i

A (4.3.6) egy elégséges de nem szükséges feltétel, vagyis átlósan nem-domináns lineáris egyenletrendszer megoldása konvergálhat a Jacobi módszerrel.

Egy szükséges és elégséges feltétel a konvergenciára a

spektrál sugár segítségével fejezhet® ki.


A ∈ Rn×n mátrixnak a sajátértékeinek a halmazát

spektrumnak nevezzük

σ (A) = λi : λi = sajátérték, i = 1, n . Spektrálsugárnak nevezzük az

A

mátrix legnagyobb sajátértékének az

abszolút értékét

ρ (A) = |λ1 | , ahol

|λ1 | ≥ |λ2 | ≥ . . . ≥ |λn | . (λ, x)

Ha

λk , x

vagyis

páros sajátérték-sajátvektora az

páros sajátérték-sajátvektora az k k

ρ A

Ak

A

mátrixnak, akkor a

mátrixnak

k = 1, 2, ...,

= (ρ (A)) .

80. Tétel.

(i) Bármilyen ||·|| szubmultiplikatív mátrix norma esetében

igaz az alábbi egyenl®tlenség:

ρ (A) ≤ ||A|| .


67

(ii) Bármilyen > 0-ra létezik egy ||·|| szubmultiplikatív mátrix norma

amelyre

ρ (A) ≤ ||A|| ≤ ρ (A) + . (i) Ha (λ, x) páros sajátérték-sajátvektora az A mátrixnak, akkor |λ| ||x|| = ||λx|| = ||Ax|| ≤ ||A|| ||x||, és mivel x 6= 0 következik, hogy |λ| ≤ ||A|| . Bizonyítás.

A (4.3.5) sorozat konvergál (bármilyen X (0) kezdeti értékre), akkor és csak akkor ha 81. Tétel.

ρ (B) < 1.

(4.3.9)

Bizonyítás. Feltételezzük, hogy

ρ (B) < 1.

Akkor létezik

> 0

ρ (B) + < 1. Az el®bbi tétel szerint létezik egy norma amelyre ||B|| ≤ ρ (B) + < 1, tehát Φ kontrakció és a sorozat konvergens. (k+1) Feltételezzük, hogy X = Φ X (k) = BX (k) + C sorozat konvergál X ∗ -hez bármilyen kezdeti megoldásra. Akkor legyen X (0) ú.h. X (0) − X ∗ sajátvektora B -nek. Következik, hogy X (k+1) −X ∗ = Φ X (k) −Φ (X ∗ ) = B X (k) − X ∗ = ... = B k+1 X (0) − X ∗ = λk+1 X (0) − X ∗

ú.h.

k → ∞ esetén zéróhoz közelít, ugyanez igaz a jobboldalra, tehát λ → 0 ami akkor igaz ha |λ| < 1. Mivel tetszöleges λ-ra igaz következik, hogy ρ (B) < 1. ( x1 − 2x2 = −1 82. Példa. , AX = b ⇔ X = BX + C −9x1 + 32x2 = 23 ! ! −1 0 2 , C = . Az A mátrix nem átlósan ahol B = 23 9 0 32 32 domináns és a B mátrix normái mind nagyobbak 1-nél, tehát a konvergenciára vonatkozó elégséges feltétel nem használható. A B mátrix sajátértékei −λ 2 9 det (B − λI) = 0 ⇔ 9 = 0, ahonnan = 0 ⇔ λ2 − 16 32 −λ |λ1,2 | = 43
k+1


−1

!

(1) X − X (0) = 1.23, X (2) = BX (1) + C = , C= 23 2 32 ! (2) 1 X − X (1) = 1.46..., X (k) −−−→ X ∗ = . Az alábbi k→∞ 2 1

68 7 16 7 16

! ,

ábrán az

els® 20 iteráció hibaváltozása látható. ÁBRA 'Lejacobihiba.eps'

4.3.2. Gauss-Seidel és SOR módszerek

A Jacobi módszer (4.3.7)

(k + 1)-ik iterációban a megoldás csak az el®bbi, (k)-ik (k+1) komponensekt®l függnek, viszont a xi komponens kiszámításánál a (k+1) (k+1) , ..., xi−1 értékek már ismertek és jobb közelítést adnak a pontos x1 képlet szerint a

megoldásnak. Tehát a

(k+1) xi

1 (k) (k+1) (k+1) (k+1) (k) = bi − ai1 x1 − ai2 x2 . . . − ai,i−1 xi−1 / − ai,i+1 xi+1 . . . − ain xn , i = 1, n aii

vagy

(4.3.10)

1 = aii

(k+1) xi

! bi −

X

(k+1) aij xj

−

j
X

(k) aij xj

, i = 1, n

j>i

iteratív eljárás gyorsabb mint a (4.3.7). A (4.3.10) algoritmus GaussSeidel néven ismert. A Jacobi és a Gauss-Seidel módszerek mátrix alakját levezethetjük ha az

A

mátrixot felbontjuk az alábbi módon:

A=L+D+U

(4.3.11) ahol

L,

illetve

U

szigorúan alsó, illetve fels® mátrixok



  0 ··· 0 0 0 a12 . . . a1n     a21 0  0 0 0  a2n ,U =  . L= .. .  ..   .. .  .   .. an1 . . . an,n−1 0 0 ... 0 0 és

D

átlós mátrix



a11 · · ·   0 a22 D=  ..  . 0 ...

0 ..

0

.

 0  0  .   ann

     

4.3. LINEÁRIS EGYENLETRENDSZEREK ITERATÍV MEGOLDÁSA Az

Ax = b

egyenlet átírható

(L + D + U ) x = b

69

alakba vagyis

Dx = b − (L + U ) x ahonnan

x = D−1 (b − (L + U ) x) és az iteratív eljárás: (4.3.12) A

x(k+1) = D−1 b − (L + U ) x(k) .

D egyszer¶ alakjának köszönhet®en a D−1 mátrixot könny¶ kiszámítani:   1 · · · 0 0  a11 1   0 a22 0  −1   D = . ..  . .  .  1 0 . . . 0 ann

A (4.3.12) képlet azonos a Jacobi (4.3.7) algoritmussal. A Gauss-Seidel módszerhez a csoportosítást a következ®képpen végezzük el

Ax = b ⇔ ((L + D) + U ) x = b ⇔ (L + D) x = b − U x x = (L + D)−1 (b − U x) és az iteratív eljárás (4.3.13)

x(k+1) = (L + D)−1 b − U x(k)

mely azonos a (4.3.10) eljárással. SOR succ over relax

ω ∈ (0, 2) ωAx = ωb ⇔ ω (L + D + U ) x = ωb ⇔

(ωL + D + (1 − ω) D + ωU ) x = ωb

(k+1)

xi

x = (ωL + D)−1 (ωb − (ωU + (1 − ω) D) x) x(k+1) = (ωL + D)−1 ωb − (ωU + (1 − ω) D) x(k) ! X X ω (k+1) (k) (k) bi − aij xj − aij xj , i = 1, n = (1 − ω) xi + aii ji

4.4. HIBABECSLÉS, KONDÍCIONÁLÁS

4.3.3. Rosszul kondicionált matrixok

1.

70 Golub, van Loan-

pg.107

Ax = b lineáris egyenletrendszerek pontosságában szerepet játszik az A mátrix, mint a b szabadtag.

Az úgy

!

83. Példa.

A=

míg a determinánsa

1 0 mátrix 0 10−6 det (A) = 10−6 .

kondiciószáma

K (A) = 106 ,

4.4. Hibabecslés, kondícionálás Az

Ax = b lineáris egyenletrendszer úgy A mint a b tagjai meggyelések

vagy számítások eredménye, tehát mindkét esetben az eredeti értékeknek csak a közelítése áll rendelkezésünkre. befolyásolja a megoldást az vektor

∂b

perturbálása az

x

A

ill

b

Felmerül a kérdés mennyire

közelit®. Feltételezzük, hogy az

megoldás egy

∂x

változását erdményezi.

Ekkor:

A (x + ∂x) = b + ∂b, =⇒ Ax + A∂x = b + ∂b, majd egyszer¶sítés után

A∂x = ∂b, vagyis

∂x = A−1 ∂b. Tehát (4.4.1)

||A|| ||A−1 || ||∂b|| ||∂x|| = A−1 ∂b ≤ A−1 ||∂b|| = , ||A||

ahonnan a relatív eltérés

||∂b|| ||∂x|| = ||A|| A−1 . ||x|| ||A|| ||x|| Felhasználva a

| |b| | = ||Ax|| ≤ ||A|| ||x|| egyenl®tlenséget következik, hogy (4.4.2)

b

||∂b|| ||∂x|| ≤ ||A|| A−1 , ||x|| ||b||

4.5. TÚLHATÁROZOTT LINEÁRIS EGYENLETRENDSZEREK vagyis az

71

x relatív hibája kisebb mint a b vektor relatív hibája szorozva

egy számmal. 84. Definíció. A

K (a) = ||A|| ||A−1 || számot az A mátrix kondíciószámának

nevezzük. Hasonló a helyzet ha az

A

mátrix helyett egy perturbált

A + ∂A

mátrixot használunk:

(A + ∂A) (x + ∂x) = b, =⇒ Ax + A∂x + ∂A (x + ∂x) = b =⇒ A∂x = −∂A (x + ∂x) =⇒ ∂x = −A−1 ∂A (x + ∂x) =⇒ ||∂x|| = A−1 ∂A (x + ∂x) ≤ A−1 ||∂A|| ||x + ∂x|| =⇒ ||∂x|| ||∂A|| ≤ A−1 ||∂A|| = A−1 ||A|| ||x + ∂x|| ||A|| ahonnan

||∂A|| ||∂x|| ≤ K (A) . ||x + ∂x|| ||A||

(4.4.3)

Habár a (4.4.3) képletben csak az

x relatív hibának a közelítése szerepel,

a kondíciószám jelent®sége jól szemléltethet®. Tehát a lineáris egyenletrendszer megoldásának relatív hibáját majorálja a kondíciószám szorzata kicsik. Ha

K (A)

A vagy b relativ hibájával amelyek rendszerint

is kis érték akkor a szorzat ebb®l kifolyolag a relatív

hiba is az. Ebben az esetben stabil vagy jól kondíciónált rendszerekr®l beszélünk. Viszont ha

K (A)

nagy, akkor a szorzat is az (lehet), ami

nagy relatív hibához vezethet. A kondiciószám legkisebb értéke 1:

1 = ||I|| = A A−1 ≤ ||A|| A−1 = K (A) . Habár a kondiciószám függ a normától, ennek nagyságrendje a fontos ezért mindegy milyen normával dolgozunk.

4.5. Túlhatározott lineáris egyenletrendszerek A gyakorlatban el®fordul, hogy olyan lineáris egyenletrendszereket kell megoldani aminek több (vagy kevesebb) egyenlete van mint ismeretlenje. Például ha egy bizonyos

n

számú paraméter meghatározására

n-nél

4.5. TÚLHATÁROZOTT LINEÁRIS EGYENLETRENDSZEREK

72

több kisérletet végzünk a keletkez® egyenletrendszer nem lesz négyzetes. Gyakran válik ilyen esetben az egyenletrendszer összeférhetetlené, annak ellenére hogy a paraméterek léteznek. Ennek egyik oka lehet például a kisérleteknek hibás végzése, leolvasása. Egy

Ax = b

(4.5.1)

A = (aij )i=1,m , x = (x1 x2 ... xn )t , b = (b1 b2 ... bm )t

(4.5.2)

j=1,n egyenletrendszer®l azt mondjuk hogy túlhatározott ha

85. Példa.

   x1 + x2 = 2 x1 + 2x2 = 3   2x1 + x2 = 4

n > m.

Az egyenletrendszernek nincs megoldása.

A feladat megoldása érdekében vezessük be a következ®

r ∈ Rm

(maradék) vektort:

r = Ax − b

(4.5.3)

Ekkor a (4.5.1) egyenletrendszernek akkor van (klasszikus értelemben vett) megoldása ha

r = θRm . Ha r 6= θRm

akkor az adott egyenletrendszer

összeférhetetlen és ebben az esetben keressük a megoldást a legkisebb négyzetek módszere szerint vagyis, az hogy az

r

x

vektort úgy határozzuk meg

vektor normája minimális legyen:

(4.5.4)

||r|| = ||Ax − b|| −→ min .

86. Definíció. Azt az

x vektort ami eleget tesz a (4.5.4) kikötésnek

az egyenletrendszer legkisebb négyzetek szerinti megoldásának nevezzük. Természetesen ha min||r||

= 0

akkor visszakapjuk a klasszikus

értelemben vett megoldást. Ha a 2-es típusú normát használjuk akkor a (4.5.4) kikötés ekvivalens a következ®vel:

 2 !2 n m n X X X (||Ax − b||2 )2 =  aij xj − bi  = aij xj − bi −→ min . j=1

2

i=1

j=1

4.5. TÚLHATÁROZOTT LINEÁRIS EGYENLETRENDSZEREK Jelöljük

f -el

73

f : Rn → R !2 m n X X f (x) = aij xj − bi

a következ® függvényt

i=1 Ahhoz hogy az

j=1

f -nek a minimumát meghatározzuk kiszámítjuk a stacionárius

pontjait:

m n X X ∂f = 2 aij xj − bi ∂xk i=1 j=1

⇔

m X

aik

i=1

n X

! aik = 0, k = 1, ..., n

! aij xj − bi

= 0, k = 1, ..., n.

j=1

ami átírva mátrix alakra:

At (Ax − b) = θRm vagy:

At Ax = At b.

(4.5.5)

Tehát, az (4.5.1) egyenletrendszer legkisebb négyzetek szerinti megoldása megegyezik a (4.5.5) egyenletrendszer klasszikus értelemben vett megoldásával.

87. Tétel.

Az (4.5.5) lineáris egyenletrendszer megoldására lesz az ||r||2 = ||Ax − b||2

kifejezés minimális. Bizonyítás. Jelöljük

rx = b − Ax ry = b − Ay ⇒ ||ry || ≥ ||rx || . ry = b − Ay = b − Ax + Ax − y = rx + A (x − y) ryt = rxt + (x − y)t At

4.5. TÚLHATÁROZOTT LINEÁRIS EGYENLETRENDSZEREK

ryt ry = rxt + (x − y)t At (rx + A (x − y)) =

(4.5.6) (4.5.7)

74

rxt rx + (x − y)t At rx + rxt A (x − y) + (x − y)t At A (x − y) = = rxt rx + (x − y)t At A (x − y)

(4.5.8) mert

At rx = At (b − Ax) = θ es

rxt A = At rx

t

= θ.

Mivel

ryt ry = ||ry ||22 a (4.5.8)-ból következik hogy:

||ry ||22 = ||rx ||22 + ||A (x − y)||22 ≥ ||rx ||22 . 88. Példa. Az el®bbi példát felhasználva, keressük az (4.5.5) egyenletrendszer megoldását:

1 1 2 1 2 1

!



 1 1    1 2  2 1

⇔ 6 5 5 6

!

⇔ x1 x2 Tehát

x1 = 18/11, x2 = 7/11

x1 x2

!

x1 x2

!

! =

=

1 1 2 1 2 1

=

13 12

18 11 7 11

!



 2    3  4

!

!

az egyenletrendszer legkisebb négyzetek

módszerének megfelel® megoldása.

5. FEJEZET

Sajátérték, sajátvektor numerikus kiszámítása Legyen

A ∈ Mn (R)


egy négyzetes mátrix.

X ∈ Cn , X 6= 0 vektort az A mátrix sajátvektorának

nevezzük ha (5.0.9)

AX = λX.

λ skalárt az A mátrix, X vektornak megfelel® sajátértéknek nevezzük. ! ! ! 3 1 1 1 90. Példa. Ha A = , akkor A =4 , tehát 2 2 1 1 ! 1 X= a mátrix sajátvektora és λ = 4 a megfelel® sajátértéke. 1 A

A (5.0.9) képletb®l k®vetkezik, hogy

AX − λX = 0 (5.0.10)

(A − λIn ) X = 0.

A (5.0.10) homogén egyenletrendszernek akkor van nullától kül®nböz® megoldása ha (5.0.11) 91. Definíció. A

det (A − λIn ) = 0. det (A − λIn )

karakterisztikus determinánsnak,

az (5.0.11) egyenletet pedig karakterisztikus egyenletnek nevezzük. A (5.0.11) karakterisztikus determináns egy

n-ed fokú polinom λ-ba

PA (λ) := det (A − λIn ) = (−1)n λn − S1 λn−1 + S2 λn−2 − ... + (−1)n Sn Si a f®átlón elhelyezked® i×i minorok összege: S1 = T r (A) , ..., Sn = det (A) . ahol

75

5.1. KRYLOV MÓDSZER

76

A Viéte összefüggésekb®l következik, hogy

n Y

λi = λ1 ...λn = det (A) ,

i=1

A szinguláris det (A) = 0, ∃k ∈ 1, n, ú.h. λk = 0.

vagyis ha zéró:

92. Tétel.

akkor (legalább) egy sajátérték

Az A mátrix kielégíti a saját karakterisztikus egyenletét An − S1 An−1 + S2 An−2 − ... + (−1)n Sn In = 0.

(5.0.12)

5.1. Krylov módszer 93. Definíció. Az

A ∈ Mn és y = (y1 , y2 , ..., yn )T ∈ Rn oszlopmátrix

által generált

Kn (A, y) = span y, Ay, ..., An−1 y

(5.1.1)

Krylov vektorrendszernek nevezzük.

A karakterisztikus polinom

Si , i = 1, n

együtthatói kiszámítása

érdekében átírjuk a (5.0.12) egyenletet

An + S1 An−1 + S2 An−2 + ... + Sn In = 0,

(5.1.2)

alakba, majd jobboldalról beszorozzuk egy (5.1.3)

y (0) ∈ Rn

vektorral

⇒

An y (0) + S1 An−1 y (0) + ... + Sn−1 Ay (0) + Sn In y (0) = 0.

Jelöljük (5.1.4)

y (1) := Ay (0) , y (2) := A2 y (0) = Ay (1) , ..., y (n) := An y (0) = Ay (n−1) , amit behelyettesítve a (5.1.3) a következ® lineáris egyenletrendszert eredményezi:

y (n) + S1 y (n−1) + ... + Sn−1 y (1) + Sn y (0) = 0, vagyis (5.1.5)

S1 y (n−1) + ... + Sn−1 y (1) + Sn y (0) = −y (n) .

5.2. HATVÁNY MÓDSZER

77

Ha a (5.1.1) vektorrendszer lineárisan független, akkor a (5.1.5) négyzetes lineáris egyenletrendszer összeférhet®, és a megoldását behelyettesítjük a

λn + S1 λn−1 + S2 λn−2 − ... + Sn = 0 karakterisztikus egyenletbe, majd megoldjuk.

94. Példa. Az

A=

3 1 2 2

! mátrix karakterisztikus egyenletét

λ2 + S1 λ + S2 = 0 alakba, ahol S1 , S2 együtthatók az S1 y (1) + S2 y (0) = −y (2) egyenletrendszer megoldásai. !Legyen y (0) ∈ R2 egy 3 (0) . A (5.1.4)-b®l következik, tetsz®leges oszlopvektor, például y = −2 ! ! 7 23 (1) hogy y = , y (2) = , tehát a (5.1.5) lineáris egyenletrendszer 2 18 ( 7S1 + 3S2 = −23 , melynek megoldása S1 = −5, S2 = 4. Innen 2S1 − 2S2 = −18 2 a karakterisztikus egyenlet λ − 5λ + 4 = 0 melynek megoldása λ1 = 4, λ2 = 1. ! ! −1 −1 (0) (1) Ha y = akkor y = y (2) = y (0) = tehát az 2 2 felírjuk

(5.1.1) vektorrendszer lineárisan összefügg® és az egyenletrendszernek nincs megoldása. Ebben az esetben más

y (0)

választunk.

5.2. Hatvány módszer Feltételezzük, hogy az

A ∈ Mn (R)

mátrix

λi , i = 1, n

sajátértékei

valósak és különböz®ek:

|λ1 | > |λ2 | > ... > |λn | . A hatvány módszer segítségével meghatározható a legnagyobb (domináns) sajátérték

λ1 ,

illetve a megfelel® (domináns) sajátvektor

X1.

6. FEJEZET

Függvény approximáció, interpoláció 6.1. Analitikus függvények numerikus kiszámítása Taylor sorral Legyen

f : [a, b] → R

egy végtelenszer deriválható függvény, és

x0 ∈ (a, b).

95. Definíció. Az

f

függvényt analitikusnak nevezzük ha hatványsorként

el®állítható: (6.1.1)

1 1 1 f (x) = f (x0 )+ f 0 (x0 ) (x − x0 )+ f 00 (x0 ) (x − x0 )2 +...+ f (n) (x0 ) (x − x0 )n +... 1! 2! n!

Tehát az megfelel®

Rn

f

függvény felírható mint

Tn

egy

n-ed

fokú polinom, és a

maradék összegeként:

f (x) = Tn (x) + Rn (x) ahol

1 1 1 Tn (x) = f (x0 )+ f 0 (x0 ) (x − x0 )+ f 00 (x0 ) (x − x0 )2 +...+ f (n) (x0 ) (x − x0 )n 1! 2! n! és a maradék tag :

1 n+1 (n+1) Rn (x) = f (x0 ) (x − x0 ) + ... . (n + 1)! Numerikusan az

f

függvényt az els®

n+1

taggal, vagyis az

Taylor polinommal közelítjük meg:

n X 1 (i) f (x) ' Tn (x) = f (x0 ) (x − x0 )i . i! i=0

78

n-ed

fokú

6.2. INTERPOLÁCIÓ

79

Az Rn marék tagot a következ® ún. Lagrange -féle alakban lehet kifejezni: 96. Tétel.

1 n+1 (n+1) Rn (x) = f (θ) (x − x0 ) , θ ∈ (x, x0 ) . (n + 1)! Az

x0 = 0 pont körüli sorfejtés egy sajátos esete a Taylor kifejtésnek

és a MacLaurin nevet viseli. 97. Példa.

1 2 x 2!

+ ... +

ex = 1 + 1!1 x + 2!1 x2 + ... + n!1 xn + ... =⇒ ex = 1 + 1!1 x +

1 n x és a hiba: n!

eθ n+1 Rn (x) = θ , θ ∈ (0, x0 ) (n + 1)! Innen ha

x = 1/2 =⇒

(vagy

θ ∈ (x0 , 0) ).

√ 2 n e = 1 + 0.5 + 0.5 + ... + 0.5 1! 2! n!

és

n-et a következ®

feltételb®l számítjuk ki:

eθ θn+1 < , θ ∈ (0, 1/2) Rn (x) = (n + 1)! ahol

>0

egy el®re megadott pontosság. Felhasználva hogy

e ∈ (2, 3)

a következ® egyenl®tlenséghez jutunk:

√ 3θ 3 1 1 n+1 (n + 1)! θ < (n + 1)! 2 < (n + 1)! < . A Taylor polinom csak lokálisan, az

x0

pont környékén alkalmas a

közelítésre.

6.2. Interpoláció Gyakran történik meg a gyakorlatban hogy bizonyos folyamatokat (csak) diszkrét mérésekkel lehet elemezni.

x1 < ... < xn

Például különböz®

x0 <

id®pontokban (nevezzük interpolációs alappontoknak

vagy osztópontok) a mérések eredményei mint a folyamatot leíró

f

y0 , y1 , ..., yn

amik nem másak

(ismeretlen) függvény értékei az alappontokban:

y0 = f (x0 ) , y1 = f (x1 ) , ..., yn = f (xn ) . F : [x0 , xn ] → alappontokban f -el:

Interpolációnak nevezzük azt az eljárást amellyel egy olyan

R

függvényt értelmezünk amely megegyezik az

(6.2.1)

F (xi ) = f (xi ) = yi , i = 0, ..., n.

6.2. INTERPOLÁCIÓ

80

Mértanilag ez azt jelenti hogy olyan meghatározott típusú

y = F (x)

(x0 , y0 ) , ..., (xn , yn )

görbét kell keresnünk, amely az adott

pontokra

illeszkedik. ÁBRA akármi interp Kézenfekv® az interpoláló

F

függvényt polinom (algebrai vagy trigonometrikus),

racionális függvény, stb. alakban keresni. Extrapoláció azt az eljárást jelenti amellyel az kívül próbáljuk az

f -et

megközelíteni.

6.2.1. Polinomiális interpoláció ... < xn

[x0 , xn ] intervallumon

Feltételezzük hogy az

x0 < x1 <

alappontok esetében ismertek a következ® adatok:

yi = f (xi ) , i = 0, ..., n. Az interpoláló

F

függvényt (algebrai) polinom alakban keressük:

F (x) = an xn + ... + a1 x + a0 ; mivel az ismeretlenek száma = egyenl® mint

n+1

a polinom fokszáma kisebb vagy

n különben megeshet hogy nem létezik a keresett polinom.

Ha az interpoláló polinom fokszáma egyenl® n-el akkor F egyértelm¶en meghatározott. 98. Tétel.

F interpoláló polinom: F (xi ) = yi , i = 0, ..., n =⇒   an xn0 + ... + a1 x0 + a0 = y0     an xn + ... + a1 x1 + a0 = y1 1 . . .  .     a xn + ... + a x + a = y

Bizonyítás. Mivel

n n

1 n

0

n

Az egyenletrendszer determinánsa Vandermonde típusú:

xn0 xn−1 ... x0 1 0 Y xn1 xn−1 ... x1 1 1 (xj − xi ) 6= 0, = ... 0≤xi <xj ≤n xnn xn−1 ... xn 1 n

6.2. INTERPOLÁCIÓ és mivel

xi 6= xj , ∀i 6= j,

81

a determináns különbözik nullától, tehát a

Cramer tétel szerint az egyenletrendszer egyértelm¶en meghatározott.

A fenti egyenletrendszer gyakran lesz rosszul kondicionált, ezért az

F

interpoláló függvény meghatározására más módszert használunk.

99. Példa. Az

 mátrix

  A= 

xi yi

0

1

4

9

0

1

2

3

0 0 0 1 1 1 64 16 4 729 81 9

1 1 1 1

interpoláló pontok esetében a Vandermonde

    

kondíciószáma

6.2.1.1. Lagrange -féle interpolációs polinom pontokat összeköt®

n-ed

cond (A) = 1.66 ∗ 103 .

(x0 , y0 ) , ..., (xn , yn ) Ln f -el. Akkor az

Az

fokú polinomot jelöljük

interpolációs feltételek a következ®k: (6.2.2) Az

Ln f

Ln f (xi ) = f (xi ) = yi , i = 0, ..., n. polinomot a következ® alakban adjuk meg:

Ln f (x) = y0 l0 (x) + y1 l1 (x) + ... + yn ln (x) ahol li

(x) , i = 0, ..., n az úgynevezett n-ed fokú Lagrange -féle alappolinomok.

Az interpoláló (6.2.2) tulajdonságot felhasználva kapjuk a következ® egyenletrendszert:

  y0 l0 (x0 ) + y1 l1 (x0 ) + ... + yn ln (x0 ) = y0     y0 l0 (x1 ) + y1 l1 (x1 ) + ... + yn ln (x1 ) = y1 . .  .     y l (x ) + y l (x ) + ... + y l (x ) = y 0 0

n

1 1

n

n n

n

aminek egyetlen megoldása: (6.2.3)

l0 (x0 ) = 1, l1 (x0 ) = 0, ..., ln (x0 ) = 0

(6.2.4)

l0 (x1 ) = 0, l1 (x1 ) = 1, ..., ln (x1 ) = 0

(6.2.5) (6.2.6)

... l0 (xn ) = 0, l1 (xn ) = 0, ..., ln (xn ) = 1

n

6.2. INTERPOLÁCIÓ vagy a

δij

82

Kronecker szimbólumot felhasználva:

( li (xj ) = δij =

1, 0,

i=j , i, j = 0, ..., n. i 6= j

ha ha

li (x) polinomnak n x0 , ..., xi−1 , xi+1 , ..., xn , tehát:

A (6.2.3)-ból következik hogy minden gyöke van, nevezetesen:

különböz®

li (x) = ci (x − x0 ) ... (x − xi−1 ) (x − xi+1 ) ... (x − xn ) . A

ci

együttható meghatározható a li

ci =

(xi ) = 1

egyenletb®l:

1 , i = 0, ..., n (xi − x0 ) ... (xi − xi−1 ) (xi − xi+1 ) ... (xi − xn )

vagyis:

ci = Q

1 , i = 0, ..., n. k6=i (xi − xk )

Ennek segítségével felírjuk az alappolinomokat:

Q k6=i (x − xk ) li (x) = Q , k6=i (xi − xk )

(6.2.7)

vagy bevezetve a következ® polinomot: aminek deriváltja

$ (x) = (x − x0 ) ... (x − xn ) ,

xi -ben:

$0 (xi ) = (xi − x0 ) ... (xi − xi−1 ) (xi − xi+1 ) ... (xi − xn ) , a (6.2.7) képlet a következ® lesz:

li (x) =

$ (x) . (x − xi ) $0 (xi )

Tehát a Lagrange interpoláló polinom a következ® alakot veszi fel:

(6.2.8)

Ln f (x) =

n X

yi · li (x) =

i=0

n X i=0

yi ·

$ (x) . (x − xi ) $0 (xi )

A (6.2.8) Lagrange interpolációs polinommal való közelítésb®l ered® hiba a következ®: 100. Tétel.

(6.2.9)

(n+1) f (c) |f (x) − Ln f (x)| = $ (x) , c ∈ (a, b) . (n + 1)! n = 1 és n = 2. Az els® esetben (x0 , y0 ) , (x1 , y1 ) pontokat összeköt®

Sajátos esetként megemlíthet® az

L1 f (x)

polinom nem más mind a

6.2. INTERPOLÁCIÓ

83

egyenes egyenlete:

L1 f (x) = y0 míg

n = 2-re

az

x − x1 x − x0 + y1 , x ∈ [x0 , x1 ] , x 0 − x1 x1 − x0

(x0 , y0 ) , (x1 , y1 ) , (x2 , y2 )

pontokat összeköt® parabola

egyenlete:

L2 f (x) = y0 Ha

(x − x0 ) (x − x2 ) (x − x0 ) (x − x1 ) (x − x1 ) (x − x2 ) +y1 +y2 , x ∈ [x0 , x2 ] . (x0 − x1 ) (x0 − x2 ) (x1 − x0 ) (x1 − x2 ) (x2 − x0 ) (x2 − x1 )

n túl nagy a magas fokú interpolációs polinom miatt a kilengések

is nagyok lesznek, ami a valóságnak a torzításához vezethet.

Ezért

célszer¶ relatív alacsony fokú polinommal dolgozni, például úgy, hogy csak egy pár jellegzetes alappontot választunk ki és ezek segítségével építjük fel a Lagrange -féle polinomot.

101. Példa. Határozzuk meg a Lagrange -féle interpoláló polinomot

xi 0 1 4 9 . yi 0 1 2 3 L3 f (x) = y0 l0 (x) + y1 l1 (x) + y2 l2 (x) + y3 l3 (x) a következ® adatok esetében:

ahol

(x − 1) (x − 4) (x − 9) (x − 0) (x − 4) (x − 9) , l1 (x) = (−1) (−4) (−9) (1) (−3) (−8) (x − 0) (x − 1) (x − 9) (x − 0) (x − 1) (x − 4) l2 (x) = , l3 (x) = (4) (3) (−5) (9) (8) (5) l0 (x) =

tehát

L3 f (x) =

x (x − 4) (x − 9) x (x − 1) (x − 9) x (x − 1) (x − 4) +2 +3 , x ∈ [0, 9] . 24 −60 360

6.2.1.2. Véges dierenciák

Feltételezzük hogy adottak a következ®

egyenköz¶ (ekvidisztáns) alappontok:

x0 < x1 < ... < xn , illetve egy

f

függvény értékei a fenti pontokban:

f (xi ) = yi , i = 0, ..., n. Mivel az alappontok ekvidisztánsak két egymásutáni különbsége konstans

h = xi+1 − xi .

A

h-t

lépéstávolságnak vagy növekménynek nevezzük.

6.2. INTERPOLÁCIÓ

x0

és

h

84

függvényében kifejezhet® az összes többi alappont:

x1 = x0 + h, x2 = x0 + 2h, ..., xi = x0 + ih, ...

102. Definíció. A

∆f (xi ) = f (xi + h) − f (xi ) = f (xi+1 ) − f (xi )

vagy

∆yi = yi+1 − yi , i = 0, ..., n − 1 kifejezést az

A fenti

f

∆

függvény els® véges dierenciájának nevezzük.

operátor lineáris, vagyis: -adítiv: -homogén :

∆ (f + g) = ∆f + ∆g ∆ (cf ) = c∆f

.

A magasabb rend¶ dierenciákat a következ®képpen értelmezzük:

∆2 f (xi ) = ∆ (∆f (xi )) = ∆ (f (xi+1 ) − f (xi )) = = ∆f (xi+1 ) − ∆f (xi ) = f (xi+2 ) − 2f (xi+1 ) + f (xi ) ... ∆n f (xi ) = ∆ ∆n−1 f (xi ) . A gyakorlatban a véges dierenciák kiszámításához táblázatokat használunk: xi yi ∆yi ∆2 yi ∆3 yi

x0 y0 ∆y0 = y1 − y0 ∆2 y0 = ∆y1 − ∆y0

x1 y1

∆3 y0 = ∆2 y1 − ∆2 y0

∆y1 = y2 − y1 ∆2 y1 = ∆y2 − ∆y1

x2 y2 ∆y2 = y3 − y2 x3 y3 . . .

. . .

6.2. INTERPOLÁCIÓ 103. Példa.

xi y i 0 0

∆yi

85

xi = i, i = 0, ..., 4, f (x) = x2 . ∆2 yi ∆3 yi

∆y0 = 1 1

∆y02 = 2

1

∆y03 = 0

∆y1 = 3 2

∆y12 = 2

4

∆y13 = 0

∆y2 = 5 3

∆y22 = 2

9 ∆y3 = 7

4

16

A fenti példából kit¶nnik hogy az

n-ed fokú polinomnak a n + 1-ed

rend¶ véges dierenciája egyenl® zéróval:

f ∈ Pn =⇒ ∆n+1 f (xi ) = 0. Ugyanakkor

∆n f (xi ) =konstans.

A fent értelmezett dierencia:

∆f (xi ) = f (xi + h) − f (xi ) = f (xi+1 ) − f (xi ) haladó-dierenciának nevezzük. Ugyanakkor értelmezhet® retrográd-:

∇f (xi ) = f (xi ) − f (xi − h) = f (xi ) − f (xi−1 ) , illetve centrális-dierencia:

h h cf (xi ) = f xi + − f xi − . 2 2 Ezek esetében hasonló képletek vezethet®k le.

6.2.1.3. Newton-féle interpoláló polinom

Tekintsük az

x0 , x1 , ..., xn

xi = x0 + ih, i = 0, ..., n, illetve egy függvény értékeit az adott pontokban: f (xi ) = yi , i = 0, ..., n. A feladat egy legfeljebb n-ed fokú polinom megszerkesztése ami interpolálja az (xi , yi ) , i = 0, ..., n pontokat. A polinomot a következ® alakban keressük: egyenköz¶ alappontokat:

(6.2.10)

Nn f (x) = a0 +a1 (x − x0 )+a2 (x − x0 ) (x − x1 )+...+an (x − x0 ) .. (x − xn−1 ) .

6.2. INTERPOLÁCIÓ Az

ai

86

együtthatokat az interpolációs feltételekb®l határozzuk meg:

x = x0 =⇒ a0 = y0 , x = x1 =⇒ a0 + a1 (x1 − x0 ) = y1 , x = x2 =⇒ a0 + a1 (x2 − x0 ) + a2 (x2 − x0 ) (x2 − x1 ) = y2 , ... x = xn =⇒ a0 + a1 (xn − x0 ) + ... + an (xn − x0 ) .. (xn − xn−1 ) = yn Mivel az alappontok egyenköz¶ek

=⇒

 a0        a0 +a1 h a0 +a1 2h +a2 2h2  .   .  .    a0 +a1 nh +a2 n(n − 1)h2 . . . +an n (n − 1) ...2 · 1hn

= y0 = y1 = y2 . = yn

A fenti háromszög alakú egyenletrendszer megoldása:

a0 = y 0 y1 − y0 ∆y0 a1 = = h h 2 ∆ y0 a2 = 2! h2 . . .

an

∆n y0 = =⇒ n!hn

1 ∆y0 1 ∆2 y0 (x − x0 ) + (x − x0 ) (x − x1 ) + ... 1! h 2! h2 1 ∆n y0 (6.2.12) + (x − x0 ) . . . (x − xn−1 ) n! hn Az |Nn f (x) − f (x)| abszolút hibát a (6.2.12) képletb®l mint a rákövetkez®

(6.2.11) (Nn f ) (x)

= y0 +

tag kapjuk meg:

|Nn f (x) − f (x)| =

1 ∆n+1 y0 ∆n+1 y0 1 (x − x ) .. (x − x ) = $ (x) . 0 n (n + 1)! hn+1 (n + 1)! hn+1

6.2. INTERPOLÁCIÓ Ha egy adott jelöljük:

x

Nn f polinomot kiszámítani, akkor 1 −x0 ) = x−x0 −(x = α − 1, stb. és a h

értékre kell a

x−x0 1 , tehát x−x h h

α =

87

(6.2.12)-os képlet a következ®képpen alakul: (6.2.13)

1 1 1 (Nn f ) (x) = y0 + α ∆y0 + α (α − 1) ∆2 y0 +...+ α (α − 1) .. (α − n + 1) ∆n y0 . 1! 2! n! A (6.2.13) képlet el®nye, hogy az interpolációs csomópontok nem szerepelnek explicit módon.

104. Példa. Felhasználva a harmadfokú Newton interpoláló polinomot számítsuk ki

sin (6o )

ha ismertek

xi

o 5

7

yi

0.087156

0.121869

Mivel

o

9

sin (x)

o

0.156434

11

függvény alábbi értékei:

o

0.190809

13

o

15

0.224951

o

0.258819.

6o az 5o és 7o értékek közé esik az els® négy értéket használjuk

fel a harmadfokú Newton polinom meghatározására.

Ehhez el®bb a

véges dierencia táblázatot szerkesztjük meg: xi

yi

5

0.087156

∆yi

∆2 yi

∆3 yi

0.034713 7

−0.000148

0.121869

−0.000042

0.034565 9

0.156434

-0.000190 0.034375

11

0.190809

A táblázatban szerepl® aláhúzott véges dierenciák segítségével megszerkesztjük a (6.2.12) képletnek megfelel® harmadfokú interpolációs polinomot:

0.034713 (−0.000148) (x − 5) + (x − 5) (x − 7) + 2 2! · 22 (−0000042) + (x − 5) (x − 7) (x − 9) . 3! · 23

N3 f (x) = 0.087156 +

6.2. INTERPOLÁCIÓ

88

x = 6 -ban használhatjuk a (6.2.13) képletet h = 2, x = 6, x0 = 5, α = 6−5 =⇒ 2

A polinom kiértékeléséhez ahol

sin 6o = 0.087156 + 0.5 0.034713 +

0.5 (0.5 − 1) (−0.000148) + 2

0.5 (0.5 − 1) (0.5 − 2) (−0.000042) = 0.104528. 3! Habár a Lagrange, illetve Newton-féle polinom egyforma adott pontokra, mégis, numerikus szempontból el®nyösebb Newton alakú polinommal dolgozni mert újabb alappontok bevezetése eseteben nem kell az összes számítást újból elvégezni (a háromszög alakú egyenletrendszernek köszönhet®en). Nem egyenköz¶ alappontokra a Newton polinom meghatározására az ú.n.

osztott dierenciákat

használjuk.

Habár a képletek aránylag egyszer¶ek ha az alappontok egyenköz¶, approximációs szempontból mégsem a legmegfelel®bbek.

Ennek oka

az, hogy magas fokú polinomok esetében a kilengések is nagyok. Ez jól szemléltethet® az ú.n. Runge-féle példán.

1 f (x) = 1+25x 2 , x ∈ [−1, 1] . Ha felbontjuk (n + 1) különböz® egyenköz¶ alappontra, majd

105. Példa. Legyen a

[−1, 1]

intervallumot

interpoláljuk azt kapjuk eredményként, hogy a központi részen az interpoláló polinom elég jól közelíti meg az

f

függvényt, viszont a peremen minél

magasabb a polinom fokszáma annál rosszabb a megközelítés.

Alapvet®en két módszerrel lehet a polinomiális interpoláció hátrányait kiküszöbölni: az egyik módszer abban áll, hogy egyenköz¶ alappontok helyett a Csebisev-féle polinom gyökeit használjuk.

Ezt a módszert

akkor alkalmazhatjuk ha az alappontokat szabadon választhatjuk meg. A másik módszer abban áll hogy globális interpoláció helyett szakaszos interpolációt használunk.

6.2.1.4. A Csebisev-féle els®fajú polinomok 106. Definíció.

n-ed fokú (els®fajú) Csebisev-féle polinomnak nevezzük

a következ® függvényt:

Tn : [−1, 1] → R, Tn (x) = cos (n arccos x) .

6.2. INTERPOLÁCIÓ Különböz®

n-re

89

a polinomok a következ®képpen néznek ki:

T0 (x) = 1, T1 (x) = x, T2 (x) = cos (2 arccos x) = 2 cos2 (arccos x) − 1 = 2x2 − 1, ... Felhasználva a

t = arccos x változócserét meghatározhatunk egy rekurzív

képletet a polinomok kiszámítására:

cos ((n + 1) t) + cos ((n − 1) t) = 2 cos (nt) cos t ⇐⇒ Tn+1 (x) + Tn−1 (x) = 2Tn (x) x ⇐⇒ Tn+1 (x) − 2xTn (x) + Tn−1 (x) = 0, x ∈ [−1, 1] .

(6.2.14)

Tehát felhasználva hogy

T0 (x) = 1, T1 (x) = x =⇒ T2 (x) = 2x2 − 1

amit az el®bb is megkaptunk. Kiszámítjuk a

Tn

polinom gyökeit:

Tn (x) (6.2.15)

= 0 ⇐⇒ cos (n arccos x) = 0 ⇐⇒ n arccos x =

(6.2.16) arccos x

=

(6.2.17) xk

2k + 1 π 2

2k + 1 π , k = 0, ..., n − 1 ⇐⇒ n 2 2k + 1 π = cos , k = 0, ..., n − 1. n 2

Mértanilag a Csebisev gyököket a következ®képpen szerkesztjük meg: a fels® félkört levetítjük az

Ox

polinom gyökeit.

n

egyenl® részre osztjuk majd ezeket a pontokat

tengelyre.

Ezek a vetületek szolgálják a Csebisev

Észrevehet® hogy ezeknek a s¶r¶sége nagyobb a

széleken mint középen. Ennek köszönhetik a jobb interpolációs tulajdonságot. Az összes

n-ed fokú interpoláló polinom közül az (6.2.17) alappontokra

felépített interpoláló polinomnak lesz a legjobb megközelítése.

Az összes [−1, 1] intervallumon értelmezett n-ed fokú interpolációs polinom közül annak van a legkisebb hibája amelynek az alappontjai megegyeznek az Csebisev-féle n-ed fokú polinom gyökeivel. 107. Tétel.

Ha a polinomok egy tetsz®leges értelmezve akkor egy

t=

x ∈ [a, b]

x−a+x−b b−a

intervallumon vannak

6.2. INTERPOLÁCIÓ változócserével visszavezethetjük az

x=

[−1, 1]

90

intervallumra. Hasonlóan

b−a b+a + t 2 2

változócserével kivetíthetjük a Csebisev gyököket egy tetsz®leges

[a, b]

intervallumra:

b+a b−a xk = + cos 2 2

2k + 1 π n 2

, k = 0, ..., n − 1.

Egy másik módszer a polinomiális interpoláció hátrányainak a kiküszöbölésére abban áll hogy globális interpolálás helyett szakaszos interpolálást alkalmazzunk.

6.2.2. Szakaszos interpoláció

A különbség a polinomiális és a

szakaszos interpolálás között abban áll hogy ez utóbbiban nem egy (globális) polinomot szerkesztünk hanem minden részintervallumon egy meghatározott fokszámú polinomot. Ezzel aránylag alacsonyan tudjuk tartani a polinom fokszámat tekintet nélkül az interpoláló alappontok számára.

6.2.2.1. Lineáris interpoláció

A legegyszer¶bb függvény ami interpolálja

az:

(xA , yA ) , (xB , yB ) pontokat a következ® els® fokú polinom:

F (x) = yA

x − xA xB − x + yB , h h

vagy:

F (x) = yA + (x − xA )

ahol

h = xB − xA ,

yB − yA . h

Ezt megkaphatjuk mint Lagrange -féle els® fokú interpoláló polinom vagy egyszer¶ ellen®rzéssel. Általánosításként tételezzük fel, hogy adottak az

(xi , yi )ni=0 interpoláló

pontok:

f (xi ) = yi , i = 0, ..., n. [xi , xi+1 ] lesz (hi =

Ebben az esetben két egymásutáni alappont között, vagyis az intervallumon, a lineáris interpoláló függvény a következ®

6.2. INTERPOLÁCIÓ

91

xi+1 − xi ): F (x) = yi

(6.2.18)

x − xi xi+1 − x + yi+1 , hi hi

ha

x ∈ [xi , xi+1 ]

vagy:

F (x) = yi + (x − xi )

(6.2.19) Valóban az

yi+1 − yi , x ∈ [xi , xi+1 ] hi

[xi , xi+1 ] intervallumon F

.

lineáris és teljesíti az interpoláló

feltételeket:

F (xi ) = yi , F (xi+1 ) = yi+1 . Összekapcsolva a lineáris szakaszokat az alappontokban egy tört vonalat

C 0 osztályú szakaszosan lineáris függvényt vagyis, az F függvény folytonos de a deriváltja már nem.

kapunk, vagyis egy az alappontokban

108. Példa. Szerkesszük meg a lineáris interpoláló függvényt az

xi yi

alábbi adatok esetében:

Az

[0, 1]

0

1

4

9

0

1

2

3

.

intervallumon az interpoláló függvény:

F (x) = 0 ·

1−x x−0 +1· = x, x ∈ [0, 1] . 1−0 1−0

Hasonlóan

4−x x−1 2 1 +2· = + x, x ∈ [1, 4] 4−1 4−1 3 3 9−x ·x − 4 6 1 F (x) = 2 · +3 = + x, x ∈ [4, 9] . 9−4 9−4 5 5

F (x) = 1 ·

ÁBRA 'linearisspline.eps'

6.2.2.2. Harmadfokú Hermite szakaszos interpoláció

A lineáris interpoláció

segítségével egy gyors képet alkothatunk magunknak az interpolációs folyamatról, a természetben viszont kevés olyan folyamat van ami modellezhet® lenne a lineáris interpolációval.

Ehelyett "simább" függvényekre lesz

szükségünk amit magasabb fokú polinomokkal érhetünk el.

Ehhez

viszont több adatra lesz szükségünk.

xA , xB alappontokban ismerjük nem csak az hanem a dA , dB iránytényez®ket is. Keressük azt

Tételezzük fel hogy az

yA , yB

ordinátákat

6.2. INTERPOLÁCIÓ a harmadfokú

H

92

polinomot ami eleget tesz a következ® feltételeknek:

(6.2.20)

H (xA ) = yA , H (xB ) = yB ,

(6.2.21)

H 0 (xA ) = dA , H 0 (xB ) = dB .

109. Definíció. A (6.2.21) feltételeknek eleget tev® polinomot harmadfokú Hermite -féle interpoláló polinomnak nevezzük.

Hasonlóan a Lagrange interpolációs polinommal, itt sem fogjuk használni a kanonikus bázist hanem az ú.n. Hermite -féle bázis függvényeket.

A harmadfokú Hermite -féle interpoláló polinomot a következ®képpen lehet megadni: 110. Tétel.

H (x) = yA h0 (x) + yB h1 (x) + dA h2 (x) + dB h3 (x)

(6.2.22)

ahol hi , i = 0, 3 a Hermite -féle alappolinomokat jelöli:

(6.2.23)

(6.2.24)

b−x x−a h0 (x) = Φ , h1 (x) = Φ h h b−x x−a h2 (x) = −h Ψ , h3 (x) = h Ψ h h

és (6.2.25)

Φ (t) = 3t2 − 2t3 , Ψ (t) = t3 − t2 , h = b − a.

Bizonyítás. Bizonyításként a

hi , i = 0, 3

alappolinomok bázis

jellegére szorítkozunk.

h0 (a) = Φ

b−a h

= Φ (1) = 1, h1 (a) = h2 (a) = h3 (a) = 0

Hasonlóan

h1 (b) = 1, h0 (b) = h2 (b) = h3 (b) = 0. h02 (a) = −h Ψ0

b−x h

|x=a

2 ! 1 b−x b−x = −h − 3 −2 |x=a = 1, h h h

h00 (a) = h01 (a) = h03 (a) = 0,

6.2. INTERPOLÁCIÓ

h03 (b) = h Ψ0

x−a h

|x=b

1 =h h

3

x−a h

93

2

−2

x−a h

! |x=b = 1,

h00 (b) = h01 (b) = h02 (b) = 0. Az el®bbi eset általánosításaként tételezzük fel hogy az

... < xn

alappontokban ismertek a következ® adatok:

(6.2.26)

f (xi ) = yi , f 0 (xi ) = di , i = 0, ..., n.

Ebben az esetben minden

x0 < x 1 <

[xi , xi+1 ] intervallumon külön megszerkesztjük

a harmadfokú Hermite -féle interpoláló polinomot.

Az [xi , xi+1 ] intervallumon a harmadfokú Hermite féle interpoláló polinomot a következ®képpen lehet megadni: 111. Tétel.

(6.2.27)

H (x) = yi h0 (x) + yi+1 h1 (x) + di h2 (x) + di+1 h3 (x)

ahol a hi , i = 0, 3 Hermite alappolinomokat a következ®képpen értelmezzük:

xi+1 − x x − xi (6.2.28)h0 (x) = Φ , h1 (x) = Φ hi hi xi+1 − x x − xi (6.2.29)h2 (x) = −h Ψ , h3 (x) = h Ψ hi hi

ahol a lépés hi = xi+1 − xi , illetve a Φ, Ψ függvények (6.2.25) képletben vannak értelmezve. di deriváltak ismeretlenek ezeknek az értékeit meg lehet közelíteni az xi , yi adatokból kiindulva. Például, hármasával Természetesen, ha a

csoportosítva, minden alappontban a derivált egyenl® a három pontra épített parabola iránytényez®jével (Bessel -féle interpoláció). Mivel az alappontokban a deriváltak megegyeznek az összetett Hermite -féle interpoláló polinom

C1

osztályú lesz.

112. Példa. Szerkesszünk egy Hermite -féle interpoláló függvényt a következ® adatokra:

6.2. INTERPOLÁCIÓ

xi 0 1 4 9 yi 0 1 2 3

.

94

Els®sorban minden ponthoz hozzá kell rendelni

egy iránytényez®t. Legyen az

(x1 , y1 )

pontban az iránytényez® egyenl®

a mellette lév® pontok által meghatározott egyenes iránytényez®jével:

d1 =

y2 − y0 2−0 1 = = . x2 − x0 4−0 2

d2 =

y3 − y1 3−1 1 = = . x3 − x1 9−1 4

Hasonlóan:

A

d0 -t

egyenl®vé tesszük az els® szakasz iránytényez®jével:

d0 =

y1 − y0 = 1. x1 − x0

d3 =

y3 − y2 1 = . x3 − x2 5

Hasonlóan

A (6.2.27) képletb®l a

[0, 1]

intervallumon

H (x) = 0h0 (x) + 1h1 (x) + 1h2 (x) + 1/2h3 (x) ahol

h1 (x) = Φ (x) = 3x2 − 2x3 , h2 (x) = −Ψ (1 − x) = − (1 − x)3 + (1 − x)2 , h3 (x) = Ψ (x) = x3 − x2 . Tehát

H (x) = 3x2 − 2x3 − (1 − x)3 + (1 − x)2 + Hasonlóan, az

[1, 4]

x3 − x2 , x ∈ [0, 1] . 2

intervallumon:

(4 − x)2 2 (4 − x)3 2 (x − 1)2 4 (x − 1)3 (4 − x)3 (4 − x)2 (x − 1)3 (x − 1)2 H (x) = − + − − + + − , 3 27 3 27 18 6 36 12 illetve a [4, 9] intervallumon: H (x) =

6 (9 − x)2 4 (9 − x)3 9 (x − 4)2 4 (x − 4)3 (9 − x)3 (9 − x)2 (x − 4)3 (x − 4)2 − + − − + + − 25 125 25 125 100 20 125 25

6.2.2.3. Harmadfokú spline interpoláció

Az el®bbi interpolációs eljárásnál

simább függvényt kapunk ha a másodrend¶ derivált folytonosságát is megköveteljük. Az így szerkesztett interpolációs polinom neve spline.

6.2. INTERPOLÁCIÓ

95

Ha az el®bbiekben a deriváltakat úgy határozzuk meg hogy a kapott függvény kétszer folytonosan deriválható legyen, akkor a harmadfokú (cubic) spline interpolációs módszert kapjuk. n Legyenek tehát (xi , yi )i=0 az interpoláló pontok. Minden szakaszon egy

[xi , xi+1 ]

si , i = 0, ..., n−1 (si : [xi , xi+1 ] → R) harmadfokú polinomot

szerkesztünk meg úgy, hogy az alappontokban az els® fokú és a másodfokú derivált folytonos legyen, tehát

4n ismeretlenünk lesz.

Az ismeretlenekhez

hozzárendelt egyenleteket a következ® feltételekb®l határozzuk meg:

• (n + 1)

egyenlet az interpoláló feltételekb®l:

si (xi ) = yi , i = 0, . . . , n, • (n − 1) egyenlet a folytonossági feltételekb®l a bels® alappontokban: si (xi ) = si+1 (xi ) , i = 1, . . . , n − 1, • (n − 1)

egyenlet az els® derivált folytonossági feltételekb®l a

bels® alappontokban:

s0i (xi ) = s0i+1 (xi ) , i = 1, . . . , n − 1, • (n − 1)

egyenlet a második derivált folytonossági feltételb®l a

bels®alappontokban:

s00i (xi ) = s00i+1 xi , i = 1, . . . , n − 1, összesen

(4n − 2) egyenlet.

kapjunk szükség van meg

Hogy egyértelm¶en meghatározott egyenletrendszert

2

feltételre.

Ezeket rendszerint a széleken

lev® alappontokhoz kötjük (perem feltételek).

A

4n × 4n-es

lineáris

egyenletrendszer megoldása szolgáltatja a szakaszokra értelmezett polinomok együtthatóit. A továbbiakban ismertetjük a leggyakrabban használt peremfeltételeket:

(1) Ha ismertek a

d0 , dn iránytényez®k a széleken, akkor a következ®

kikötéseket használjuk:

s00 (x0 ) = d0 , s0n−1 (xn ) = dn . Az így felépített függvényt teljes ( complete) spline -nak nevezzük.

6.2. INTERPOLÁCIÓ

96

(2) Ha ismertek a másodrend¶ deriváltak értékei

D0 , Dn a széleken,

akkor a következ® kikötéseket használjuk:

s000 (x0 ) = D0 , s00n−1 (xn ) = Dn . (3) Az egyik legrégebbi perem feltétel az ú.n. természetes feltétel:

s000 (x0 ) = 0, s00n−1 (xn ) = 0. A natural spline név a rajzoló szerkezetre (spline, splain) vezethet® vissza, ugyanis a végpontokban a fémrúd görbülete nulla (lásd a gyakorlati kivitelezését). (4) Ha a deriváltak ismeretlenek egy másik peremfeltétel az ún. not-a-knot , vagyis az els® és utolsó bels® interpolációs pont

s000 folytonos s0 = s1 illetve sn−1 =

inaktív. Ehhez kikötjük hogy a harmadrend¶ derivált legyen

x1 ,

és

xn−1 -ben.

Eredményként

sn−2 . Az összes szakaszonként értelmezett harmadfokú polinomok közül a köbös (cubic) spline függvénynek van maximális simasági foka, nevezetesen

C 2.

Ennél magasabb simasági fokkal csak a globális harmadfokú polinom

rendelkezik. A harmadfokú spline függvény általánosítása a következ®:

(xi )ni=0

(2n + 1)ed fokú spline-nak azt a függvényt nevezzük, amely minden egyes [xi , xi+1 ] szakaszon 2n + 1-ed fokú polinom, továbbá az x0 , ..., xn pontokban a deriváltak 2n-ed rendig bezárólag eleget tesznek a folytonossági feltételeknek. 113. Definíció. Az

alappontokhoz hozzárendelt

Az egyik legegyszer¶bb spline a már ismertet lineáris interpoláló függvény (a tört vonal) aminek neve lineáris spline. Minden szakaszon els®fokú polinom a folytonossági rendje pedig

C 0.

Bármilyen fokú sline függvény értelmezhet®, de a páratlan fokú spline-ok esetében az az el®ny hogy a peremfeltételeket szimmetrikusan lehet kezelni. 114. Példa. Szerkesszük meg a köbös spline interpolációs függvényt az alábbi adatok esetében különböz® peremfeltételek esetében:

6.3. KÉTVÁLTOZÓS LINEÁRIS (BILINEÁRIS) INTERPOLÁCIÓ

xi 0 1 4 9 yi 0 1 2 3

97

.

   s1 (x) , x ∈ [0, 1) s (x) = s2 (x) , x ∈ [1, 4)   s3 (x) , x ∈ [4, 9] ahol

si (x) = ai x3 + bi x2 + ci x + di , i = 0, 1, 2 A 12 együtthatót a következ® lineáris egyenletrendszerb®l számítjuk ki:

                                                

0a0 + 0b0 + 0c0 + d0 = 0 a1 + b 1 + c 1 + d 1 = 1 64a2 + 16b2 + 4c2 + d2 = 2 279a2 + 81b2 + 9c2 + d2 = 3 a0 + b0 + c0 + d0 = a1 + b1 + c1 + d1 64a1 + 16b1 + 4c1 + d1 = 64a2 + 16b2 + 4c2 + d2 3a0 + 2b0 + c0 = 3a1 + 2b1 + c1 48a1 + 8b1 + c1 = 48a2 + 8b2 + c2 6a0 + 2b0 = 6a1 + 2b1 24a1 + 2b1 = 24a2 + b2 2b0 = 0 54a2 + 2b2 = 0

ahol az utolsó két egyenletben a "természetes" peremfeltételeket használtuk. A spline függvények esetében is a kanonikus bázis felíráson kívül van más alkalmasabb bázis az ú.n. B-spline függvények.

6.3. Kétváltozós lineáris (bilineáris) interpoláció Hasonlóan az egyváltozós esettel, a kétváltozós függvények esetében felmerül a probléma hogy bizonyos diszkrét pontokból kiindulva egy interpoláló felületet kell szerkeszteni. A legegyszer¶bb módszer, a bilineáris interpoláció, akkor alkalmazható ha az adott diszkrét pontok rácsszer¶en helyezkednek el.

6.3. KÉTVÁLTOZÓS LINEÁRIS (BILINEÁRIS) INTERPOLÁCIÓ Legyen

98

n f (x, y) egy kétváltozós függvény amelynek csak (xi , yj )m, i=0,j=0

pontokban ismerjük az értékeit:

zij = f (xi , yj ) , i = 0, ..., m, j = 0, ..., n. Az

(xi , yj ) párokról azt mondjuk hogy egy m×n-es háló (rács) csomópontjait

alkotják.

x \ y y0 x0 z00

... yj z0j

. . . yn z0n

zi0

zij

zin

zm0

zmj

zmn .

. . .

xi . . .

xm

F (x, y) függvény megszerkesztése amely pontjaiban f -el és megközelíti egy tetsz®leges (x, y)

A feladat egy kétváltozós megegyezik a háló

pontban. A legegyszer¶bb módszer természetesen a (bi)lineáris interpoláció amit úgy valósítunk meg hogy mind a két tengely irányában alkalmazzuk a lineáris interpolációt i.e., egyszer az egyik, majd a másik változót tekintjük konstansnak.

(x, y) szomszédságában lév® pontokat, vagyis közrefogni a x pontot két egymásutáni Ox tengelyen lev® alapponttal: Az els® lépés meghatározni az

xi < x < xi+1 , majd hasonlóan:

yj < y < yj+1 . Az így meghatározott pontok illetve az

f

függvény értékei ezekben a

pontokban az alábbi táblázatban szerepelnek:

xy yj y yj+1 xi f (xi , yj ) = zij f (xi , yj+1 ) = zi,j+1 x f (x, y) =? xi+1 f (xi+1 , yj ) = zi+1,j f (xi+1 , yj+1 ) = zi+1,j+1. y = yj változót, ekkor az F függvény egyváltozós válik (x változó szerinti): F (x, yj ) . Az így kapott egyváltozós függvényre Rögzítjük például az

6.3. KÉTVÁLTOZÓS LINEÁRIS (BILINEÁRIS) INTERPOLÁCIÓ

99

alkalmazzuk a (6.2.19) lineáris interpolációs képletet:

F (x, yj ) = f (xi , yj ) + (x − xi )

f (xi+1 , yj ) − f (xi , yj ) , x ∈ [xi , xi+1 ] , xi+1 − xi

ahonnan (6.3.1)

F (x, yj ) = f (xi , yj ) + (x − xi )

Hasonlóan rögzített

y = yj+1

f (xi+1 , yj ) − f (xi , yj ) . xi+1 − xi

azt kapjuk hogy:

(6.3.2)

f (xi+1 , yj+1 ) − f (xi , yj+1 ) . xi+1 − xi

F (x, yj+1 ) = f (xi , yj+1 ) + (x − xi )

F (x, yj ) , F (x, yj+1 ) adatokra újból alkalmazzuk a (6.2.19) rögzített x = x-re azt kapjuk hogy:

A már ismert képletet és (6.3.3)

F (x, y) = F (x, yj ) + (y − yj )

F (x, yj+1 ) − F (x, yj ) , y ∈ [yj , yj+1 ] yj+1 − yj

vagyis: (6.3.4)

F (x, y) = F (x, yj ) + (y − yj )

F (x, yj+1 ) − F (x, yj ) . yj+1 − yj

A (6.3.4) értéknél jobb megközelítést kapunk ha a bilineáris interpoláció helyett a bicubic sémát alkalmazzuk, vagyis mindkét tengely mentén harmadfokú polinommal interpolálunk.

115. Példa. Az

f (θ, ϕ) =

´ϕ 0

√

dϕ 1−sin2 θ sin2 ϕ

elliptikus integrálra

ismert a következ® táblázatban szerepl® értékek:

θϕ 50◦ 60◦ 70◦ 80◦

50◦ 0.9401 0.9647 0.9876 1.0044

55◦ 1.05 1.0848 1.1186 1.1444

60◦ 1.1643 1.2125 1.2619 1.3014

Számítsuk ki bilineáris interpolációval

65◦ 1.2833 1.3489 1.4199 1.4810

f (55, 53)

70◦ 1.4068 1.4944 1.5959 1.6918

közelít® értékét.

6.4. FÜGGVÉNY APPROXIMÁCIÓ A LEGKISEBB NÉGYZETEK MÓDSZERÉVEL 100

Bizonyítás. Az (55,53) pont szomszédságában lév® adatok:

xy 50 50 0.9401 55 0.9524 60 0.9647

Alkalmazzuk a (6.3.2) ill. (6.3.3) képletet:

F (60, 50) − F (50, 50) = 60 − 50 0, 9401 + 5 ∗ (0, 9647 − 0, 9401)/10 = 0.9524

F (55, 50) = F (50, 50) + (55 − 50)

F (55, 55) = f (50, 55) + (55 − 60) majd rögzítjük a

x = 55

F (60, 55) − F (50, 55) = 1.0675 60 − 50

változót és alkalmazzuk a (6.3.4) képletet:

F (55, 55) − F (55, 50) = 55 − 50 0.9524 + 3 ∗ (1.0675 − 0.9524)/5 = 1.0214.

F (55, 53) = F (55, 50) + (53 − 50)

ABRA bilinterp1.eps

6.4. Függvény approximáció a legkisebb négyzetek módszerével Amint az interpolációnál láttuk, ha ismert egy

f

függvény értékei az

xi , i = 0, ..., n csomópontokban akkor felépíthet® egy F approximáció úgy, hogy F (xi ) = f (xi ) = yi . Ha viszont az f (xi ) értékek nem pontosak (például hozzávet®leges mérések alapján számított értékek) akkor a pontok interpolációja nem Az adott pontok elhelyezkedését®l függ®en a közelítés különböz® fokszámú polinommal történhet. A lineáris függvénnyel való közelítés:

F (x) = ax + b.

Az interpolációs felt következik, hogy

  ax0 + b = y0     ax1 + b = y1 . .  .     ax + b = y n n

53

55 1.05 1.0675 1.0848

6.4. FÜGGVÉNY APPROXIMÁCIÓ A LEGKISEBB NÉGYZETEK MÓDSZERÉVEL 101 vagy mátrix alakban

  x0 1 !      x1 1  a   .  =   .  b   .  xn 1    x0 1 y0     x1 1   y  , Y =  .1 X= .  .   .  .   . xn 1 yn 

ahol

 y0  y1   .  ⇐⇒ X . .  yn    .  

a b

! = Y,

A túlhatározott le megoldására

használhatjuk a 4.5 alfejezetben tárgyalt módszert .....

X tX

a b

! = X tY

vagyis

(6.4.1)

( P Pn Pn n 2 xi yi i=1 xi · b = i=1 xi · a + . Pi=1 Pn Pn n = i=1 yi i=1 1 · b i=1 xi · a +

116. Példa. A

xi 0 1 4 9 pontok esetében az (6.4.1) egyenletrendszer yi 0 1 2 3 ( 98a + 14b = 36 , 14a + 4b = 6

15 , b = 37 , a = 49 15 egyenlete: y = x + 73 . 49

és a megoldása: egyenes

tehát a pontokra legjobban illeszked®

ABRA LNM

Lineáris közelítés helyett használhatunk egy tetsz®leges a pontokat közelíthetjük egy tetsz®leges egyenletrendszer általánosítható egy

k ≤ n

(x, y)i

fokú polinommal A ()

k−ad fokú közelít® polinomra F (x) =

k

a0 x + . . . + ak−1 x + ak Az xi pontok közelítésére használhatunk egy tetsz®leges

k

ad fokú

polinomot, ebben az esetben a () egyenletrendszer a következ®képpen

6.4. FÜGGVÉNY APPROXIMÁCIÓ A LEGKISEBB NÉGYZETEK MÓDSZERÉVEL 102 alakul:

  M2k · a0 + M2k−1 a1 + . . . + Mk · ak = Vk     M2k−1 · a0 + M2k−2 a1 + . . . + Mk−1 · ak = Vk−1 . .  .     Mk · a0 + Mk−1 a1 + . . . + M0 · ak = V0 ahol

Mj =

n X

xji , j = 0, . . . , 2k,

i=1

Vj =

n X

xji · yi , j = 0, . . . . , k.

i=1 117. Példa. Az el®bbi adatokkal

 Pn Pn 2 Pn 2 Pn 3 4   Pi=1 xi · a0 + Pi=1 xi · a1 + Pi=1 xi · a2 = Pi=1 xi yi n n n n 2 3 xi y i i=1 xi · a2 = i=1 xi · a1 + i=1 xi · a0 +  Pi=1 Pn Pn  Pn 2 n = i=1 yi i=1 1 · a2 i=1 xi · a1 + i=1 xi · a0 + ⇐⇒    6818a0 + 794a1 + 98a2 = 276 794a0 + 98a1 + 14a2 = 36   98a0 + 14a1 + 4a2 = 6 a megoldás

a0 =

−37 , 1086

ABRA LNM2.eps

a1 =

673 , 1086

a2 =

30 . 181

6.4. FÜGGVÉNY APPROXIMÁCIÓ A LEGKISEBB NÉGYZETEK MÓDSZERÉVEL 103

6.2.1. ábra.

7. FEJEZET

Bézier görbék, Bézier felületek A Bézier görbék a CAGD tantárgy alapjait képezik. Habár elméleti szinten már a 30-as években jelentek meg dolgozatok e téren, az igazi fejl®dése a 60-as években kezd®dött P. de Casteljau (Citroen), P. Bézier (Renault), S. Coons (Ford), W. Gordon (General Motors), J. Ferguson (Boeing) munkásságával. Az interpolációtól eltér®en, a Bézier görbék szerkesztése a "free form" kategóriába sorolható vagyis, a formatervez® szabadon választ meg bizonyos pontokat amelyeket összekötve egy megközelít® képet adnak a görbe alakjáról.

7.1. Bézier görbe szerkesztése "divide et impera" algoritmussal Tekintsük a

P0 , P1 , P2

pontokat. Egy

P0

és

P2

pontok közé húzódó

görbét szeretnénk szerkeszteni, ugyanakkor a görbe alakját a

P1 ponttal

fogjuk befolyásolni. A

Pi

pontokat kontroll pontoknak, a

P0 P1 P2

kontroll poligonnak

nevezzük. Az algoritmus abban áll hogy a görbe pontjait sorozatos felezéssel hozzuk létre. Legyen

(1)

P0

a

P0 P 1

szakasz,

(1)

P1

pedig a

104

P 1 P2

szakasz felez®pontja.

7.1. BÉZIER GÖRBE SZERKESZTÉSE "DIVIDE ET IMPERA" ALGORITMUSSAL 105

(2) szakaszt felezzük a P0 ponttal. (2) Az így megszerkesztett P0 pont, rajta lesz a görbén.

A

(1)

(1)

P0 P1

Hogy a görbe többi pontjait megkapjuk a kapott pontokat: átnevezzük (1) (2) majd folytatjuk az osztási algoritmust. Jelöljük a P0 , P0 , P0 pontokat

P0 , P1 , P2 -vel

majd alkalmazzuk az el®bb ismertetett eljárást. (2) Az újonnan kapott P0 pont úgyszintén rajta lesz a görbén. (2) (1) Hasonlóan járunk el az els® felezésb®l származó P0 , P1 , P2 pontokkal.

ÁBRA'bezier2felezesB.eps'; (2) A három P0 pont illetve az algoritmus további alkalmazása nyomán (2) létrejöv® P0 pontok alkotják a négyzetes Bézier görbét.

7.2. NÉGYZETES BÉZIER GÖRBÉK

106

7.2. Négyzetes Bézier görbék Az el®bb említett algoritmusnál egy sajátos esetet használtunk nevezetesen a felezést; ezt most általánosítani fogjuk olyan értelemben hogy bevezetünk egy

t ∈ [0, 1]

meg.

paramétert ami a szakaszok felosztásának arányait adja

Minden egyes

például a

t = 3/4

t

értékre kapunk egy pontot a görbén.

Vegyük

értéket. Minden egyes lépésben a további pontokat

a meglév® pontok konvex kombinációjaként szerkesztjük meg:

(1)

(7.2.1)

P0

= (1 − t) P0 + tP1

(7.2.2)

(1) P1

= (1 − t) P1 + tP2 .

7.2. NÉGYZETES BÉZIER GÖRBÉK

A görbén elhelyezked®

(2)

P0

107

pontot a már meglév®

(1)

(1)

P0 , P1

pontok

konvex kombinációjaként számítjuk ki:

(2)

P0

(7.2.3) Módosítva

(1)

(1)

= (1 − t) P0 + tP1 .

t paraméter értékeit 0-tól 1-ig az el®bbi konvex kombinációk

megkapjuk a görbe összes pontját.

t paramétert®l, a további jelöléseknél (2) hogy a P0 (t) pontok írják le a görbét a

Mivel a pontok helyzete függ a ezt gyelembe vesszük. Azt következ®képpen jelöljük:

(2)

P0 (t) = P (t) , t ∈ [0, 1] . Felhasználva a (7.2.3),(7.2.2) képleteket azt kapjuk, hogy:

(2)

(1)

(1)

P (t) = P0 (t) = (1 − t) P0 (t) + tP1 (t) = (1 − t) [(1 − t) P0 + tP1 ] + t [(1 − t) P1 + tP2 ] = (1 − t)2 P0 + 2t (1 − t) P1 + t2 P2

7.3. HARMADFOKÚ BÉZIER GÖRBÉK

108

tehát a görbe egyenlete:

P (t) = (1 − t)2 P0 + 2t (1 − t) P1 + t2 P2 , t ∈ [0, 1] .

(7.2.4)

A görbe egyenletéb®l kiolvasható, hogy ez egy másodfokú görbe. A görbe legfontosabb tulajdonságai a következ®k:

• P0 , P2

a görbe két végpontja:

P (0) = P0 , P (1) = P2 ; • •

A polinom alakjából következik a görbe folytonossága; A görbe tangense a

P0 , P2 végpontokban megegyezik a P0 P1 , P1 P2

szakaszokkal:

d P (t) = −2 (1 − t) P0 + 2 (1 − 2t) P1 + 2tP2 dt ahonnan:

d d P (0) = 2 (P1 − P0 ) , P (1) = 2 (P2 − P1 ) ; dt dt •

A kontroll pontok által alkotott konvex burkoló (a mi esetünkben a

P0 P1 P2

háromszög) tartalmazza a görbét.

Az (7.2.4) egyenletben az együtthatók összege:

(1 − t)2 + 2t (1 − t) + t2 = 1 vagyis a görbe konvex kombinációja a kontroll pontoknak. Összeillesztési szempontokat tartva szem el®tt hatékonyabb a harmadfokú Bézier görbék használata.

7.3. Harmadfokú Bézier görbék A (7.2.4) képletben szerepl® együtthatókat:

(1 − t)2 = b20 (t) , 2t (1 − t) = b21 (t) , t2 = b22 (t) Bernstein-féle másodfokú alap-polinomoknak nevezzük. A fenti jelöléseket felhasználva a másodfokú Bézier görbe egyenlete felírható:

2

2

P (t) = (1 − t) P0 + 2t (1 − t) P1 + t P2 =

2 X i=0

b2i (t) Pi , t ∈ [0, 1] .

7.3. HARMADFOKÚ BÉZIER GÖRBÉK

109

Hasonlóan értelmezhet®k a harmadfokú Bernstein-féle alap polinomok:

b30 (t) = (1 − t)3 , b31 = 3t (1 − t)2 , b32 = 3t2 (1 − t) , b33 (t) = t3 , és általánosan az

n-ed

fokú Bernstein-féle alap-polinomok:

bni (t) = Cni ti (1 − t)n−i , i = 0, ..., n, ahol

t ∈ [0, 1] .

Ezeknek az alap-polinomoknak a f®tulajdonsága, hogy

az egységnek egy partícióját alkotják, vagyis:

bni

(t) ≥ 0,

n X

bni (t) = 1.

i=0 A harmadfokú Bézier -féle görbék esetében ugyanúgy járhatunk el mint a másodfokúak esetében, nevezetesen értelmezhetjük ®ket rekurrens képlettel vagy analitikusan. Legyen

P0 , P1 , P2 , P3

kontroll pontok.

118. Definíció. Mértanilag a harmadfokú Bézier görbe a (3) P3 (t) képlettel szerkeszthet® meg ahol:

( (k) Pi

és

(t) =

(k−1)

(1 − t) Pi

(k−1)

(t) + tPi+1

P (t) =

(t) , ha k = 1, 2, 3, i = 0, .., 3 − k Pi , ha k = 0

t ∈ [0, 1] . 119. Definíció. Analitikusan a harmadfokú

P (t)

Bézier görbe a

következ®képpen értelmezhet®:

(7.3.1) P (t)

=

3 X

b3i (t) Pi = b30 (t) P0 + b31 P1 + b32 P2 + b33 (t) P3

i=0 (7.3.2) ahol

= (1 − t)3 P0 + 3t (1 − t)2 P1 + 3t2 (1 − t) P2 + t3 P3

t ∈ [0, 1] .

A másodfokú görbék esetében említett tulajdonságok itt is érvényesek. A Bézier görbék el®nye hogy a

görbe alakja követi a kontroll poligon

alakját, ebb®l kifolyólag egy görbe szerkesztéséhez elégséges a kontroll pontok módosítása.

7.3. HARMADFOKÚ BÉZIER GÖRBÉK A gyakorlatban, a görbén lév® -

110

t paraméternek megfelel®- pontnak

a megszerkesztéséhez, a következ® táblázatot használjuk (de Casteljau algoritmus):

P0 (1) P1 P0 (t) (1) (2) P2 P1 (t) P0 (t) (1) (2) (3) P3 P2 (t) P1 (t) P0 (t) A táblázat jobb-alsó sarkában lév® érték adja meg a pont helyzetét a görbén

t

paraméternek megfelel®en.

t = 1/3 P0 (1, −1) , P1 (2, 2) , P2 (−2, 1) , P3 (3, 0) .

120. Példa. Szerkesszük meg az alábbi pontok esetében nek megfelel® pontot a harmadfokú Bézier görbén:

Oszlopként írva a koordinátákat a de Casteljau algoritmus a következ®képpen alakul:

P0 P1 P2 P3

1 −1 2 2 −2 1 3 0

! ! (1)

P0

= 32 P0 + 31 P1 =

! (1)

P1

= 23 P1 + 31 P2 =

! (1)

P2

= 23 P2 + 13 P3 =

4 3

!

0 2 3 5 3 −1 3 2 3

! (2)

P0

=

! (2)

P1

=

10 9 5 9 3 9 12 9

! ! (3)

P0

Az alábbi ábrán a harmadfokú Bézier görbe látható, illetve a (3) nak megfelel® P0 (t) pont.

23 27 22 27

=

t=

1 3

! .

7.4. BÉZIER FELÜLETEK

111

Komplexebb görbék esetében több pontra van szükségünk. Ebben az esetében magasabb fokú Bézier görbét szerkeszthetünk vagy- és ez a gyakoribb- összetett görbéket használunk.

Ez abban áll hogy a

pontokat csoportosítjuk, például négyesével ha köbös görbét szeretnénk szerkeszteni, majd minden csoportra megszerkesztünk egy görbét (természetesen ebben az esetben a

[0, 1] intervallum átalakítható tetsz®leges [a, b] intervallumra).

Ezeket összetett Bézier görbéknek nevezzük.

Ezen görbék esetében

P3n összeilleszkedési pontokban a simasági fokot. 1 Ha a P2 , P3 , P4 pontok kollineárisak G mértani (geometriai) folytonosságról 1 beszélünk. Ahhoz hogy a görbe C osztályú legyen még szükséges hogy a P2 P3 , illetve P3 P4 szakaszok hossza arányos legyen a két görbe gyelembe kell venni a

értelmezési tartományával.

7.4. Bézier felületek A Bézier felületek szerkesztésére -úgy mint a görbék esetében- két mód van, egy mértani szerkesztés, illetve analitikus képleten alapuló. A felület generálását a

tenzoriális szorzat esetben fogjuk vizsgálni ami azt

jelenti hogy az elemzés lebontható a két (f®)irányra. Ennek érdekében a felületeket a következ® formális denícióval értelmezzük: felületnek nevezzük a

térben mozgó és alakját változtató görbe pontjait.

7.4. BÉZIER FELÜLETEK A mozgó görbét konstans,

m

P (u) =

(7.4.1)

m-ed

m X

112

fokú Bézier görbének fogjuk tekinteni:

P i bm i (u) , u ∈ [0, 1] .

i=0

u-ra) a kontroll pontok határozzák minden eredeti Pi kontroll pont úgyszintén ezúttal egy n-ed fokún:

Ezt a görbét minden pozícióban (minden meg. Feltételezzük, hogy egy Bézier görbén mozog,

Pi = Pi (v) =

(7.4.2)

n X

Pij bnj (v) , v ∈ [0, 1] .

j=0 Összekombinálva a két képletet megkapjuk a

(u, v)

P m,n

felület képletét az

paraméternek megfelel®en:

P

(7.4.3)

m,n

(u, v) =

m X n X

Pij bnj (v) bm i (u) , (u, v) ∈ [0, 1] × [0, 1] .

i=0 j=0 A (7.4.3) képlet mátrix alakja a következ®: (7.4.4)

 P m,n (u, v) =

m bm 0 (u) . . . bm (u)

 

P00 . . .

. . . P0n . . .

Pm0 . . . Pmn Konkrétan, ha egy

(2, 3)

  

bn0 (v) . . .

  .

bnn (v)

fokú (másodfokú az egyik irányban, illetve

harmadfokú a másik irányban) Bézier felületet szeretnénk szerkeszteni akkor a

P00 , P01 , P02 , P03 , P10 , P11 , P12 , P13 , P20 , P21 , P22 , P23 kontroll pontok

segítségével meghatározzuk a felület analitikus képletét

P

2,3

(u, v) =

2 X 3 X

Pij b2i (v) b3j (u) , (u, v) ∈ [0, 1] × [0, 1] .

i=0 j=0 Az alábbi ábrán a

(2, 3)

fokú felület látható:


Rögzített Ezt nevezik

v -re

a felület

u-tól

v (=konstans)-nek

függ®,

m-ed

113

fokú Bézier görbe lesz.

megfelel® izoparametrikus görbének.

Amellett hogy aránylag egyszer¶ képleteket eredményez, a tenzoriális szorzat el®nye hogy a felület vizsgálata akármelyik iránnyal kezdhet®.

0 0 0 , P01 2 0 0 , P02 4 0 0 , P10 0 2 0 , P11 2 2 0 , P12 4 2 2 , P20 0 4 0 , P21 2 4 4 , P22 4 4 4 . Számítsuk ki az (u, v) = 1 1 , paraméternek megfelel® pontot a (2, 2) bikvadratikus Bézier felületen. 2 2 121. Példa. Adottak az alábbi pontok:

P00


114

A (7.4.4) képletet felhasználva

  2 (1 − v) P P P 00 01 02    2u (1 − u) u2  P10 P11 P12   2v (1 − v)  v2 P20 P21 P22 

P m,n

1 1 , 2 2

=



ahol

(1 − u)2

 2   =  2  1 (u, v) = 21 , 12 .

8. FEJEZET

Dierenciálegyenletek 8.1. Els®rend¶ dierenciálegyenletek A gyakorlatban gyakran szükséges olyan függvények meghatározása

y)

(jel.

amelyeknek csak a változósát ismerjük.

általában id®ben történik, független változóként

t-t

Mivel a változás használunk.

Els®rend¶ dierenciálegyenletnek nevezzük a

dy (t) = f (t, y (t)) , vagy egyszer¶sítve y˙ = f (t, y) dt típusú egyenletet, ahol t a független változó, f pedig egy el®re megadott függvény. A dierenciálegyenlet Ha

f

y = y (t) megoldását integrál görbének nevezzük.

sajátos alakú akkor elemi úton is megoldható, de gyakran ez

olyan bonyolult, hogy sokkal kényelmesebb közelít® módszerekkel egy partikuláris megoldás el®állítása. csak numerikusan lehet megoldani

Ugyanakkorbizonyos t2 Pl. y˙ = e .

egyenleteket

122. Példa. Az

y˙ (t) = 2ty = f (t, y) dierenciálegyenlet megoldása

2

y (t) = cet (c=konstans).

Az els®rend¶ közönséges dierenciálegyenleteknek végtelen sok megoldása van, amelyek egy konstansban különböznek egymástól. Ha egy plusz feltételt rendelünk a feladathoz, akkor a megoldás egyértelm¶vé válik. Ezt nevezzük Cauchy-féle feladatnak:

( (8.1.1)

y˙ (t) = f (t, y) . y (t0 ) = y0

Mértanilag ez azt jelenti, hogy a görbeseregb®l azt a görbét választjuk amelyik keresztül halad a

P0 (t0 , y0 )

ponton.

115

8.1. ELSREND DIFFERENCIÁLEGYENLETEK 123. Példa. Ha az el®bbi feladathoz hozzárendeljük a kezdeti feltételt, akkor a

c

konstans egyenl®

(

Cauchy feladatnak

2

y (t) = et

1-gyel,

116

y (0) = 1

vagyis a

y˙ (t) = 2ty y (0) = 1 a megoldása.

A dierenciálegyenletek numerikus megoldásának alapelve, hogy az integrál görbét pontszer¶en közelítjük meg

y (ti ) = yi , ahol

ti

a tanulmányozott intervallum

[t0 , tf ] 3 t

egy felosztása.

A módszereket feloszthatjuk egylépéses, illetve többlépéses módszerekre,

yi+1

attól függ®en hogy az

ordináta kiszámításához egy vagy több

el®zetes ordinátát használunk. Az egy lépéses módszerek közé tartozik az Euler, Runge-Kutta, Taylor, fokozatos közelítések módszere, míg a prediktor-korrektor többlépéses módszer.

8.1.1. Az Euler-féle (törtvonal) módszer

Tekintsük a Cauchy-

féle feladatot:

( (8.1.2)

illetve az

[t0 , tf ]

y˙ (t) = f (t, y) , t ∈ [t0 , tf ] y (t0 ) = y0

intervallum egy ekvidisztáns felosztását:

t0 < t1 < ... < tn = tf , ti+1 = ti + h ahol

h=

tf −t0 a lépést jelöli. n

Az Euler módszer lényege, hogy az elméleti görbét -pontról pontra haladva- lineáris szakaszokkal közelítjük meg és eredményül egy

P0 P1 ...Pn

törtvonalat kapunk; innen származik a módszer elnevezése. ABRA eulerde.eps Kiinduló pontként a

P0 (t0 , y0 ) használjuk, majd rendre megszerkesztjük

Pi (ti , yi ) , i = 1, n pontokat. Az y1 ordináta kiszámításához gyelembe vesszük a t0 id®pontban ismert y0 függvény értékét és feltételezzük hogy a [t0 , t1 ] intervallumban a függvény növekedése konstans, vagyis y˙ (t0 ) = f (t0 , y0 ): a

8.1. ELSREND DIFFERENCIÁLEGYENLETEK

P0 P1 :

(8.1.3)

117

y1 − y0 = f (t0 , y0 ) ⇒ y1 = y0 + h · f (t0 , y0 ) . t1 − t0

Hasonlóan szerkesztjük meg az összes többi

yi

ordinátát:

yi+1 = yi + h · f (ti , yi ) i = 1, ..., n − 1.

(8.1.4) Habár

h → 0 -ra a módszer konvergens, a konvergenciarend lineáris.

A (8.1.4) relációt a Taylor képletb®l kapjuk els® két tagjára korlátozva:

y (t + h) = y (t) +

(8.1.5) majd

1 y˙ (t) · h + O h2 , 1!

t = ti -re y (ti+1 ) = y (ti + h) = y (ti ) + y˙ (ti ) · h + O h2 yi+1 = yi + h · f (ti , yi ) + O h2 .

Tehát minden iterációban az elkövetett hiba négyzetes

n

lépés után (ennyi lépés szükséges

Pn

⇔

O (h2 ),

de

kiszámításához) a hibarend

lineárissá válik:

n × O h2 = O h−1 × O h2 = O (h) . 124. Példa. Oldjuk meg a következ® Cauchy-féle feladatot az Euler módszerrel (h

= 0.1) (

Bizonyítás.

y˙ = 2ty . y (0) = 1

f (t, y) = 2ty, t0 = 0, y0 = 1 ⇒ P0 (0, 1) .

ti = i ∗ 0.1, y1 = y0 + hf (t0 , y0 ) = y0 + h (2t0 y0 ) = 1 + 0.1 · 2 · 0 · 1 = 1 =⇒ P1 (0.1, 1) y2 = y1 + hf (t1 , y1 ) = 1 + 0.1 · 2 · 0.1 · 1 = 1.02 =⇒ P2 (0.2, 1.02) , ...

8.1.2. Javított módszerek 8.1.2.1. Taylor módszer

Az Euler módszert - mint lineáris tagú

Taylor sorfejtés - pontosíthatjuk további tagok gyelembevételével, például


118

négyzetes tag gyelembevételével:

y(t + h) = y(t) +

(8.1.6) Az

y¨(t)

1 1 h · y(t) ˙ + + h2 · y¨(t) + O h3 . 1! 2!

kiszámításához használjuk a (8.1.2) képletet:

(8.1.7)

y¨ (t) =

∂f ∂f ∂f ∂f ∂ y˙ (t) = (t, y)+ (t, y)·y(t) ˙ = (t, y)+f (t, y) (t, y) . ∂t ∂t ∂y ∂t ∂y

Tehát (8.1.8)

yi+1

1 2 ∂f ∂f = yi +h·f (ti , yi )+ h (ti , yi ) + f (ti , yi ) (ti , yi ) , i = 0, . . . , n. 2 ∂t ∂y

és a módszer hibarendje

n × O (h3 ) = O (h2 ) .

Az eljárás általánosítható viszont egyre komplikáltabb kápleteket kapunk az ...

y (t) =

y

deriváltakra, például a harmadrend¶ derivált:

∂ 2f ∂ 2f ∂f ∂ 2f (t, y)+2 (t, y)· y ˙ (t)+ (t, y)·(y˙ (t))2 + (t, y)·¨ y (t) . 2 2 ∂t ∂t∂y ∂y ∂y

8.1.2.2. Téglalap módszer

Az Euler módszer pontosságának javítására

használjuk újból a Taylor sorfejtést integrál maradékkal:

ˆ

y (t + h) = y (t) +

(8.1.9)

h

y˙ (t + s) ds, 0

majd az integrálra különböz® numerikus közelítést használunk. Az integrál közelítésére használjuk a (

ˆ

módszert:

h

0

??)

h y˙ (t + s) ds = hy˙ t + 2

(középponti) téglalap

.

A (8.1.2)-b®l

h h h y˙ t + = f t + ,y t + , 2 2 2 és gyelembe véve a Taylor képletet:

h y t+ 2

h h = y (t) + y˙ (t) = y (t) + f (t, y (t)) , 2 2


119

következik, hogy

h y˙ t + 2

h h = f t + , y (t) + f (t, y (t)) . 2 2

Visszahelyettesítve a (8.1.9) képletbe a kapott integrált azt kapjuk, hogy:

y (t + h) = y (t) + hf

(8.1.10) vagyis

t = ti

-re

y (ti+1 ) = y (ti ) + hf

(8.1.11) ami

y (ti ) = yi

h h ti + , y (ti ) + f (ti , y (ti )) , 2 2

közelítést használva a következ® képlethez vezet:

yi+1

(8.1.12)

h h t + , y (t) + f (t, y (t)) , 2 2

h h = yi + hf ti + , yi + f (ti , yi ) , i = 0, n − 1. 2 2

8.1.2.3. Trapéz módszer (Heun módszer)

A trapéz módszert használva

a (8.1.9) képletben szerepl® integrál kiszámításához :

ˆ

h

y˙ (t + s) ds = 0

h (y˙ (t) + y˙ (t + h)) , 2

azt kapjuk, hogy:

h h [y˙ (t) + y˙ (t + h)] = y (t) + [f (t, y (t)) + f (t + h, y (t + h))] = 2 2 h = y (t) + [f (t, y (t)) + f (t + h, y (t) + hf (t, y (t)))] , 2

y (t + h) = y (t) +

és

t = ti -re

(8.1.13)

y (ti+1 ) = y (ti ) +

h [f (ti , y (ti )) + f (ti + h, y (ti ) + hf (ti , y (ti )))] , 2

vagyis (8.1.14)

yi+1 = yi +

h [f (ti , yi ) + f (ti + h, yi + hf (ti , yi ))] , i = 0, n − 1. 2

A (8.1.14) képlet átírható (8.1.15)

yi+1 = yi +

1 [k1 + k2 ] , i = 0, n − 1 2


120

alakra, ahol:

k1 = hf (ti , yi ) , k2 = hf (ti + h, yi + k1 ) , és Heun vagy Runge-Kutta2 néven ismert. A két javított módszer: (8.1.12) és (8.1.14) hibarendje

8.1.2.4. Runge-Kutta4 módszer

O (h2 ) .

Az eddigi ismertetett módszereknél

igényesebb, ám pontosabb az úgy nevezett Runge-Kutta4 módszer: (8.1.16)

yi+1 = yi +

1 (k1 + 2k2 + 2k3 + k4 ) , 6

ahol

k1 k3

k1 h , = hf (ti , yi ) , k2 = hf ti + , yi + 2 2 k2 h , k4 = hf (ti + h, yi + k3 ) . = hf ti + , yi + 2 2 O (h5 ), tehát n lépés után a hibarend: n × O h5 = O h−1 × O h5 = O h4 .

Minden részintervallumon a hiba

125. Példa. Az alábbi táblázatban a

?? példa különböz® módszereknek

megfelel® összehasonlító adatai szerepelnek (négy tizedes pontossággal):

ti yi 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

(Euler) 1.0000 1.0000 1.0200 1.0608 1.1244 1.2144 1.3358 1.4961 1.7056 1.9785 2.3346

(RK2) 1.0000 1.0100 1.0407 1.0940 1.1732 1.2835 1.4324 1.6306 1.8934 2.2426 2.7091

(RK4) 1.0000 1.0101 1.0408 1.0942 1.1735 1.2840 1.4333 1.6323 1.8965 2.2479 2.7183

(pontos) 1.0000 1.0101 1.0408 1.0942 1.1735 1.2840 1.4333 1.6323 1.8965 2.2479 2.7183

8.2. ELSREND, DIFFERENCIÁL EGYENLETRENDSZEREK

121

8.1.3. Fokozatos közelítések módszere yi+1 (t) = F yi (t) , ˆ

ahol

F yi (t) =

f (t, y (s)) ds...

8.2. Els®rend¶, dierenciál egyenletrendszerek A dierenciál egyenletrendszerek tanulmányozása több szempontból is fontos.

Egyfel®l bizonyos gyakorlati folyamatok leírásánál jutunk

ilyen matematikai modellhez, másfel®l az els®rend¶nél magasabb dierenciálegyenletek visszavezethet®k dierenciál egyenletrendszerekre. Tekintsük az alábbi,

x = x (t)

és

y = y (t)

függvényeket tartalmazó

dierenciál egyenletrendszert:

( (8.2.1)

x˙ (t) = f (t, x, y) , t ∈ [t0 , tf ] . y˙ (t) = g (t, x, y)

Az adott egyenleteken kívül a függvények eleget tesznek az alábbi kezdeti feltételeknek:

( (8.2.2)

Az

x, y

x (t0 ) = x0 . y (t0 ) = y0

függvények megszerkesztéséhez bármelyik ismert módszerhez

folyamodhatunk: Euler, Runge-Kutta, Taylor, stb. Az alábbiakban az Euler módszert mutatjuk be. A tanulmányozott intervallumot felosztjuk egyenköz¶ osztópontokra:

t0 < t1 < ... < tn = tf , ti+1 − ti = h. A (8.2.1) és (8.2.2) feltételekb®l következik, hogy:

x˙ (t0 ) = f (t0 , x (t0 ) , y (t0 )) = f (t0 , x0 , y0 ) y˙ (t0 ) = g (t0 , x (t0 ) , y (t0 )) = g (t0 , x0 , y0 ) .

8.2. ELSREND, DIFFERENCIÁL EGYENLETRENDSZEREK Ezekben az osztópontokban megszerkesztjük az közelít® értékét

122

x illetve y függvény

xi = x (ti ) , yi = y (ti ) : x1 − x0 = x˙ (t0 ) = f (t0 , x0 , y0 ) , t1 − t0 y1 − y0 = y˙ (t0 ) = g (t0 , x0 , y0 ) , t1 − t0

ahonnan:

x1 = x0 + hf (t0 , x0 , y0 ) , y1 = y0 + hg (t0 , x0 , y0 ) . Hasonlóan kapjuk általánosan:

xi+1 = xi + hf (ti , xi , yi ) ,

(8.2.3)

yi+1 = yi + hg (ti , xi , yi ) , i = 1, n − 1. Többfüggvényes egyenletrendszerek esetében hasonlóan járunk el.

 x˙ (t) = −y (t) (= f (t, x, y))     y˙ (t) = x(t)(= g (t, x, y)) 126. Példa. , t ∈ [0, 2π] dierenciál  x (0) = 1    y (0) = 0 ( x (t) = cos(t) egyenletrendszer feladat pontos megoldása a trigonometrikus kör: , t∈ y (t) = sin(t) [0, 2π] . A (8.2.3) képleteket használva a következ® (diszkrét) közelítést kapjuk az x illetve y függvényekre (h = 0.1): ti xi yi 0 1 0 0.1 1 0.1 0.2 0.99 0.2 0.3 0.97 0.299 0.4 0.9401 0.396 . . .

. . .

. . .

8.2. ELSREND, DIFFERENCIÁL EGYENLETRENDSZEREK

123

A dierenciál egyenletrendszerek egyik legegyszer¶bb alakja a lineáris alak:



      x˙ 1 a11 a12 · · · x1 b1 (t)  x˙   a   x   b (t)   2  =  21 a22  2  +  2 , . . .

. . .

..

. . .

.

. . .

vagy mátrix alakban

X˙ = AX + B,

(8.2.4) ahol

A

konstansokat tartalmaz, a

B

oszlopvektor pedig

t

függvénye.

127. Példa. Az el®bbi példa felírható az alábbi alakba (x1

:=

x, x2 := y ):

0 −1 1 0

A=

vagyis

Ha az

x˙ 1 x˙ 2 !

! =

, B

0 0

0 −1 1 0 !

!

x1 x2

! ,

.

A mátrix sajátvektorai lineárisan függetlenek, akkor A diagonalizálható

és ebben az esetben a dierenciál egyenletrendszer által összekapcsolt függvények leválaszthatóak egymástól egyszer¶ dierenciálegyenleteket eredményezve.



Legyen

 λ1 0 · · ·   Λ =  0 λ2  . . .

Λ=S

−1

AS

ahol

..

A

az

mátrix átlós alakja, vagyis

.

S = [S1 S2 . . .]

az

A

Si (lineárisan független)  y1   Y =  y2  úgy, hogy mátrix



sajátvektoraiból alkotott mátrix. Ha

. . .

(8.2.5)

X = SY

akkor a (8.2.4) dierenciál egyenletrendszer átírható:

S Y˙ = ASY + B,

8.3. MAGASABB-REND DIFFERENCIÁLEGYENLETEK

124

vagyis

Y˙ = S −1 ASY + S −1 B = ΛY + B,

(8.2.6) ahol



 b1 (t)   B = S −1 B =  b2 (t)  . . . .

A (8.2.6) relációt komponensekre lebontva a következ® dierenciálegyenletekhez jutunk:

y˙ i (t) = λi yi (t) + bi (t) , i = 1, 2, ... aminek általános megoldása:

yi (t) = e

λi t

ˆ −λi t ci + e bi (t) dt .

Visszahelyettesítve a (8.2.5) képletbe megkapjuk az

128. Példa. Az el®bbi példát folytatva

√

2 2

1 1 −i i

X

Λ =

megoldást.

i 0 0 −i

! , S =

! és (8.2.6)-nek megfelel®en

y˙ 1 y˙ 2

! =

i 0 0 −i

!

y1 y2

! ⇔

y˙ 1 = iy1 , y˙ 2 = −iy2 aminek megoldása

y1 (t) = c1 eit , y2 (t) = c2 e−it . Visszahelyettesítve a (8.2.5) képletbe és felhasználva a kezdeti feltételeket √ 2 it −it következik, hogy c1 = c2 = és x1 = (e + e ) /2 = cos (t) , x2 = 2 it −it

(e − e

) /2 = sin (t) .

8.3. Magasabb-rend¶ dierenciálegyenletek Tekintsük a következ® harmadrend¶ dierenciálegyenletet: (8.3.1)

...

y = f (t, y, y, ˙ y¨) ,


125

amihez az alábbi kezdeti feltételeket társítjuk:

   y (t0 ) = y0 y˙ (t0 ) = y˙0 .   y¨ (t0 ) = y¨0

8.3.1. Visszavezetés dierenciál egyenletrendszerekre

Magasabb-

rend¶ dierenciálegyenletek megoldása visszavezethet® els®rend¶ dierenciál egyenletrendszerek megoldására. Ennek érdekében a dierenciálegyenlet rendjének megfelel® számú függvényt vezetünk be:

y1 := y, y2 := y˙ 1 , y3 := y˙ 2 . Akkor a (8.3.1) dierenciálegyenlet a következ® dierenciál egyenletrendszerrel ekvivalens:

   y˙ 1 = y2 y˙ 2 = y3   y˙ 3 = f (t, y1 , y2 , y3 )

és a kezdeti feltételek:

y1 (t0 ) = y (t0 ) = y0 := y10 y2 (t0 ) = y˙ 1 (t0 ) = y˙ (t0 ) = y˙0 := y20 y3 (t0 ) = y˙ 2 (t0 ) = y¨1 (t0 ) = y¨0 := y30 . 129. Példa. Vizsgáljuk az alábbi tompítatlan, linearizált, harmonikus mozgást:

 1   y¨ + 10 y = 0 , t ∈ [0, 10] . y (0) = 1   y˙ (0) = 2 Bevezetve az

y1 := y, y2 := y˙ 1 függvényeket, az adott dierenciálegyenlet

átalakítható az alábbi (lineáris) dierenciál egyenletrendszerré:

 y˙ 1 = y2     y˙ = − 1 y 2 10 1  y1 (0) = 1    y2 (0) = 2 amit a már ismertetett módszerrel oldunk meg.


126

8.3.2. Másodrend¶ perem dierenciálegyenlet megoldása véges dierenciákkal Tekintsük az általános másodrend¶ dierenciálegyenletet: y¨ (t) + p (t) y˙ (t) + q (t) y (t) = r (t) , t ∈ [a, b]

(8.3.2)

p, q, r ∈ C [a, b] . A lényeg újból megközelíteni a tényleges y (t) n értéket bizonyos (ti )i=0 ekvidisztáns csomópontokban: y (ti ) = yi .

ahol

Ennek érdekében felhasználjuk a numerikus deriválásból ismert képleteket, például:

y (ti+1 ) − y (ti−1 ) yi+1 − yi−1 = 2h 2h y (ti+1 ) − 2y (ti ) + y (ti−1 ) yi+1 − 2yi + yi−1 y¨ (ti ) = = h2 h2 ahol yi := y (ti ) jelölést használtuk. Behelyettesítve a (8.3.2) képletbe és a pi = p (ti ) , qi = q (ti ) , ri = r (ti ) jelöléseket használva a következ® y˙ (ti ) =

lineáris egyenletrendszert kapjuk:

yi+1 − yi−1 yi+1 − 2yi + yi−1 + pi + qi yi = ri , i = 1, ..., n − 1 2 h 2h ami átrendezve: (8.3.3)

pi h 1− 2

2

yi−1 − 2 − qi h

pi h yi + 1 + 2

yi+1 = ri h2 , i = 1, ..., n−1.

8.3. MAGASABB-REND DIFFERENCIÁLEGYENLETEK Figyelembe véve, hogy

y0 , yn

127

adott, a (8.3.3) lineáris egyenletrendszer

mátrix alakja:

          

− (2 − q1 h2 ) 1 + p12h 0 p2 h 2 1− 2 − (2 − q2 h ) 1 + p22h 0 1 − p32h − (2 − q3 h2 )

0 0 1 + p32h

0 0 0 ..

0

0

0 

y1   y2     . ·  ..      yn−2 yn−1

0 

0

.

0

0



            =           

0 r0 h2 − 1 − p12h y0 r1 h2 . . .

rn−2 h2 h rn−1 h2 − 1 − pn−1 yn 2

0 0 0

1−

pn−2 h 2

0        .      

2

− (2 − qn−2 h ) 1 pi h 1− 2 − (2

Irodalomjegyzék [1] Bahvalov, N.Sz., A gépi matematika numerikus módszerei, M¶szaki Könyvkiadó, Budapest, 1977. [2] Bajcsay, P., Numerikus analízis, Tankönyvkiadó, Budapest, 1978. [3] Bjezikovics, Ja.Sz., Közelít® számítások, Tankönyvkiadó, Budapest, 1952. [4] Bjoerck A., Dahlquist G., Numerical mathematics and scientic computation, vol.1, SIAM, 2008. [5] de Boor, C.: A practical guide to splines, Springer-Verlag New York Heidelberg Berlin, 1978. [6] Conte S.D., de Boor C., Elementary numerical analysis. An algorithmic approach, 3ed., McGraw-Hill, 1980. [7] Coman, Gh.: Analiz numeric , Editura Libris, Cluj-Napoca, 1995. [8] Du, I.S., Erisman, A.M., Reid, J.K.: Direct methods for sparse matrices, Oxford Science Publications, 1989. [9] Farin, G., Hansford, D.: The essentials of CAGD, A.K. Peters, 2000. [10] Forsythe G.E., Malcolm M.A., Moler C.B., Computer methods for mathematical computations, Prentice-Hall, 1977. [11] Golub, G.H., van Loan, C.F.: Matrix computations, The Johns Hopkins University Press, (3ed.) 1996. [12] Homan J.D.: Numerical methods for engineers and scientists, (2ed.), M.Dekker, 2001. [13] Kahaner, D., Moler, C., Nash, S.: Numerical methods and software, Prentice Hall, 1988. [14] Laz r, I.: Metode numerice cu funcµii în C++, Presa universitar clujean , Cluj-Napoca, 2001. [15] Mathews, J.H, Fink, K.D., Numerical methods using Matlab, (3,ed.) PrenticeHall, 1999. [16] Moler, C.: Numerical Computing with Matlab, SIAM, 2008. [17] Obádovics, J. Gy.: Gyakorlati számítási eljárások, Gondolat Kiadó, 1972. [18] Popper, Gy., Csizmás, F.: Numerikus módszerek mérnököknek, Akadémiai Kiadó, Typotex, Budapest, 1993. [19] Press W., Teukolsky S., Vetterling W., Flannery B., Numerical recipes in C, Cambridge Univ. Press, 1992. 128

IRODALOMJEGYZÉK

129

[20] Singiresu, R.: Applied numerical methods for engineers and scientists, New Jersey Prentice-Hall, 2002. [21] Stoer J., Bulirsch, R., Introduction to Numerical Analysis, Springer New York, 2002. [22] Stoyan, G., Takó, G.: Numerikus módszerek, Typotex, 2005. [23] Strang G.: Introduction to applied mathematics, Wellesley-Cambridge, 1986. [24] Van Loan, Ch.F., Introduction to Scientic Computing, Prentice-Hall, Upper Sad- dle River, NJ, 2000. [25] Virágh, J.: Numerikus matematika, Szeged JATE Press, 2003. [26] Vladislav, T., Ra³a, I.: Analiz numeric , Editura Tehnic , Bucure³ti, 1997. [27] ***, Analiz numeric -Lucr ri de laborator, Univ. Babe³-Bolyai, Fac. de Matematic ³i Informatic , Cluj Napoca, 1994.

Numerikus módszerek - jegyzet. Kupán Pál

Recommend Documents