Numerieke Analyse [version 26 Dec 2003 ed.] [PDF]

Zitiervorschau

VE NI U

RSIT EIT

BR

S US EL

VRIJ E

Vrije Universiteit Brussel Faculteit Toegepaste Wetenschappen

V IN

N

NT

EB

RAS

SCIE

IA

CERE

TE

Numerieke Analyse

S. Caenepeel

Syllabus bij de cursus “Numerieke Analyse”(later “Numerieke Algoritmen”) Tweede Kandidatuur Burgerlijk Ingenieur Eerste Jaar Brugprogramma Burgerlijk Ingenieur

1991-1999

Inhoudsopgave 1

2

3

4

Inleiding

4

1.1

Probleemstelling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4

1.2

Fouten en foutenbronnen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

7

1.3

Voorstelling van re¨ele getallen en de totale rekenfout . . . . . . . . . . . . . . . .

7

1.4

Het conditiegetal en de onvermijdbare fout . . . . . . . . . . . . . . . . . . . . . . 11

Het oplossen van vergelijkingen

13

2.1

Algoritmen voor het oplossen van vergelijkingen . . . . . . . . . . . . . . . . . . 13

2.2

Convergentie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20

2.3

De orde van een iteratieve methode . . . . . . . . . . . . . . . . . . . . . . . . . . 23

2.4

Versnelling van de convergentie: de methode van Aitken . . . . . . . . . . . . . . 25

Het oplossen van stelsels

29

3.1

Algoritmen voor het oplossen van stelsels . . . . . . . . . . . . . . . . . . . . . . 29

3.2

Metrieken op Rn

3.3

Convergentie van de iteratieve formules . . . . . . . . . . . . . . . . . . . . . . . 34

3.4

De methode van Bairstow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31

Lineaire stelsels

44

4.1

De iteratieve methode van Gauss . . . . . . . . . . . . . . . . . . . . . . . . . . . 44

4.2

Permutatiematrices en Frobeniusmatrices . . . . . . . . . . . . . . . . . . . . . . 45

4.3

Gauss eliminatie en LU-decompositie . . . . . . . . . . . . . . . . . . . . . . . . 48

4.4

Tridiagonale stelsels . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53

4.5

Singuliere waarden ontbinding . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 1

4.6 5

6

7

8

9

Veralgemeende inversen van matrices . . . . . . . . . . . . . . . . . . . . . . . . 58

Eigenwaarden en eigenvectoren

62

5.1

Algemene eigenschappen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62

5.2

De methode van de machten . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64

5.3

De methode van Jacobi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68

Interpolatie

73

6.1

De formule van Lagrange . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74

6.2

Het algoritme van Neville . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76

6.3

De methode der gedeelde verschillen . . . . . . . . . . . . . . . . . . . . . . . . . 77

6.4

Differentierekening . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78

6.5

De minimaxbenadering . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82

Benadering volgens de methode van de kleinste kwadraten

87

7.1

Algemene methode . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87

7.2

Benadering met Forsythe veeltermen . . . . . . . . . . . . . . . . . . . . . . . . . 88

Numeriek Afleiden

94

8.1

Afleiden van de interpolatieformule van Lagrange . . . . . . . . . . . . . . . . . . 94

8.2

Praktische formules . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95

Numerieke Integratie

98

9.1

De Newton-Cotes formules . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98

9.2

De rest in de formule van Newton-Cotes . . . . . . . . . . . . . . . . . . . . . . . 103

9.3

Rombergintegratie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104

9.4

De integratieformules van Gauss . . . . . . . . . . . . . . . . . . . . . . . . . . . 106

10 Differentiaalvergelijkingen

114

10.1 De methode van Euler-Heun . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114 10.2 De methode van Runge-Kutta . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115 10.3 De orde van een eenstapsmethode . . . . . . . . . . . . . . . . . . . . . . . . . . 118 10.4 Meerstapsmethodes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120 2

10.5 Differentiaalvergelijkingen met randvoorwaarden . . . . . . . . . . . . . . . . . . 121 10.6 Parti¨ele differentiaalvergelijkingen . . . . . . . . . . . . . . . . . . . . . . . . . . 124

3

Hoofdstuk 1 Inleiding 1.1

Probleemstelling

In de cursussen “Analyse¨en “Algebra en Meetkunde”hebben we technieken gezien die toelaten om allerlei problemen op te lossen: • oplossen van algebraische vergelijkingen; • berekenen van afgeleiden en bepaalde integralen; • oplossen van lineaire stelsels; • berekenen van eigenwaarden en eigenvectoren; • integreren van differentiaalvergelijkingen; • ··· Wanneer we deze methoden in de praktijk brengen kunnen allerlei moeilijkheden optreden: 1) Het kan zijn dat ons probleem niet analytisch oplosbaar is. De vergelijking tg x = x heeft een oneindig aantal oplossingen, maar we kunnen enkel de triviale oplossing x = 0 exact uitrekenen. Vele integralen kunnen we niet oplossen. Sommige ervan zijn in de praktijk wel belangrijk. Als voorbeeld geven we de integraal Φ(x) =

Z x

exp(−t 2 /2)dt

0

We kunnen deze integraal niet uitrekenen; in de cursus “Waarschijnlijkheidsrekening en Statistiek”hebben we wel het belang ervan gezien. 2) Voor sommige problemen bestaan er oplossingsmethoden die in theorie altijd een oplossing bieden. Voor een lineair stelsel van n vergelijkingen met n onbekenden, met een reguliere matrix, hebben we de regel van Cramer. Zoals we reeds benadrukten in de cursus “Algebra en Meetkunde¨ıs de regel van Cramer totaal onbruikbaar als n groot wordt, omdat het zelfs met grote computers 4

onmogelijk is om determinanten van grote omvang (bijvoorbeeld 1000 × 1000) uit te rekenen. 3) Zelfs indien we een exacte oplossingsmethode hebben die niet al te veel rekenwerk vraagt, dan kunnen nog andere problemen optreden. Om dit te illustreren geven we een eenvoudig voorbeeld. Voorbeeld 1.1.1 Gevraagd wordt om de bepaalde integraal Z 1

I10 =

x10 ex−1 dx

0

te berekenen. Theoretisch is er geen enkel probleem. Met parti¨ele integratie vinden we Z 1

In =

xn ex−1 dx

0

Z 1  n x−1 1 = x e −n xn−1 ex−1 dx 0 0

= 1 − nIn−1

(1.1)

Het is duidelijk dat I0 = 1 − 1/e, en we kunnen de In recursief uitrekenen. L. Van Hamme voerde de berekening uit met behulp van een programmeerbare zakrekenmachine, en vond I0 = 0.632 121 I2 = 0.264 242 I4 = 0.170 904 I6 = 0.127 120 I8 = 0.118 720

I1 = 0.367 879 I3 = 0.207 274 I5 = 0.145 480 I7 = 0.110 160 I9 = −0.068 480

Hier gaat duidelijk wat mis, want het is onmogelijk dat I9 < 0. Als we dezelfde berekening uitvoeren op een Macintosh Performa 5300 computer, met behulp van het programma Matlab, versie 4.2c, dan vinden we I0 = 0.63212055882856 I1 = 0.36787944117144 I2 = 0.26424111765712 I3 = 0.20727664702865 I4 = 0.17089341188538 I5 = 0.14553294057308 I6 = 0.12680235656152 I7 = 0.11238350406936 I8 = 0.10093196744509 I9 = 0.09161229299417

I10 = 0.08387707005829 I11 = 0.07735222935878 I12 = 0.07177324769464 I13 = 0.06694777996972 I14 = 0.06273108042387 I15 = 0.05903379364190 I16 = 0.05545930172957 I17 = 0.05719187059731 I18 = −0.02945367075154

We vinden duidelijk andere resultaten. Na 18 stappen vinden we echter weer een negatieve waarde. De verklaring hiervoor is de volgende: beide rekenmachines werken slechts met een eindig aantal decimalen. Bij elke stap wordt dus een afrondingsfout gemaakt. Bij elke stap wordt deze fout 5

vermenigvuldigd met n, zodat de fout na n stappen gelijk is aan n! maal de afrondinsfout op I0 . De fout wordt dus snel groot. Om betere resultaten te verkrijgen moeten we (1.1) op een andere manier toepassen. We beweren dat In =

1 1 1 1 − + − + · · · (1.2) n + 1 (n + 1)(n + 2) (n + 1)(n + 2)(n + 3) (n + 1)(n + 2)(n + 3)(n + 4)

We bewijzen (1.2) per inductie op n. Voor n = 0 hebben we ∞

I0 =

(−1)n ∑ n! = 1 − e−1 n=1

Als (1.2) klopt voor n − 1, dan hebben we  1 1 1 + − +··· n n(n + 1) n(n + 1)(n + 2) n(n + 1)(n + 2)(n + 3) 1 1 1 1 + − + −··· = 1− n + 1 (n + 1)(n + 2) (n + 1)(n + 2)(n + 3) (n + 1)(n + 2)(n + 3)(n + 4)

1 − In = nIn−1 = n

1

−

zodat (1.2) ook klopt voor n. (1.2) is een alternerende reeks, en de som kan benaderd worden door de parti¨ele sommen uit te rekenen. We vinden zelfs een boven- en een ondergrens voor de reekssom. Als we (1.2) toepassen voor n = 10, dan vinden we voor de vijfde parti¨ele som s5 = 0.08387723387723 De eerste zes decimalen zijn correct! We zien hier dus dat een benaderde methode (reekssom vervangen door een parti¨ele som) een beter resultaat geeft dan een exacte methode! Een alternatieve manier bestaat erin om (1.1) als volgt te herschrijven: In−1 =

1 − In n

(1.3)

Bij elke stap wordt de fout nu gedeeld door n. We vertrekken van I18 = 0. Dit is natuurlijk niet correct, maar de absolute fout die we maken is klein, aangezien In klein wordt voor n groot. Als we terugrekenen tot I10 , dan wordt de absolute fout op I10 klein. We vinden volgende tabel I18 = 0.000 000 000 I16 = 0.055 555 556 I14 = 0.062 731 481 I12 = 0.071 773 250 I10 = 0.083 877 070

I17 = 0.055 555 556 I15 = 0.059 027 778 I13 = 0.066 947 751 I11 = 0.077 352 229

Als we vertrekken van I15 = 0, dan vinden we I10 = 0.08387723387723 en I9 = 0.09161227661228 6

1.2

Fouten en foutenbronnen

Het oplossen van een practisch probleem (bijvoorbeeld in de natuurkunde, of in de economie) gebeurt in verschillende stappen. 1) Eerst wordt er een wiskundig model gemaakt. Dit model is uiteraard slechts een benadering van de werkelijkheid. De fouten die optreden ten gevolge hiervan noemen we modelfouten. Indien we de modelfouten willen verkleinen, dan moeten we het model vervangen door een meer verfijnd model. Een bekend voorbeeld is de (algemene) relativiteitstheorie, die een verfijning is van de klassieke Newtoniaanse mechanica. Voor de numerieke analyst zijn de modelfouten in zoverre van belang dat bij een grof model een erg nauwkeurige berekening weinig zin heeft. 2) Bij iedere berekening worden gegevens verwerkt om resultaten te verkrijgen. Deze gegevens bevatten fouten, omdat ze het resultaat zijn van metingen of van vorige berekeningen. Dit leidt tot een fout op het eindresultaat die men de onvermijdbare fout noemt. 3) Van het wiskundig model wordt een numeriek model gemaakt, bijvoorbeeld door het probleem te discretiseren of te lineariseren. De fout die hierbij optreedt noemen we de methodefout. In een verder hoofdstuk zullen we zien dat een bepaalde integraal kan benaderd worden door de volgende formule (de regel van Simpson): Z x2 x0

h5 (4) h f (x)dx = ( f (x0 ) + 4 f (x1 ) + f (x2 )) − f (ξ) 3 90

waarbij x2 − x1 = x1 − x0 = h, en ξ ∈ (x0 , x2 ). De methodefout is h5 (4) R = − f (ξ) 90 4) In berekeningen worden re¨ele getallen steeds voorgesteld door een eindig aantal cijfers. Dit veroorzaakt een fout die men de afrondingsfout noemt. Bij elke nieuwe bewerking wordt een nieuwe afrondingsfout ingevoerd. Deze afrondingsfouten speelden een cruciale rol in voorbeeld 1.1.1. We gaan hier verder op in in de volgende paragraaf. De gecumuleerde fout op het eindresultaat die ontstaat ten gevolge van de methodefouten en de afrondingsfouten noemen we de totale rekenfout.

1.3

Voorstelling van re¨ele getallen en de totale rekenfout

De moeilijkheden in voorbeeld 1.1.1 treden op tengevolge van de afrondingsfouten. Als we zouden beschikken over een rekenmachine die met een oneindig aantal cijfers werkt (en dus re¨ele getallen exact kan voorstellen), dan zouden we deze problemen niet hebben. Onze rekenmachines werken meestal electronisch; men kan zich echter ook andere soorten rekenmachines inbeelden, bijvoorbeeld mechanische. In de limiet beschouwd is een telraam een mechanische rekenmachine. Eigenlijk komen de moeilijkheden hierop neer dat een rekenmachine werkt met slechts een eindig aantal getallen. Deze getallen noemt men machinegetallen. Voor een telraam is dit aantal klein (bijvoorbeeld 100), voor de hedendaagse electronische rekenmachines is dit aantal natuurlijk veel groter. 7

Electronische rekenmachines werken steeds in het tweetallig stelsel. Ook hier kan men zich machines inbeelden die in andere talstelsels werken. In de navolgende beschouwingen zullen we ervan uitgaan dat onze machine in het tientallig stelsel werkt, omdat we daar nu eenmaal zelf mee werken. Bij de allereerste zakrekenmachines werd dikwijls met een vaste komma gewerkt. Dit betekent dat elke machinegetal bestaat uit een vast aantal cijfers (bijvoorbeeld 8), waar de komma ergens in dit getal staat. Een machinegetal is dan van de vorm m = ±i1 i2 i3 i4 .i5 i6 i7 i8 waarbij i j ∈ {0, · · · , 9}, en waar de komma ook op een andere plaats kan staan. In moderne rekenmachines gebruikt men vrijwel steeds de drijvende komma representatie (floating point representation). De machinegetallen zien er dan als volgt uit: m = ±.i1 i2 . . . it × 10e

(1.4)

waarbij i j ∈ {0, · · · , 9}, i1 6= 0 (om de representatie uniek te maken) en e een geheel getal dat in absolute waarde kleiner is dan 10q . Men noemt m = ±0.i1 i2 . . . it de mantisse, en e de exponent. We zeggen dat we een drijvende komma representatie hebben met t cijfers in de mantisse en q cijfers in de exponent. Het getal 0 valt niet onder (1.4), maar we noemen het per definitie ook een machinegetal. Neem bijvoorbeeld t = 3 en q = 2. In de praktijk is t natuurlijk veel groter, maar voor de eenvoud beperken we ons tot een kleine waarde van t. We hebben dan bijvoorbeeld 1 = .100 × 101 1/2 = .5 × 100 De som x1 + x2 van twee machinegetallen is niet noodzakelijk een machinegetal, en hetzelfde geldt voor product en quoti¨ent. Er zijn twee mogelijke problemen. 1) Merk op dat het grootste machinegetal xmax gegeven wordt door de formule q −1

xmax = .99 . . . 9 × 1010

= (1 − 10−t )1010

q −1

Het is duidelijk dat xmax + xmax geen machinegetal is. Dit verschijnsel noemen we overflow. Analoge problemen hebben we met het kleinste strikt positieve machine getal xmin : q

−q

xmin = .100 · · · 0 × 101−10 = 10−10

en duidelijk is xmin xmin geen machinegetal. Dit verschijnsel wordt underflow genoemd. 2) Het is ook mogelijk dat x1 + x2 tussen twee machinegetallen invalt. Neem bijvoorbeeld t = 3 en q = 2, en stel x1 = .333 × 100 = .333 x2 = .100 × 10−3 = .0001 x1 + x2 = .3331 × 100 is geen machinegetal. Neem een re¨eel getal x, en onderstel dat x binnen het bereik van de machine ligt, dit wil zeggen dat 8

xmin ≤ |x| ≤ xmax . We noteren fl(x) voor het machinegetal dat het dichts bij x ligt; indien x net op de helft tussen twee opeenvolgende machinegetallen ligt, dan spreken we bijvoorbeeld af om naar boven af te ronden. Welke fout maken we als we x vervangend door fl(x)? Onderstel dat x = m × 10e , met .1 ≤ |m| < 1 en |e| < 10q . Dan is 1 |x − fl(x)| ≤ 10−t 10e 2 Aangezien x ≥ 10e−1 hebben we voor de relatieve fout x − fl(x) 1 −t e 1−e 101−t ≤ 10 10 10 = 2 x 2 Het getal η = 101−t /2 noemt men de machinenauwkeurigheid. We kunnen nu de som,het verschil, het product en het quoti¨ent van twee machinegetallen x1 en x2 defini¨eren: x1 ⊕ x2 x1 x2 x1 ∗ x2 x1 ÷ x2

= = = =

fl(x1 + x2 ) fl(x1 − x2 ) fl(x1 x2 ) fl(x1 /x2 )

Aangezien (a ⊕ b) − (a + b) ≤η a+b hebben we

a ⊕ b = (a + b)(1 + ε)

(1.5)

waarbij |ε| ≤ η. Op analoge wijze hebben we a b = (a − b)(1 + ε1 ) a ∗ b = ab(1 + ε2 ) a ÷ b = (a/b)(1 + ε3 )

(1.6) (1.7) (1.8)

waarbij |εi | ≤ η. De bewerkingen met machinegetallen voldoen niet langer aan de gebruikelijke eigenschappen. Zo is de som van machinegetallen niet associatief, zoals blijkt uit volgend voorbeeld. Voorbeeld 1.3.1 We werken met 3 cijfers in de mantisse. Stel x = .333 × 100 = .333 ; y = z = .400 × 10−3 = .0004 Als we x + y + z op de gebruikelijke manier uitrekenen, dan vinden we x + y + z = .3338

9

We vinden achtereenvolgens x ⊕ y = .333 × 100 en (x ⊕ y) ⊕ z = .333 × 100 en y ⊕ z = .800 × 10−3 en x ⊕ (y ⊕ z) = .334 × 100 Deze laatste benadering is beter! Een vuistregel is de volgende. Wanneer men een groot aantal getallen moet optellen, dan is het aan te bevelen om de getallen volgens grootte te rangschikken, en ze dan van klein naar groot op te tellen. Voorbeeld 1.3.2 Gevraagd wordt om y = 1 − x2 te berekenen, voor x dicht bij 1 gelegen. De voor de hand liggende manier om y te berekenen is de volgende. y = 1 x∗x

(1.9)

Gebruik makende van (1.6) en (1.7) vinden we y = (1 − x2 (1 + ε1 ))(1 + ε2 ) = (1 − x2 ) + (1 − x2 )ε2 − x2 ε1 (1 + ε2 ) ≈ y + yε2 − ε1 We gebruikten het feit dat ε1 ε2 verwaarloosbaar is ten opzichte van ε1 (omdat η klein is), en ook dat x2 ≈ 1. Nu volgt dat y − y ε1 y ≈ ε2 − y Hoe dichter x bij 1 ligt, hoe dichter y bij 0 ligt, en hoe groter de relatieve fout op y! Een beter algoritme om y te benaderen is het volgende. Bereken y = (1 x) ∗ (1 ⊕ x)

(1.10)

We vinden nu dat y = = = ≈ en

(1 x) ∗ (1 ⊕ x) ((1 − x)(1 + ε1 )(1 + x)(1 + ε2 ))(1 + ε3 ) y(1 + ε1 )(1 + ε2 )(1 + ε3 ) y(1 + ε1 + ε2 + ε3 )

(1.11) (1.12) (1.13) (1.14)

y − y y ≈ |ε1 + ε2 + ε3 | ≤ 3η

en het blijkt dat de relatieve fout op y van dezelfde grootte orde is als die op x. Het algoritme (1.10) is dan ook beter dan (1.9); we zeggen dat (1.10) numeriek stabiel is. Laten we alles verifi¨eren aan de hand van een numeriek voorbeeld. We werken met 4 cijfers in de mantisse. Neem x = 1.002 = .1002 × 10. Dan is x ∗ x = .1004 × 10 en y = 1 x ∗ x = −.0004 × 10 = −.4000 × 10−3 10

Op dezelfde manier vinden we dat 1 x = .2000 × 10−2 , 1 ⊕ x = .2002 × 10 en y = (1 x) ∗ (1 ⊕ x) = −.4004 × 10−3 en we zien dat y een betere benadering is van y dan y. Wat er mist gaat met het algoritme (1.9) is het feit dat we twee getallen die ongeveer aan mekaar gelijk zijn (1 en x2 ) van elkaar aftrekken. Dit gaat gepaard met een groot verlies in relatieve nauwkeurigheid. In (1.10) vermijden we deze aftrekking, en dit is de reden waarom we een beter resultaat verkrijgen. Een vuistregel is dat we bij het opstellen van een numeriek algoritme steeds trachten om het van elkaar aftrekken van twee ongeveer aan elkaar gelijke getallen te vermijden. De analyse die we verrichtten in voorbeeld 1.3.2 noemen we een voorwaartse foutenanalyse: we kennen een bovengrens voor de fouten op de gegevens, en trachten hieruit een bovengrens te vinden voor de fouten op het eindresultaat. Het spreekt vanzelf dat zulk een foutenanalyse voor meer ingewikkelde algoritmen niet altijd eenvoudig is. Een analoog vraagstuk is dat van de achterwaartse foutenanalyse: men wenst de eindresultaten met een zekere nauwkeurigheid te bepalen, en men vraagt welke fouten men mag dulden op de gegevens, en met welke nauwkeurigheid de tussenresultaten moeten berekend worden. Dit probleem is in het algemeen moeilijker dan het vorige.

1.4

Het conditiegetal en de onvermijdbare fout

Beschouw het probleem waarin een numerieke functie y = f (x) dient berekend te worden. We onderstellen dat ons probleem goed gesteld is. Hiermee bedoelen we dat y eenduidig bepaald is op een omgeving van x, en ook dat y continu afhangt van x op een omgeving van x. Onderstel een kleine afwijking δx op de gegevens. Zoals we hierboven gezien hebben is de onvermijdbare fout de afwijking δy die tengevolge hiervan ontstaat op het eindresultaat y. |δy| = | f (x + δx) − f (x)| ≈ | f 0 (x)δx| voor de relatieve onvermijdbare fout vinden we 0 δy x f (x) δx ≈ y y x Het getal c = |x f 0 (x)/y| noemen we het conditiegetal. Hoe groter c, hoe groter de onvermijdbare fout. Als we onderstellen dat de relatieve fout op de gegevens enkel ontstaat door de afronding bij het invoeren van de gegevens, dan wordt deze relatieve fout begrensd door de machinenauwkeurigheid η. In dit geval vinden we voor de onvermijdbare fout δy ≤ cη y We noemen ons probleem goed geconditionneerd als c niet te groot is (bijv. c ≤ 100). Anders zeggen we dat het probleem slecht geconditionneerd is. 11

Voorbeeld 1.4.1 We keren terug naar de situatie in voorbeeld 1.3.2. Voor het conditiegetal vinden we 2 x c = 2 1 − x2 en we zien dat het probleem slechter geconditionneerd is naarmate x dichter bij 1 komt te liggen. Neem bijvoorbeeld x = 1.0023 Dan is fl(x) = 1.002, en 1 − x2 = −.460529 × 10−2 1 − fl(x)2 = −.400400 × 10−2 De relatieve onvermijdbare fout bedraagt 13 %. Voorbeeld 1.4.2 Bekijk het lineaire stelsel  x+y = 2 1.01x + y = 2.01 De oplossing van dit stelsel is x = y = 1. Als we de co¨effici¨enten van het stelsel lichtjes wijzigen, dan heeft dit aanzienlijke gevolgen voor de oplossing:  x+y = 2 1.001x + y = 2.01 heeft als oplossing x = 10, y = −8. In het algemeen heeft 

x+y = 2 (1 + ε)x + y = 2.01

als determinant −ε. Voor ε klein is de determinant ook klein, zodat we moeilijkheden kunnen verwachten. De algemene oplossing is x = .01/ε en y = 2 − .01/ε. Een kleine wijziging in ε brengt een groot verschil in de oplossing teweeg. Ons probleem is dan ook slecht geconditionneerd. Over de bibliografie [2] is het meest elementaire werk uit de referentielijst; de inhoud ervan is grotendeels parallel met die van deze syllabus. Andere - meer uitgebreide - referentiewerken voor numerieke analyse zijn [3] and [4]. [1] is het standaardwerk voor numerieke lineaire algebra.

12

Hoofdstuk 2 Het oplossen van vergelijkingen 2.1

Algoritmen voor het oplossen van vergelijkingen

Gegeven is een continue functie f : [a, b] → R, en gevraagd wordt om de nulpunten van de functie f te bepalen. Onze methode bestaat er in om een rij (xn ) te zoeken die naar een oplossing van f (x) = 0

(2.1)

convergeert. De rij (xn ) wordt recursief gedefinieerd: Men vertrekt van een startwaarde x0 , en men berekent xn+1 uit xn met behulp van een formule van de vorm xn+1 = φ(xn ) In sommige gevallen heeft men meer ingewikkelde recursieformules; soms is de recursieformule van de vorm xn+1 = φ(xn , xn−1 ) Volgende vragen stellen zich: 1) Hoe vinden we een gepaste iteratieve rij (xn )? 2) Is de gevonden rij (xn ) convergent? 3) Convergeert de rij (xn ) naar een oplossing x van (2.1)? 4) Hoe groot is de fout |xn − x|. In deze paragraaf geven we een antwoord op de eerste vraag. Een algemene formule Neem een willekeurige continue functie k : [a, b] → R die nergens nul wordt. (2.1) is dan equivalent met x = x + k(x) f (x) (2.2) De oplossingen van van (2.2) zijn de dekpunten van de operator O : [a, b] → R : x 7→ x + k(x) f (x) 13

De dekpuntstelling vertelt ons wanneer een operator van een metrische ruimte naar zichzelf een uniek dekpunt heeft (zie [WISA4, hfst. 1]). We komen op de dekpuntstelling verder meer uitgebreid terug, en merken hier alleen op dat het unieke dekpunt kan gevonden worden door te vertrekken van een punt x0 en daarop successievelijk de operator O toe te passen. In ons geval is O een gewone numerieke functie, en stellen als algoritme voor xn+1 = xn + k(xn ) f (xn )

(2.3)

Als de rij (xn ) convergeert, dan is de limiet x zeker een oplossing van (2.2) (en (2.1)). Het volstaat in (2.3) de limiet te nemen voor n → ∞ van beide leden. De methode van Newton Onderstel dat f differentieerbaar is over [a, b]. Onderstel dat we een benadering xn van een wortel van (2.1) kennen. We benaderen de kromme y = f (x) door de raaklijn in het punt (xn , f (xn )). We nemen als nieuwe benadering xn+1 van de oplossing van (2.1) het snijpunt van deze raaklijn en de x-as. y

y = f(x) f(x)n

xn+1

xn

x

x

Figuur 2.1: De methode van Newton De vergelijking van de raaklijn in (xn , f (xn )) is y − f (xn ) = f 0 (xn )(x − xn ) Stel y = 0, en los op naar x. Dit geeft de formule van Newton xn+1 = xn −

f (xn ) f 0 (xn )

(2.4)

(2.4) is het speciaal geval van (2.3), waarbij k(x) = −1/ f 0 (x). De methode van Newton is de belangrijkste methode om algebraische vergelijkingen op te lossen. We zullen zien dat de methode relatief snel convergeert naar een oplossing in zeer veel gevallen. Het nadeel is dat men bij elke iteratieslag de afgeleide f 0 (x) moet berekenen. Daarom gebruikt men soms de volgende vereenvoudiging. 14

De methode van de vaste richting Om te vermijden dat men de afgeleide bij elke iteratieslag moet berekenen vervangt men f 0 (xn ) door f 0 (x0 ). Men krijgt dan de methode van de vaste richting. xn+1 = xn −

f (xn ) f 0 (x0 )

(2.5)

(2.5) is weer een speciaal geval van (2.3) De methode van de koorde Dit is een iteratieve methode waarbij xn+1 berekend wordt uit xn en xn−1 . Men vervangt de kromme y = f (x) door de koorde die (xn−1 , f (xn−1 )) en (xn , f (xn )) verbindt (zie Figuur ??). y y = f(x)

x xn–1

xn+1

x xn

Figuur 2.2: De methode van de koorde De vergelijking van deze koorde is y − f (xn ) f (xn−1 ) − f (xn ) = x − xn xn−1 − xn Het snijpunt xn+1 van de koorde met de x-as wordt gegeven door de formule xn+1 =

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

(2.6)

De methode van de koorde is de belangrijkste methode om vergelijkingen op te lossen na die van Newton. Een van de voordelen is dat men geen afgeleiden hoeft te berekenen.

15

Regula Falsi Dit is een verfijning van de methode van de koorde. Onderstel bij de methode van de koorde dat xn+1 berekend werd uitgaand van xn en xn−1 . Bij het berekenen van xn+2 kan men dan ofwel xn+1 en xn gebruiken, ofwel xn+1 en xn−1 . Bij de regula falsi (methode van de valse positie van het nulpunt) gaat men als volgt tewerk. Onderstel dat x0 en x1 kunnen gevonden worden zodat f (x0 ) en f (x1 ) verschillend teken hebben. Men bepaalt x2 met behulp van de methode van de koorde, en berekent f (x2 ). Als f (x0 ) en f (x2 ) verschillend teken hebben, dan werkt men verder met x0 en x2 . Als f (x1 ) en f (x2 ) verschillend teken hebben, dan werkt men verder met x1 en x2 . In sommige gevallen kan men een van de uiteinden van het interval kiezen. Dit doet zich bijvoorbeeld voor indien f (x0 ) > 0, f (x1 ) < 0 en f 00 > 0 ( f is concaaf). Men spreekt dan van de methode van Lin. Oplossing door dichotomie Men gaat te werk zoals bij de regula falsi, maar in plaats van de methode van de koorde toe te passen, deelt men gewoon het interval in twee. Onderstel dat x0 en x1 kunnen gevonden worden zodat f (x0 ) en f (x1 ) verschillend teken hebben. Stel x2 = (x0 + x1 )/2, en werk dan voort met x0 en x2 of met x1 en x2 alnaargelang het teken van f (x2 ). Een voordeel van de methode is dat men steeds een bovengrens voor de fout kent: de fout wordt bij elke iteratie door twee gedeeld. Een nadeel is dat de convergentie nogal traag verloopt: men heeft drie tot vier iteratieslagen nodig om een beduidend cijfer meer van de oplossing te kennen. √ √ Voorbeeld 2.1.1 Gegeven is een positief getal a. Er wordt gevraagd om a te berekenen. a is een wortel van de vergelijking x2 − a = 0 De methode van Newton levert als algoritme a xn2 − a 1  = xn + xn+1 = xn − 2xn 2 xn We beweren dat lim xn = n→∞

√

(2.7)

a als x0 > 0. Om dit in te zien tonen we eerst aan dat √ n √ n √ (x0 + a)2 + (x0 − a)2 √ n √ n xn = a (x0 + a)2 − (x0 − a)2

(2.8)

Dit gaat per inductie op n. Voor n = 0 is de formule duidelijk. Onderstel dat (2.8) waar is voor n. Dan is 1 a xn+1 = xn + 2 xn √  √ n √ n √ n √ n a (x0 + a)2 + (x0 − a)2 (x0 + a)2 − (x0 − a)2 √ n √ n+ √ n √ n = 2 (x0 + a)2 − (x0 − a)2 (x0 + a)2 + (x0 − a)2 16

√ n √ n
2 √ n
2 √ n a (x0 + a)2 + (x0 − a)2 + (x0 + a)2 − (x0 − a)2 = √ n+1 √ n+1 2 (x0 + a)2 − (x0 − a)2 √ n+1 √ n+1 √ (x0 + a)2 + (x0 − a)2 a = √ n+1 √ n+1 (x0 + a)2 − (x0 − a)2 √

zoals gewenst. Stel nu

√ x0 − a √ y= x0 + a √ √ Als x0 > a, dan is 0 < y < 1, en als 0 < x0 < a, dan is −1 < y < 0. In elk geval vinden we dat lim xn =

n→∞

We berekenen als voorbeeld

√

√

n

√ 1 + y2 a lim a n = 2 n→∞ 1 − y

3. Een mogelijk MATLAB programma is het volgende: √ % berekening van a a = 3; x(1) = 2; f or i = 1 : 10 x(i + 1) = (x(i) + (a/x(i)))/2; end x

Startwaarde was x(1) = 2; na drie iteraties vinden we reeds x(4) = 1.73205080756888; alle decimalen zijn correct. Voorbeeld 2.1.2 Men kan√de methode van Newton ook gebruiken om andere algoritmen op te schrijven die toelaten om a te berekenen. Het volstaat de vergelijking onder een andere vorm te √ herschrijven. a is een oplossing van de vergelijking a 1− 2 = 0 x De methode van Newton levert het volgende algoritme 1 (3axn − xn3 ) (2.9) 2a √ √ We beweren dat dit algoritme convergeert naar a als 0 < x < a. √ Om dit in √te zien merken we eerst op dat, indien (xn ) convergeert naar x, dat dan x gelijk is aan 0, a of − a. Het volstaat om 1 3 de limiet √van (2.9) √ voor n → ∞ te nemen. We krijgen x = 2a (3ax − x ), en de wortels hiervan zijn juist 0, a en − a. We bestuderen nu eerst het teken van 1 xn+1 − xn = (axn − xn3 ) 2a √ xn √ = ( a − xn )( a + xn ) 2a xn+1 =

17

√ Als 0 < xn < a, dan is dus xn+1 > xn . Bekijk nu de functie f (x) = (3ax − x3 )/2a. De afgeleide hiervan is f 0 (x) =

√ 1 3 √ (3a − 3x2 ) = ( a − x)( a + x) 2a 2a

√ Voor x ∈ (0, a), is f 0 (x) > 0. We hebben dus een monotoon stijgende bijectie √ √ f : (0, a)−→(0, a) √ Als x0 ∈ (0, a), dan geldt dus dat √ 0 < x0 < x1 < x2 < x3 < · · · < a √ (xn ) is stijgend en begrensd, en dus convergent. De limiet kan enkel a zijn. √ √ 0 Voor x ∈ ( a, 3a), is f (x) < 0. We hebben dus een monotoon dalende bijectie √ √ √ f : ( a, 3a)−→(0, a) √ √ √ Als x0 ∈ ( a, 3a), dan is x1 ∈ (0, a) en √ 0 < x1 < x2 < x3 < · · · < a √ en (xn ) convergeert naar a. √ Bij wijze van voorbeeld berekenen we weer 3. Een eenvoudig MATLAB programma is √ % berekening van a a = 3; x(1) = 2; f or i = 1 : 10 x(i + 1) = (3 ∗ a ∗ x(i) − x(i)3 )/(2 ∗ a); end x We namen als startwaarde x(1) = 2, zoals in het voorgaande voorbeeld. Na vijf iteraties vinden we de correcte waarde x(6) = 1.73205080756888. √ √ Ga zelf na dat het algorimte convergeert naar − a als x0 ∈ (− 3a, 0). De methode van Newton voor een veelterm Gegeven is een re¨ele veelterm f (x), en gevraagd wordt om deze in factoren te ontwikkelen. We weten dat de lineaire factoren overeenstemmen met de nulpunten van de veelterm. Deze kunnen we bepalen met behulp van de methode van Newton. Bij elke iteratieslag moeten we f (xn ) en f 0 (xn ) berekenen. Uiteindelijk moeten we ook de co¨effici¨enten 18

van het quoti¨ent bepalen. Ook hiervoor trachten we bij elke iteratieslag benaderingen te vinden. Schrijf f (x) = a0 xN + a1 xN−1 + · · · + aN en neem een re¨eel getal a. Rest en quoti¨ent bij deling van f (x) door x − a kunnen bepaald worden met behulp van de regel van Horner. Laten we deze (opnieuw) opstellen. Delen van f (x) door x − a geeft f (x) = (x − a)q(x) + f (a) (2.10) Schrijf q(x) = b0 xN−1 + b1 xN−2 + · · · + bN−1 Gelijkstellen van co¨effici¨enten van gelijke machten in (2.10) levert  a0 = b0     a1 = −ab0 + b1     .. . ak = −abk−1 + bk     .  ..    aN = −abN−1 + f (a) Stel bN = f (a). b0 , . . . , bN kunnen berekend worden met de formules b0 = a0 en bk = abk−1 + ak

(2.11)

voor k = 1, . . . , N. Dit is inderdaad niets anders dan de regel van Horner. Om f 0 (a) te bepalen leiden we (2.10) af naar x. f 0 (x) = q(x) + (x − a)q0 (x) en hieruit volgt dat f 0 (a) = q(a). Men berekent q(a) opnieuw met de regel van Horner, maar nu toegepast op de veelterm q. Dit geeft de volgende formules. c0 = b0 en ck = ack−1 + bk

(2.12)

voor k = 1, . . . , N − 1. Nu is cN−1 = q(a) = f 0 (a). Bij elke iteratieslag gaat men nu als volgt te werk: onderstel dat men een benadering xn voor een wortel heeft. Men bepaalt f (xn ), f 0 (xn ) en de co¨effici¨enten van het quoti¨ent qn (x) zoals hierboven (met xn = a), en men past dan de formule van Newton toe om xn+1 te berekenen. We kunnen het algoritme als volgt samenvatten: 1) Onderstel dat we een benadering xn voor a hebben; 2) stel b0 = c0 = a0 ; 3) bereken b1 , b2 , . . . , bN met de recursieformule (2.11): bk = xn bk−1 + ak 19

4) bereken c1 , c2 , . . . , cN−1 met de recursieformule (2.12): ck = xn ck−1 + bk 5) bepaal xn+1 met behulp van de formule van Newton: xn+1 = xn −

2.2

bn cN−1

Convergentie

De dekpuntstelling Onder welke voorwaarden convergeert een iteratieve methode? Om dit algemeen te kunnen behandelen herhalen we eerst de contractiestelling , die een speciaal geval is van de dekpuntstelling. Voor meer details verwijzen we naar volume 4 van de cursus “Wiskundige Analyse”. Een metrische ruimte is een verzameling E uitgerust met een afstand d. Een afstand d is een afbeelding d : E × E−→R+ die voldoet aan de volgende axiomas. 1) d(x, y) = d(y, x); 2) d(x, y) = 0 ⇐⇒ x = y; 3) d(x, y) ≤ d(x, z) + d(z, y), voor elke x, y, z ∈ E. Een rij (xn ) in een metrische ruimte E convergeert naar x ∈ E als ∀ε > 0, ∃N : n > N =⇒ d(x, xn ) < ε We noemen (xn ) een Cauchy rij als ∀ε > 0, ∃N : n, m > N =⇒ d(xm , xn ) < ε Men bewijst gemakkelijk dat elke convergente rij een Cauchy rij is. Indien elke Cauchy rij convergent is, dan zeggen we dat de metrische ruimte E volledig is. Een afbeelding O : E → E (soms ook operator genoemd) is een contractie of een samentrekking indien er een c ∈ [0, 1) bestaat zodat d(O(x), O(y)) ≤ cd(x, y) voor elke x, y ∈ E. De contractiestelling kan nu als volgt geformuleerd worden. Stelling 2.2.1 Onderstel dat E, d een volledige metrische ruimte is, en dat O : E → E een contractie is. Dan hebben we volgende eigenschappen. 1) O heeft juist e´ e´ n dekpunt x. Dit betekent dat er juist e´ e´ n x ∈ E bestaat die een oplossing is van de vergelijking O(x) = x 20

2) Voor elke x ∈ E geldt: lim On (x) = x

n→∞

3) We hebben bovendien: d(x, On (x)) ≤

cn d(x, O(x)) 1−c

Convergentie van de algemene formule Beschouw de vergelijking f (x) = 0, en onderstel dat k : [a, b] → R een continue functie is zonder nulpunten in [a, b]. We zagen in vorige paragraaf de volgende iteratieve formule xn+1 = xn + k(xn ) f (xn ) = g(xn )

(2.13)

waarbij g(x) = x + k(x) f (x). We zullen altijd impliciet onderstellen dat g continu differentieerbaar is. We zouden graag de contractiestelling toepassen in de volgende situatie: E = [a, b], d(x, y) = |x−y| en O = g. Merk op dat [a, b] volledig is als metrische ruimte. Dit volgt uit het feit dat R volledig is, en het feit dat het beschouwde interval gesloten is. Stelling 2.2.2 Onderstel dat voldaan is aan de volgende twee voorwaarden: 1) g(x) ∈ [a, b], voor elke x ∈ [a, b]; 2) Er bestaat een c ∈ [0, 1) zodat |g0 (x)| ≤ c voor elke x ∈ (a, b); dan bezit de vergelijking f (x) = 0 juist e´ e´ n wortel in [a, b]. Deze kan gevonden worden als de limiet van de rij (xn ), waarbij x0 een willekeurige startwaarde in [a, b] is. Bewijs. De eerste voorwaarde garandeert dat g : E → E een operator is. Uit de stelling van Lagrange volgt dat |g(x) − g(y)| = |g0 (ξ)(x − y)| ≤ c|x − y| voor elke x, y ∈ [a, b]. Hierin is ξ ∈ (x, y) ⊂ (a, b), zodat |g0 (ξ)| ≤ c. g is dus een samentrekking, en we kunnen de contractiestelling toepassen.  Opmerkingen 2.2.3 1) De eerste voorwaarde in de stelling is in het algemeen degene die het moeilijkste te verifi¨eren is. Een goede keuze van het interval [a, b] is hier cruciaal. [a, b] mag niet te groot zijn: het kan zijn dat f (x) = 0 meer dan e´ e´ n oplossing heeft, en in dit geval moet het interval [a, b] zo gekozen worden dat er maar e´ e´ n wortel in ligt. 2) De contractiestelling geeft ons ook een bovengrens voor de fout. |xn − x| ≤

cn |x0 − x1 | 1−c

Hoe dichter de startwaarde x0 bij de gezochte wortel ligt, hoe kleiner |x0 − x1 | is, en hoe beter de convergentie zal zijn. Hoe dichter c bij 1 ligt, hoe trager de convergentie zal verlopen; dit wordt veroorzaakt door twee factoren: cn in de teller gaat trager naar nul, en de noemer 1 − c ligt dichter bij nul. 21

Als het interval symmetrisch is rond de gezochte wortel, dan kunnen we voorwaarde 1) uit de voorgaande stelling elimineren. Stelling 2.2.4 Onderstel dat x een oplossing is van de vergelijking f (x) = 0. Als er een c ∈ [0, 1) bestaat zodat |g0 (x)| ≤ c voor elke x ∈ (x − δ, x + δ), dan is voldaan aan de voorwaarden van stelling 2.2.2, en we kunnen besluiten dat de rij (xn ) naar x convergeert. Bewijs. We moeten alleen voorwaarde 1) uit stelling 2.2.2 verifi¨eren. Neem x ∈ (x − δ, x + δ). Dan is |g(x) − x| = |g(x) − g(x)| = |g0 (ξ)||x − x| ≤ c|x − x| < δ zodat g(x) ∈ (x − δ, x + δ). g is (links)continu in x − δ, zodat g(x − δ) = lim g(x − δ + h) ∈ [x − δ, x + δ] h→0+



en hetzelfde geldt voor het andere eindpunt.

Gevolg 2.2.5 Als x een wortel is van de vergelijking f (x) = 0, g0 continu is in x, en |g0 (x)| < 1, dan convergeeert de iteratieve methode (2.13) zodra de startwaarde voldoende dicht bij x ligt. Bewijs. g is continu differentieerbaar in x. Stel c = |g0 (x)|, en neem ε > 0 zodat c + ε < 1. Er bestaat een δ > 0 zodat |x − x| < δ =⇒ |g0 (x) − g0 (x)| < ε =⇒ −c − ε ≤ g0 (x) − ε < g0 (x) < g0 (x) + ε ≤ c + ε =⇒ |g0 (x)| < c + ε 

We kunnen nu stelling 2.2.4 toepassen. Convergentie van de methode van Newton In dit bijzonder geval is g(x) = x − en g0 (x) = 1 −

f (x) f 0 (x)

f (x) f 00 (x) f 0 (x) f (x) f 00 (x) = + f 0 (x) f 0 (x)2 f 0 (x)2

Onderstel dat f 0 (x) 6= 0. Dan is g0 (x) = 0 zodat aan de voorwaarden van gevolg 2.2.5 voldaan is. We hebben dus volgende eigenschap bewezen. Stelling 2.2.6 Als x een nulpunt is van de vergelijking f (x) = 0, en f 0 (x) 6= 0, dan convergeert de methode van Newton als de startwaarde x0 voldoende dicht bij x gekozen wordt. 22

Opmerkingen 2.2.7 1) We zien ook dat lim g0 (x) = g0 (x) = 0

x→x

Dit betekent het volgende: hoe dichter de benadering xn bij x, hoe kleiner g0 (x), en hoe kleiner c. Dit betekent dat de convergentie sneller en sneller verloopt naarmate xn dichter bij de wortel x komt. Men zegt dat de methode van Newton zelfverbeterend is. 2) De methode van Newton werkt niet als f 0 (x) = 0. Dit betekent in feite dat x geen enkelvoudige wortel van de vergelijking is. Als men van tevoren weet dat er een dubbele wortel x bestaat (d.w.z. dat f (x) = (x − x)2 f1 (x), met f1 (x) 6= 0), dan kan men de vergelijking f (x) = 0 vervangen door p | f (x)| = 0 Deze vergelijking heeft dan een enkelvoudig nulpunt. Toepassing van de methode van Newton levert volgend algoritme. f (xn ) xn+1 = xn − 2 0 f (xn ) Convergentie van de methode van de vaste richting Nu is g(x) = x −

f (x) f 0 (x0 )

g0 (x) = 1 −

f 0 (x) f 0 (x0 )

en

De voorwaarde |g0 (x)| ≤ c < 1 wordt 1−c ≤

f 0 (x) ≤ 1+c f 0 (x0 )

met c ∈ [0, 1). Dit betekent dat de richting van de raaklijn aan de kromme y = f (x) ongeveer constant is over het beschouwde interval.

2.3

De orde van een iteratieve methode

Onderstel dat x een wortel is van de vergelijking f (x) = 0, en dat (xn ) een iteratieve rij is, waarvan we verwachten dat ze naar x convergeert. De fout die we maken bij de n-de iteratie is dan xn − x. Onderstel dat we een rij bovengrenzen (mn ) voor de fout kunnen bepalen: |xn − x| ≤ mn

23

Definitie 2.3.1 Als er een rij bovengrenzen (mn ) bestaat, en constanten k en α zodat mn+1 = kmαn dan zeggen we dat de iteratieve methode van orde α is. Stelling 2.3.2 Als α ≥ 1 en kmα−1 n0 < 1 voor een zekere index n0 , dan is de iteratieve methode convergent. Bewijs. De convergentie van de rij (xn ) verandert niet als we de eerste n0 termen van de rij weglaten. We kunnen dus gerust onderstellen dat n0 = 0. We stellen c = kmα−1 0 Dan is en

m1 = kmα0 = kmα−1 0 m0 = cm0 1+α m2 = kmα1 = kcα mα0 = kcα mα−1 m0 0 m0 = c

en

m3 = kmα2 = kcα+α mα0 = c1+α+α m0 2

2

In het algemeen vinden we mn = c1+α+α

2 +···+αn−1

m0

Bewijs de algemene formule zelf, met behulp van inductie op n. Omdat α ≥ 1 en 0 ≤ c < 1 is lim c1+α+α

2 +···+αn−1

n→∞

=0

en lim mn = 0

n→∞



zodat er convergentie is.

Een iteratieve methode is dus convergent als de orde minstens 1 is, en als de startwaarde voldoende dicht bij de wortel ligt. Indien er convergentie is, dan gaat de convergentie des te sneller naarmate de orde van de methode hoger is. Voorbeelden 2.3.3 1) De methode door dichotomie is van orde 1: we hebben gezien dat de fout bij elke iteratie door twee gedeeld wordt. 1 mn+1 = mn 2 2) De methode van Newton is van orde 2. Herinner dat g(x) = x − 24

f (x) f 0 (x)

Onderstel dat g00 (x) begrensd is in een omgeving van x. Dan is zeker f 0 (x) 6= 0. Met behulp van de formule van Taylor vinden we xn+1 − x = g(xn ) − g(x) = (xn − x)g0 (x) + (xn − x)2 = (xn − x)2

g00 (ξ) 2

g00 (ξ) 2

Neem een gesloten omgeving van x waarop g00 begrensd is. Over deze omgeving kunnen we schrijven: |g00 (x)| < k. Als |xn − x| ≤ mn , dan geldt dus dat k k |xn+1 − x| ≤ |xn − x|2 ≤ m2n 2 2 en we kunnen dus mn+1 = 2k m2n nemen. 3) De methode van de vaste richting is van orde 1. Het bewijs is analoog met dat voor de methode van Newton, maar men schrijft de Taylor ontwikkeling slechts tot op orde 1 op. Als de methode van de vaste richting convergeert, dan convergeert ze dus trager dan de methode van Newton. Dit is net wat we verwachten. √ 4) Men kan aantonen dat de methode van de koorde asymptotisch van orde α = (1+ 5)/2 ≈ 1.618 is. Dit betekent dat er majoranten mn voor de fout na de n-de iteratie bestaan zodat mn+1 lim =c n→∞ mα n Het bewijs maakt gebruik van de getallen van Fibonacci. 5) De volgende iteratieve methode is van orde 3: f 00 (xn ) f (xn )2 f (xn ) − xn+1 = g(xn ) = xn − 0 f (xn ) 2 f 0 (xn )3 Het volstaat om aan te tonen dat g0 (x) = g00 (x) = 0. Redeneer vervolgens zoals bij de methode van Newton, maar ontwikkel nu tot op orde 3. Het nadeel bij bovenstaande formule is dat ze ingewikkeld is. In de praktijk wordt ze weinig gebruikt.

2.4

Versnelling van de convergentie: de methode van Aitken

We beschouwen weer de vergelijking f (x) = 0, en we onderstellen dat x een wortel is. Beschouw een iteratieproces xn+1 = g(xn ) zodanig dat g0 continu is op een omgeving van x. Als 0 < |g0 (x)| < 1, dan convergeert het iteratieproces naar x als de startwaarde x0 dicht genoeg bij x ligt (zie gevolg 2.2.5), en bovendien is de methode van orde 1 (net zoals de methode van de vaste richting). Dit betekent dat de convergentie relatief traag verloopt. Onderstel nu dat x0 − x 6= 0 25

Aangezien xn+1 − x = g(xn ) − g(x) = g0 (ξn )(xn − x) volgt dat ook xn+1 − x 6= 0. We krijgen ook dat xn+1 − x = lim g0 (ξn ) = g0 (x) n→∞ xn − x n→∞ lim

want g0 is continu. In benadering hebben we dus xn+1 − x ≈ g0 (x)(xn − x) xn+2 − x ≈ g0 (x)(xn+1 − x)

(2.14) (2.15)

We kunnen dit benaderde stelsel oplossen naar g0 (x) en x, en we verwachten dat dit een betere benadering geeft van x dan xn . (2.15) aftrekken van (2.14) geeft xn+2 − xn+1 g0 (x) ≈ xn+1 − xn Uit (2.14) volgt nu xn+1 − g0 (x)xn ≈ x(1 − g0 (x)) en xn+1 − g0 (x)xn 1 − g0 (x) xn+1 − xn xn − g0 (x)xn = + 1 − g0 (x) 1 − g0 (x) xn+1 − xn = xn + 1 − g0 (x) xn+1 − xn ≈ xn + −xn+1 1 − xn+2 xn+1 −xn

x ≈

= xn −

(xn+1 − xn )2 = xn0 xn+2 − 2xn+1 + xn

We verwachten dat xn0 een betere benadering is van x dan xn . Dit wordt aangetoond door de volgende stelling. Stelling 2.4.1 Beschouw een convergente numerieke rij (xn ), met limiet x. Als xn 6= x voor elke index n, en xn+1 − x lim =a n→∞ xn − x met 0 < |a| < 1, dan convergeert de rij (xn0 ) sneller naar x dan (xn ). Hiermee bedoelen we dat xn0 − x =0 n→∞ xn − x lim

(2.16)

xn0 wordt gegeven door de formule xn0 = xn −

(xn+1 − xn )2 xn+2 − 2xn+1 + xn 26

(2.17)

Bewijs. We noteren fn = xn − x, en εn = fn+1 / fn − a. Dan is fn+1 = fn (a + εn ) en lim εn = 0 n→∞

We vinden nu dat xn+2 − 2xn+1 + xn = (xn+2 − x) − 2(xn+1 − x) + (xn − x) = fn+2 − 2 fn+1 + fn = (a + εn+1 )(a + εn ) fn − 2 fn (a + εn ) + fn   2 0 = fn (a − 1) + εn met lim ε0n = 0. n→∞

De noemer in (2.17) is nooit nul omdat |a| < 1 en fn 6= 0. Verder is xn+1 − xn = fn+1 − fn = fn (a − 1 + εn ) en xn0 − x = xn − x −

fn2 (a − 1 + εn )2
fn (a − 1)2 + ε0n

zodat   (a − 1 + εn )2 xn0 − x = lim 1 − lim n→∞ n→∞ xn − x (a − 1)2 + ε0n ε0 − 2(a − 1)εn + ε2n =0 = lim n n→∞ (a − 1)2 + ε0n 

en dit bewijst de stelling.

Opmerking 2.4.2 Deze methode om de convergentie te versnellen wordt de methode van Aitken genoemd. Men kan de convergentie opnieuw versnellen door op de rij (xn0 ) opnieuw de methode van Aitken toe te passen. √ Voorbeeld 2.4.3 Stel xn = g(xn−1 ) = 1 − xn−1 . Dan is √ −1 + 5 x= = .618033989 2 Merk op dat √ 1 1 1+ 5 0 g (x) = − √ =− =− 6= 0 2x 4 2 1−x zodat stelling 2.4.1 van toepassing is. We berekenen eerst xn met behulp van het MATLAB programma x(1) = .5 f or i = 2 : 20 x(i) = sqrt(1 − x(i − 1)); end x 27

We vinden x1 x2 x3 x4 x5 x6 x7 x10 x20

= = = = = = = = =

.5 .7071068 .5411961 .6773506 .5680223 .6572501 .5854484 .6349981 .62008333

We passen nu de Aitken methode toe en berekenen xn0 , met behulp van het MATLAB programma f or i = 1 : 18 y(i) = x(i) − ((x(i + 1) − x(i))2 )/(x(i + 2) − 2 ∗ x(i + 1) + x(i)); end y We vinden x10 x20 x30 x40 x50 x60 0 x18

= = = = = = =

.6149898 .6159796 .6167128 .6171526 .6174642 .6176566 .6180317

Pas Aitken nogmaals toe, met MATLAB programma f or i = 1 : 16 z(i) = y(i) − ((y(i + 1) − y(i))2 )/(y(i + 2) − 2 ∗ y(i + 1) + y(i)); end z We vinden x100 x200 x300 x900 00 x12 00 x16

= = = = = =

.6188085 .6178119 .6182211 .6180374 .6180336 .6180340 28

Hoofdstuk 3 Het oplossen van stelsels 3.1

Algoritmen voor het oplossen van stelsels

Beschouw een stelsel van n vergelijkingen met n onbekenden   f1 (x1 , x2 , . . . , xn ) = 0    f2 (x1 , x2 , . . . , xn ) = 0 ..  .    f (x , x , . . . , x ) = 0 n 1 2 n

(3.1)

We zullen dit soms verkort opschrijven in vectorvorm. Stellen we     f1 x1      f2   x2     F =  ..  en X =  ...   .    xn

fn dan kunnen we (3.1) herschrijven als F(X) = 0

(3.2)

We gaan tewerk als in het voorgaand hoofdstuk, en trachten een rij vectoren (Xm ) te construeren die naar de oplossing convergeert. (3.2) is equivalent met X = X + M(X)F(X)

(3.3)

indien M(X) een reguliere n × n-matrix is, voor elke waarde van X. Naar analogie met (2.2) stellen we volgende iteratieve formule voor Xm+1 = Xm + M(Xm )F(Xm )

29

(3.4)

De methode van Newton-Raphson Dit is het hoger dimensionale analogon van de methode van Newton. Onderstel dat we een bena(m) (m) dering Xm = (x1 , . . . , xn ) van een oplossing van (3.2) gevonden hebben. We benaderen nu de functies fi door hun Taylor ontwikkeling tot op orde 1: ∂ fi (m) (m) (m) (x1 , . . . , xn )(x j − x j ) ∂x j j=1 n

fi (x1 , . . . , xn ) ≈ fi (x1 , . . . , xn ) + ∑ (m)

(m)

of, in matrixvorm F(X) ≈ F(Xm ) + J(Xm )(X − Xm ) waarbij

 ∂ f1

(X)

 ∂x1  ∂ f2  (X)  J(X) =  ∂x1  ..  .  ∂ fn (X) ∂x1 Een benaderde oplossing van F(X) = 0 is dus

∂ f1 (X) · · · ∂x2 ∂ f2 (X) · · · ∂x2 .. . ∂ fn (X) · · · ∂x2

 ∂ f1 (X) ∂xn   ∂ f2 (X)   ∂xn   ..  .  ∂ fn (X) ∂xn

Xm+1 = Xm − J(Xm )−1 F(Xm )

(3.5)

Dit is het algoritme van Newton-Raphson. Merk op dat men bij elke iteratie een lineair stelsel moet oplossen. Met behulp van de methode van Cramer kunnen we de formules uitschrijven voor kleine waarden van n. Neem n = 2, en schrijf X = (x, y) en Xm = (xm , ym ). (3.5) wordt  ∂ f1 ∂ f2   (xm , ym ) − f2 (xm , ym ) (xm , ym ) f1 (xm , ym )   ∂y ∂y   xm+1 = xm −    ∂ f1 ∂ f2 ∂ f1 ∂ f2   (xm , ym ) (xm , ym ) − (xm , ym ) (xm , ym )   ∂x ∂y ∂y ∂x (3.6)   ∂ f ∂ f  1 2   f2 (xm , ym ) (xm , ym ) − f1 (xm , ym ) (xm , ym )   ∂x ∂x  y = y −  m m+1  ∂ f1 ∂ f2 ∂ f1 ∂ f2    (xm , ym ) (xm , ym ) − (xm , ym ) (xm , ym ) ∂x ∂y ∂y ∂x De methode van Morrey Dit is het analogon van de methode van de vaste richting. Xm+1 = Xm − J(X0 )−1 F(Xm )

30

(3.7)

De methode van de diagonaalterm J(X)−1 is gemakkelijk op te lossen als J(X) een diagonaalmatrix is. We stellen D(X) de matrix die ontstaat door in J(X) alle elementen die niet op de diagonaal staan gelijk aan nul te stellen   ∂ f1 (X) 0 ··· 0   ∂x1   ∂ f 2   0 (X) · · · 0   ∂x2 D(X) =     . . .   .. .. ..   ∂ fn (X) 0 0 ··· ∂xn en men benadert het stelsel F(X) = 0 door het stelsel F(X) ≈ F(Xm ) + D(Xm )(X − Xm ) = 0 De nieuwe iteratieve formule is nu Xm+1 = Xm − D(Xm )−1 F(Xm ) of, in componenten uitgeschreven (m)

(m+1) xi

(m) = xi −

(m)

fi (x1 , . . . , xn ) ∂ fi (m) (m) (x1 , . . . , xn ) ∂xi

(3.8)

(3.8) is natuurlijk eenvoudiger dan (3.5). Zoals we verder zullen zien is de convergentie van (3.8) minder goed dan die van (3.5).

3.2

Metrieken op Rn

We gebruiken hier dezelfde methode als in het voorgaande hoofdstuk. Met behulp van de dekpuntstelling trachten we aan te tonen dat onze iteratieve methode naar een oplossing convergeert. Onze oplossingen zijn nu vectoren in Rn , zodat we op zoek moeten gaan naar metrieken op Rn . We herhalen eerst even de definitie van een afstand of metriek op een verzameling. Definitie 3.2.1 Een metrische ruimte (E, d) is een verzameling E, tesamen met een afbeelding d : E × E → R+ : (x, y) 7→ d(x, y) die voldoet aan de volgende voorwaarden, voor elke x, y, z ∈ E: 1. d(x, y) = 0 als en alleen als x = y; 2. d(x, y) ≤ d(x, z) + d(z, y); 31

3. d(x, y) = d(y, x). We noemen d(x, y) de afstand van x tot y. Herhaal ook dat een norm op een vectorruimte als volgt gedefinieerd wordt. Definitie 3.2.2 Een re¨ele vectorruimte V die is uitgerust met een afbeelding k • k : V −→R+ die voldoet aan de volgende eigenschappen, voor alle ~x,~y ∈ V en α ∈ R: 1. k~xk = 0 ⇐⇒ ~x = ~0; 2. kα~xk = |α|k~xk; 3. k~x +~yk ≤ k~xk + k~yk. wordt een genormeerde ruimte genoemd. Een genormeerde ruimte is altijd een metrische ruimte. Eigenschap 3.2.3 Als V, k • k een genormeerde re¨ele vectorruimte is, dan definieert de formule d(~x,~y) = k~x −~yk een afstand op V . Elke genormeerde ruimte is dus een metrische ruimte. Op Rn kan men de volgende normen - en hun bijhorende metrieken - bekijken. q k~xk1 = x12 + x22 + · · · + xn2

(3.9)

Deze norm noemt men de Euclidische norm op Rn . De bijhorende afstand noteren we d1 . k~xk2 = max{|x1 |, |x2 |, . . . , |xn |}

(3.10)

Deze norm noemt men de maximum norm op Rn . De bijhorende afstand noteren we d2 . k~xk3 = |x1 | + |x2 | + · · · + |xn |

(3.11)

Op Rn kan men dus verschillende normen (en dus verschillende afstanden) defini¨eren. Er is een verband tussen deze drie afstanden. We hebben namelijk k~xk2 = max{|x1 |, |x2 |, . . . , |xn |} q x12 + x22 + · · · + xn2 = k~xk1 ≤ q ≤ (|x1 | + |x2 | + · · · + |xn |)2 = k~xk3 ≤ n max{|x1 |, |x2 |, . . . , |xn |} = nk~xk2 32

Samengevat: k~xk2 ≤ k~xk1 ≤ k~xk3 ≤ nk~xk2

(3.12)

voor alle ~x ∈ Rn . Voor n = 1 vallen de drie normen samen. Voor de convergentie speelt het geen rol met welk van de drie afstanden we werken. Stelling 3.2.4 Onderstel dat (~xn ) een rij vectoren in Rn is. De volgende uitspraken zijn equivalent: 1. d3 - lim ~xn = ~a; n→∞

2. d1 - lim ~xn = ~a; n→∞

3. d2 - lim ~xn = ~a. n→∞

Hierbij zijn d1 , d2 en d3 de afstanden die behoren bij de normen k • k1 , k • k2 en k • k3 . Een analoge eigenschap geldt voor Cauchyrijen. Bewijs. 1. =⇒ 2. Uit 1. volgt ∀ε > 0, ∃N : n > N =⇒ k~xn −~ak3 < ε Uit (3.12) volgt dat k~xn −~ak1 ≤ k~xn −~ak3 en dus

∀ε > 0, ∃N : n > N =⇒ k~xn −~ak1 < ε

en 2. volgt. Op juist dezelfde manier bewijzen we dat 2. =⇒ 3. =⇒ 1..



Het groot voordeel van stelling 3.2.4 is dat we, om de convergentie te onderzoeken, we de “gemakkelijkste”van de drie normen kunnen kiezen. In de praktijk is dit bijna nooit de Euclidische norm, maar wel dikwijls de maximumnorm of de norm k • k3 . Stelling 3.2.5 Onderstel dat G een gesloten deel van Rn is. Dan is (G, di ) een volledige metrische ruimte (i = 1, 2, 3). Bewijs. Uit stelling 3.2.4 blijkt dat het voldoende is om de stelling te bewijzen voor een van de drie hierboven ingevoerde metrieken. We werken met de metriek d2 . Ondersel dat (Xm ) een d2 -Cauchy rij in G is. ∀ε > 0, ∃N ∈ N : p, m > N =⇒ kX p − Xm k2 < ε (p)

(m)

(m)

Omdat |xi − xi | ≤ kX p − Xm k2 volgt hieruit dat (xi ) een numerieke Cauchyrij is, voor elke i = 1, . . . , n. We weten dat elke numerieke Cauchyrij convergeert, en stellen (m)

xi = lim xi m→∞

33

en X = (x1 , x2 , . . . , xn ). Voor elke ε > 0 hebben we (m)

∃Ni ∈ N : m > Ni =⇒ |xi

− xi | < ε

Stel N = max{N1 , N2 , . . . , Nn }. Dan geldt voor elke m > N dat (m)

kXm − Xk2 = max{|xi

− xi | : i = 1, . . . , n} ≤ ε

en d2 - lim Xm = X m→∞

Omdat G gesloten is, hebben we dat X ∈ G. Dus is G volledig.

3.3



Convergentie van de iteratieve formules

Om de notaties wat te vereenvoudigen werken we in R3 . Alles kan gemakkelijk uitgebreid worden tot Rn . Beschouw een functie O : G ⊂ R3 → R3 . Neem (x, y, z) ∈ G, en schrijf (u, v, w) = O(x, y, z). We zeggen dat O voldoet aan de voorwaarde (V ) indien er αi , βi , γi ∈ R3 bestaan zodat   |u2 − u1 | ≤ α1 |x2 − x1 | + β1 |y2 − y1 | + γ1 |z2 − z1 | |v2 − v1 | ≤ α2 |x2 − x1 | + β2 |y2 − y1 | + γ2 |z2 − z1 | (3.13)  |w2 − w1 | ≤ α3 |x2 − x1 | + β3 |y2 − y1 | + γ3 |z2 − z1 | voor alle (x1 , y1 , z1 ), (x2 , y2 , z2 ) ∈ G. We zeggen dat O voldoet aan het kolomcriterium als O voldoet aan de voorwaarde (V ), met   α1 + α2 + α3 ≤ θ < 1 β1 + β2 + β3 ≤ θ < 1 (3.14)  γ1 + γ2 + γ3 ≤ θ < 1 We zeggen dat O voldoet aan het rijcriterium als O voldoet aan de voorwaarde (V ), met αi + β i + γi ≤ θ < 1

(3.15)

for i = 1, 2, 3. Stelling 3.3.1 Als O voldoet aan het rijcriterium, dan is O een contractie voor de metriek d2 . Als O voldoet aan het kolomcriterium, dan is O een contractie voor de metriek d3 . Bewijs. We schrijven P1 = (x1 , y1 , z1 ) ; P2 = (x2 , y2 , z2 ) O(P1 ) = (u1 , v1 , w1 ) ; O(P2 ) = (u2 , v2 , w2 ) Onderstel O voldoet aan het rijcriterium. We hebben bijvoorbeeld dat max{|u2 − u1 |, |v2 − v1 |, |w2 − w1 |} = |v2 − v1 | 34

Dan is d2 (O(P1 ), O(P2 )) = = ≤ ≤ ≤

max{|u2 − u1 |, |v2 − v1 |, |w2 − w1 |} |v2 − v1 | α2 |x2 − x1 | + β2 |y2 − y1 | + γ2 |z2 − z1 | (α2 + β2 + γ2 )d2 (P1 , P2 ) θd2 (P1 , P2 )

zodat O een contractie is voor de metriek d2 . Onderstel vervolgens dat O voldoet aan het kolomcriterium. d3 (O(P1 ), O(P2 )) = |u2 − u1 | + |v2 − v1 | + |w2 − w1 | ≤ (α1 + α2 + α3 )|x2 − x1 | +(β1 + β2 + β3 )|y2 − y1 | +(γ1 + γ2 + γ3 )|z2 − z1 | ≤ θd3 (P1 , P2 ) 

zodat O een contractie is voor de metriek d3 . Convergentie van de algemene iteratieve formule We keren terug naar het stelsel F(X) = 0, en de iteratieve formule Xm+1 = Xm + M(Xm )F(Xm ) = O(Xm )

(3.16)

Om de notaties te vereenvoudigen werken we in R3 . Voor X ∈ R3 schrijven we     u φ1 (x, y, z)     O(X) = U =  v  =  φ2 (x, y, z)  w

φ3 (x, y, z)

We geven nu twee convergentiestellingen die alletwee een uitbreiding zijn van stelling 2.2.4. We ∂φi continu zijn. zullen in het vervolg steeds onderstellen dat de parti¨ele afgeleiden ∂x j Stelling 3.3.2 Onderstel dat X een oplossing is van het stelsel F(X) = 0, en dat er αi , βi , γi bestaan zodat ∂φ i (x, y, z) ≤ αi ∂x ∂φ i (3.17) (x, y, z) ≤ βi ∂y ∂φ i (x, y, z) ≤ γi ∂z voor i = 1, 2, 3 en voor alle (x, y, z) ∈ G2 = {X ∈ R3 : d2 (X, X) ≤ δ}. Als αi + βi + γi ≤ θ < 1 dan convergeert het iteratieproces naar X voor elke startwaarde X0 ∈ G2 . 35

(3.18)

Bewijs. We nemen (x1 , y1 , z1 ), (x2 , y2 , z2 ) ∈ G2 , en schrijven O(x j , y j , z j ) = (u j , v j , w j ) voor j = 1, 2. Vanwege de middelwaardestelling bestaat er een (ξ, η, ζ) ∈ [(x1 , y1 , z1 ), (x2 , y2 , z2 )] ⊂ G2 (G2 is convex) zodat |u2 − u1 | = |φ1 (x2 , y2 , z2 ) − φ1 (x1 , y1 , z1 )| ∂φ1 ∂φ ∂φ 1 1 = (ξ, η, ζ)(x2 − x1 ) + (ξ, η, ζ)(y2 − y1 ) + (ξ, η, ζ)(z2 − z1 ) ∂x ∂y ∂z ≤ α1 |x2 − x1 | + β1 |y2 − y1 | + γ1 |z2 − z1 | Dit is de eerste regel uit voorwaarde (V). Op analoge wijze zijn ook de twee andere regels voldaan, en uit (3.18) volgt dat O voldoet aan het rijcriterium. O is dus een contractie voor de metriek d2 . Bovendien is O : G2 → G2 . Immers, als X ∈ G2 , dan is d2 (O(X), X) = d2 (O(X), O(X)) ≤ θd(X, X) ≤ θδ < δ en O(X) ∈ G2 . We kunnen nu de contractiestelling toepassen, en concluderen dat (On (X0 )) naar X convergeert voor elke X0 ∈ G2 .  Stelling 3.3.3 Onderstel dat X een oplossing is van het stelsel F(X) = 0, en dat er αi , βi , γi bestaan zodat ∂φ i (x, y, z) ≤ αi ∂x ∂φ i (x, y, z) ≤ βi ∂y ∂φ i (x, y, z) ≤ γi ∂z voor i = 1, 2, 3 en voor alle (x, y, z) ∈ G3 = {X ∈ R3 : d3 (X, X) ≤ δ}. Als   α1 + α2 + α3 ≤ θ < 1 β1 + β2 + β3 ≤ θ < 1  γ1 + γ2 + γ3 ≤ θ < 1

(3.19)

dan convergeert het iteratieproces naar X voor elke startwaarde X0 ∈ G3 . Bewijs. Volledig analoog aan dat van stelling 3.3.2. We zeggen dat een matrix A ∈ Mnn (R) voldoet aan het rijcriterium als n

∑ |ai j | < 1

j=1

36



voor elke rijindex i. Analoog voldoet A aan het kolomcriterium als n

∑ |ai j | < 1

i=1

voor elke kolomindex j. We noteren F voor de Jacobiaanse matrix van φ1 , . . . , φn :  ∂φ1 ∂φ1 (X) (X) · · · ∂x2  ∂x1  ∂φ2 ∂φ2  (X) (X) · · ·  ∂x2 F (X) =  ∂x1  .. ..  . .  ∂φn ∂φn (X) (X) · · · ∂x1 ∂x2

 ∂φ1 (X) ∂xn   ∂φ2 (X)   ∂xn   ..  .  ∂φn (X) ∂xn

Stelling 3.3.4 Onderstel dat X een wortel is van het stelsel F(X) = 0, en dat F (X) voldoet aan het rij- of kolomcriterium Dan convergeert de algemene iteratieve methode van zodra de startwaarde X0 voldoende dicht bij X ligt. Bewijs. Onderstel bijvoorbeeld dat F (X) voldoet aan het rijcriterium. Dit betekent dat ∂φi (X) + ∂φi (X) + ∂φi (X) < 1 ∂x ∂y ∂z

(3.20)

Omdat de parti¨ele afgeleiden van de φi continu zijn, geldt (3.20) over een gesloten d2 -omgeving G2 van X. Stelling 3.3.2 is van toepassing, en de iteratieve methode convergeert naar X zodra X0 ∈ G2 .  Convergentie van de methode van Newton-Raphson We voeren volgende notatie in: voor een vectorwaardige functie O(X) = (φ1 (X), φ2 (X), . . . , φn (X)) schrijven we dO(X) = (dφ1 , dφ2 , . . . , dφn ) = F (X)dX Analoog defini¨eren we voor een matrix A(X) de differentiaal dA(X). We kunnen dan gemakkelijk nagaan dat d(AU) = AdU + (dA)U Onderstel nu dat X een oplossing is van F(X) = 0, en dat J(X) een reguliere matrix is. De methode van Newton-Raphson gebruikt volgende operator: O(X) = X − J(X)−1 F(X) 37

en we vinden dO(X) = = = =

F (X)dX dX − d(J(X)−1 )F(X) − J(X)−1 dF(X) dX − d(J(X)−1 )F(X) − J(X)−1 J(X)dX −d(J(X)−1 )F(X)

Omdat F(X) nul is vinden we dat dO(X) = F (X)dX = 0 en F (X) = 0. F (X) voldoet dus zeker aan het rij- en kolomcriterium, en uit stelling 3.3.4 volgt dat de methode van Newton-Raphson convergeert als de startwaarde X0 dicht genoeg bij X ligt. Opmerking 3.3.5 Net zoals in het geval van e´ e´ n veranderlijke kan men aantonen dat de methode van Newton-Raphson voor een stelsel niet-lineaire vergelijkingen een methode van orde 2 is. Convergentie van de methode van de diagonaalterm We onderstellen weer dat F(X) = 0, en dat ∂ fi (X) 6= 0 ∂xi voor i = 1, 2, · · · , n. De iteratieve formule is nu φi (X) = xi −

We vinden

fi (X) ∂ fi (X) ∂xi

∂2 f i (X) ∂xi2 ∂φi (X) = 1 − 1 +  2 ∂xi ∂ fi (X) ∂xi fi (X)

en

Voor i 6= j vinden we

∂φi (X) = 0 ∂xi ∂ fi ∂2 f i (X) fi (X) (X) ∂x j ∂xi ∂x j ∂φi (X) = − +  2 ∂ fi ∂x j ∂ f i (X) (X) ∂xi ∂xi 38

en

∂ fi (X) ∂x j ∂φi (X) = − ∂ fi ∂x j (X) ∂xi

We besluiten dat de methode van de diagonaalterm naar X genoeg bij X kiezen, en als de matrix  ∂ f1 (X) ∂x  2 0  ∂ f1  (X)  ∂x  1   ∂ f2  (X)  ∂x1  0 ∂ f2 F (X) =   (X)  ∂x2   . .  .. ..  ∂ fn  ∂ fn  (X) (X) ∂x  ∂x1 2  ∂ f ∂ f n n (X) (X) ∂x ∂x n n

convergeert als we een startwaarde dicht  ∂ f1 (X) ∂x  n  · · · ∂ f1 (X)  ∂x  1  ∂ f2   ∂x (X)   · · · n  ∂ f2 (X)  ∂x   2   ..  .      ··· 0

voldoet aan het rij- of kolomcriterium. Voorbeeld 3.3.6 We bepalen de snijpunten van de ellipsen x2 y2 x2 y2 + = 1 en + =1 4 9 9 4 Er zijn vier snijpunten. Het snijpunt in het eerste kwadrant is 6 x = y = √ = 1.66410058867569 13 De methode van de diagonaalterm levert xn+1 =

xn 2y2n 2 − + 2 9xn xn

yn 2xn2 2 − + 2 9yn yn We nemen als startwaarden bijvoorbeeld x0 = y0 = 2. Een eenvoudig MATLAB programma yn+1 =

x(1) = 2; y(1) = 2; f or i = 1 : 20 39

x(i + 1) = x(i)/2 − (2 ∗ y(i)2 )/(9 ∗ x(i)) + 2/x(i); y(i + 1) = y(i)/2 − (2 ∗ x(i)2 )/(9 ∗ y(i)) + 2/y(i); end [x; y]0 levert 20 iteraties. We vinden x10 = y10 = 1.66417940061712 x15 = y15 = 1.66409922203013 x20 = y20 = 1.66410061237543 Opmerking 3.3.7 Het is natuurlijk onontbeerlijk om een goede startwaarde te hebben. Een mogelijke strategie is de volgende: men maakt een schets van de meetkundige situatie, en men tracht zo de ligging van de wortels te bepalen. In bovenstaand voorbeeld kan men de twee ellipsen eerst schetsen, en zo een ruwe benadering voor x0 en y0 bepalen. Een eenvoudig MATLAB programma is het volgende: t = 0 : .05 : 2 ∗ pi; x1 = 2 ∗ cos(t); y1 = 3 ∗ sin(t); x2 = 3 ∗ cos(t); y2 = 2 ∗ sin(t); axis(0 square0 ) plot(x1, y1, x2, y2) grid Voor stelsels met meer dan twee vergelijkingen en onbekenden wordt dit uiteraard ingewikkelder. Indien het gegeven stelsel een storing is van een eenvoudiger stelsel, dan kan men als startwaarde de oplossing van het eenvoudige stelsel proberen. Om de snijpunten van de ellipsen y2 x2 y2 x2 + = 1 en + =1 4.2 9.3 9.5 4.1 te bepalen, kunnen we vertrekken van de oplossingen van voorgaand voorbeeld.

3.4

De methode van Bairstow

Gegeven is een re¨ele veelterm f (x) = a0 xn + a1 xn−1 + · · · + an Gevraagd wordt om de veelterm te ontbinden in lineaire en kwadratische factoren. In Hoofdstuk 2 hebben we al gezien hoe we met behulp van de methode van Newton de lineaire factoren kan 40

bepalen. Om de kwadratische factoren te vinden hebben we de complexe nulpunten nodig. Deze zijn twee aan twee complex toegevoegd. Een eerste manier bestaat erin op te merken dat de methode van Newton ook werkt voor complexe analytische functies f : C → C. Men gebruikt dan als startwaarde bij de methode van Newton een complex getal. Deze methode is van toepassing als men beschikt over een rekenmachine waarin ook de bewerkingen met de complexe getallen geprogrammeerd zijn. We zullen hier een methode bespreken die enkel berekeningen in R gebruikt. Men deelt de veelterm f (x) door een kwadratische veelterm x2 + px + q. Dit geeft een quotient van graad n − 2 en een rest van graad 1: f (x) = (x2 + px + q)(b0 xn−2 + b1 xn−3 + · · · + bn−2 ) + Rx + S

(3.21)

De co¨effici¨enten b0 , b1 , . . . , bn−2 , R en S zijn in feite functies van p en q. x2 + px + q is een factor van f (x) als  R(p, q) = 0 (3.22) S(p, q) = 0 We moeten het stelsel (3.22) oplossen. We doen dit met de methode van Newton-Raphson. Om (3.6) te kunnen toepassen moeten we R, S,

∂S ∂R ∂R ∂S , , en ∂p ∂q ∂p ∂q

berekenen. Om R en S te bepalen voeren we gewoon de deling uit. Dit komt erop neer de regel van Horner te veralgemenen tot de situatie waarbij we delen door veeltermen van graad 2. Gelijkstellen van de co¨effici¨enten in (3.21) levert  a0 = b0     a1 = b1 + pb0     a2 = b2 + pb1 + qb0     .. . (3.23) ak = bk + pbk−1 + qbk−2     ..   .     a = R + pbn−2 + qbn−3   n−1 an = S + qbn−2 Stel nu b−1 = 0, bn−1 = R, bn = S − pbn−1 De co¨effici¨enten bk , en dus ook R en S worden dan gevonden uit de betrekkingen bk = ak − pbk−1 − qbk−2

(3.24)

met startwaarden b−1 = 0 en b0 = a0 . De twee laatste waarden bn−1 en bn leveren R en S. Om de parti¨ele afgeleiden van R en S naar p te bepalen leiden we (3.24) af naar p. Dit geeft ∂bk−1 ∂bk−2 ∂bk = −bk−1 − p −q ∂p ∂p ∂p 41

Stel nu ck−1 = −

∂bk ∂p

We krijgen ck−1 = bk−1 − pck−2 − qck−3 of ck = bk − pck−1 − qck−2 Omdat c−1 = −

(3.25)

∂a0 ∂b1 ∂b0 =− = 0 en c0 = = b0 = a0 ∂p ∂p ∂p

laat (3.25) ons toe om recursief c1 , c2 , . . . , cn−1 te berekenen. De getallen cn−2 en cn−1 leveren en

∂S want ∂p

∂R ∂p

∂bn−1 ∂R = = −cn−2 ∂p ∂p ∂S ∂bn ∂bn−1 = + bn−1 + p = bn−1 − cn−1 − pcn−2 ∂p ∂p ∂p De parti¨ele afgeleiden van R en S naar q worden op analoge manier berekend. We leiden (3.24) af naar q, en krijgen ∂bk−1 ∂bk−2 ∂bk = −p − bk−2 − q ∂q ∂q ∂q We stellen nu c0k−2 = −

∂bk ∂q

zodat c0k−2 = bk−2 − pc0k−3 − qc0k−4 of c0k = bk − pc0k−1 − qc0k−2 We merken op dat c0−1 = −

(3.26)

∂b1 ∂b2 = 0 = c−1 en c00 = − = a0 = c0 ∂q ∂q

De beginwaarden en de recursieformules ter berekening van ck en c0k zijn dezelfden, en we kunnen dus concluderen dat ck = c0k . Verder is ∂R ∂bn−1 = = −cn−3 ∂q ∂q ∂S ∂bn ∂bn−1 = +p = −cn−2 − pcn−3 ∂q ∂q ∂q 42

Als we de gevonden waarden voor R, S en hun parti¨ele afgeleiden vervangen in de formules van Newton-Raphson (3.6) ∂S ∂R R −S ∂q ∂q pm+1 = pm − ∂R ∂S ∂R ∂S − ∂p ∂q ∂q ∂p ∂R ∂S −R ∂p ∂p qm+1 = qm − ∂R ∂S ∂R ∂S − ∂p ∂q ∂q ∂p S

dan vinden we pm+1 = pm +

bn−1 cn−2 − bn cn−3 2 cn−2 + cn−3 (bn−1 − cn−1 )

(3.27)

qm+1 = qm +

bn cn−2 + bn−1 (bn−1 − cn−1 ) c2n−2 + cn−3 (bn−1 − cn−1 )

(3.28)

en we kunnen ons algoritme als volgt samenvatten. 1) Onderstel dat pm en qm bepaald zijn. 2) Stel b−1 = c−1 = 0 en b0 = c0 = a0 . 3) Bereken b1 , · · · , bn uit de recursieformule (3.24) bk = ak − pm bk−1 − qm bk−2 4) Bereken c1 , · · · , cn−1 uit de recursieformule (3.25) ck = bk − pm ck−1 − qm ck−2 5) Bepaal pm+1 en qm+1 met behulp van (3.27) en (3.28). Bij elke iteratie vindt men een betere benadering voor de kwadratische factor x2 + px + q, en ook voor het quoti¨ent b0 xn−2 + · · · + bn−2 .

43

Hoofdstuk 4 Lineaire stelsels Om een lineair stelsel vergelijkingen op te lossen zijn er twee klassen methodes. Vooreerst heeft men iteratieve methodes, waarbij men een rij kolomvectoren iteratief berekent, in de hoop dat deze naar de oplossing convergeren. Er zijn ook de directe methodes, zoals de eliminatiemethode van Gauss.

4.1

De iteratieve methode van Gauss

Dit is de methode van de diagonaalterm toegepast op een lineair stelsel. Beschouw, om de gedachten te vestigen, een 3 × 3-stelsel ai x + bi y + ci z = di (4.1) voor i = 1, 2, 3. Als we de methode van de diagonaalterm toepassen, dan krijgen we volgende iteratieve formules c1 d1 b1 − ym − zm a1 a1 a1 d2 a2 c2 = − xm − zm b2 b2 b2 d3 a3 b3 = − xm − ym c3 c3 c3

xm+1 =

(4.2)

ym+1

(4.3)

zm+1

(4.4)

Merk op dat we de formules gemakkelijk terugvinden door (4.1) op te lossen naar x, y en z. Een voldoende voorwaarde voor convergentie is dat de matrix  b1 c1  0 a1 a1   a c2   2 0 b2   b2 a3 b3 0 c3 c3 voldoet aan het rij- of kolomcriterium. 44

De variante van Seidel De computer werkt sequentieel en berekent eerst xm+1 en nadien ym+1 en zm+1 . Op het ogenblik dat ym+1 berekend wordt, is xm+1 reeds gekend. Aangezien we verwachten dat xm+1 een betere benadering voor x is dan xm , kunnen we in (4.3) xm beter vervangen door xm+1 . Om dezelfde reden vervangen we in (4.4) xm en ym door xm+1 en ym+1 . Dit levert de formules van Gauss-Seidel .

c1 d1 b1 − ym − zm a1 a1 a1 d2 a2 c2 = − xm+1 − zm b2 b2 b2 b3 d3 a3 − xm+1 − ym+1 = c3 c3 c3

xm+1 =

(4.5)

ym+1

(4.6)

zm+1

(4.7)

We verwachten dat formules van Gauss-Seidel sneller convergeren dan die van Gauss. In de praktijk blijkt dit ook het geval te zijn. Zonder bewijs vermelden we verder de volgende eigenschap. Stelling 4.1.1 Indien de matrix A van het stelsel AX = B symmetrisch en positief definiet is, dan convergeren de formules van Gauss-Seidel voor elke beginwaarde.

4.2

Permutatiematrices en Frobeniusmatrices

Zij A een n × m-matrix. We zullen soms volgende notatie gebruiken.   a11 a12 A13   A =  a21 a22 A23  A31

A32

A33

Hierbij zijn de Ai j matrices van de geschikte afmetingen: A33 is een (n − 2) × (m − 2)-matrix, A13 is een rijvector, enzovoort.

45

Permutatiematrices De n × n-matrix 

1 0 .  ..  i 0 Pi j =  .. .  j 0 .  .. 0

0 1 .. . 0 .. . 0 .. . 0

i 0 0 .. . 0 .. . 1 .. .

··· ··· ··· ···

··· 0

··· ··· ··· ···

j 0 0 .. . 1 .. . 0 .. .

··· 0

 ··· 0 ··· 0 ..  .  ··· 0 ..   .  ··· 0 ..  . ··· 1

noemen we een permutatiematrix. Pi j A is de matrix A met de i-de en j-de rij verwisseld. Voor i = 1 en j = 2 bijvoorbeeld, vinden we      0 1 0 a11 a12 A13 a21 a22 A23      P12 A =  1 0 0   a21 a22 A23  =  a11 a12 A13  0

0

I

A31

A32

A33

A31

A32

A33

Op dezelfde manier is APi j de matrix A met de i-de en de j-de kolom met elkaar verwisseld. Merk ook op dat Pi−1 j = Pi j . Een product van permutatiematrices noemen we nog steeds een permutatiematrix. Een matrix A linksvermenigvuldigen met een permuatiematrix komt er dus op neer de rijen van A te permuteren. Rechtsvermenigvuldigen komt erop neer de kolommen te permuteren. Frobeniusmatrices Een matrix van de vorm 

1 0 .  ..  0 i Fi =  i + 1 0  .. . 0

0 1 .. . 0 0 .. .

··· ···

0

···

i 0 0 .. . 1

··· · · · ai+1 .. .

46

an

 ··· 0 ··· 0 ..  .  ··· 0  ··· 0  ..  . ··· 1

noemen we een Frobeniusmatrix. De inverse van een Frobeniusmatrix is opnieuw Frobenius: 

1 0 .  ..  0 Fi−1 = i  i + 1 0  .. . 0

0 1 .. . 0 0 .. .

··· ···

0

···

i 0 0 .. . 1 −ai+1 .. .

··· ···

−an

 ··· 0 ··· 0 ..  .  ··· 0  ··· 0  ..  . ··· 1

Een Frobeniusmatrix is een benedendriehoeksmatrix, en dus is het product van Frobeniusmatrices opnieuw een benedendriehoeksmatrix. Als Fi een Frobeniusmatrix is, en j, k > i, dan is Pjk Fi = Fi0 Pjk

(4.8)

met Fi0 een nieuwe Frobeniusmatrix. We gaan dit na voor i = 1, j = 2 en k = 3.      1 0 0 1 0 0 1 0 0       0 0 1   a2 1 0  =  a3 0 1  0

1

0

a3

0

1 

a2 1

 =  a3 a2

1 0  0 0 1  1 00

0

0

0

 1

0

1

0

1

0



Voor twee Frobeniusmatrices Fi en Fj met i < j kunnen we het product gemakkelijk uitrekenen: 

Fi Fj =

1  ...   i 0  i+1 0 .  ..  0 j  j + 1 0  .. . 0

ai+1 .. .

j ··· 0 .. . ··· 0 ··· 0 .. .

··· aj · · · a j+1 .. .

··· 1 ··· 0 .. .

··· ···

···

··· 0

···

··· ··· ···

i 0 .. . 1

an

··· ··· ···

47

 0 1 ..   .. . .  0  0  0   0. ..    .   .. 0 0  0 0  . ..   .   .. 0 1

··· ··· ··· ··· ···

i 0 .. . 1 0 .. . 0 0 .. .

··· 0

···

j 0 .. . 0 0 .. . 1

··· · · · b j+1 .. .

 ··· 0 ..  . ··· 0  ··· 0 ..  .  ··· 0  ··· 0  ..  .

···

··· 1

··· ···

bn



=

1  ...   i 0  i+1 0 .  ..  0 j  j + 1 0  .. . 0

4.3

··· ··· ···

i 0 .. . 1 ai+1 .. .

··· ··· ···

j 0 .. . 0 0 .. .

··· aj · · · a j+1 .. .

··· 1 · · · b j+1 .. .

···

···

an

 ··· 0 ..  .  ··· 0  ··· 0 ..  .  ··· 0  ··· 0  ..  .

(4.9)

··· 1

bn

Gauss eliminatie en LU-decompositie

Beschouw het lineair stelsel AX = B. Door geschikte elementaire rijoperaties toe te passen op de matrix (A, B) verkrijgen we een rijequivalente matrix (A0 , B0 ) die in rijechelonvorm staat (zie de cursus “Algebra en Meetkunde”). Het stelsel A0 X = B0 is equivalent met het oorspronkelijk stelsel, en is gemakkelijk te bespreken en op te lossen. We bekijken vanaf nu het bijzonder geval waarin A een reguliere vierkante matrix is. We schrijven A0 = U en B0 = C. U is nu een bovendriehoeksmatrix, en het equivalente stelsel UX = C kan opgelost worden: men lost eerst de laatste vergelijking op, substitueert de oplossing in de voorlaatste, en gaat zo verder. Laten we kort herhalen hoe een matrix in rijechelonvorm kan gebracht worden. Stap 1 Men zorgt ervoor dat a11 6= 0, door eventueel de eerste rij te verwisselen met een andere. Om de numerieke stabiliteit te verhogen passen we deze eerste stap nog een beetje aan: we zorgen ervoor dat het element a11 in absolute waarde het grootste op de eerste kolom is. Kies dus de index r zo dat |ar1 | = max{|ai1 | : i = 1, . . . , n} en verwissel de eerste en de r-de rij. De nieuwe matrix noemen we A. Stap 2 Voor i = 2, . . . , n trekken we ai1 /a11 maal de eerste rij af van de i-de. We verkrijgen aldus een matrix van de vorm ! a A b 11 12 1 (A(1) , B(1) ) = e 0 A Be e B), e waarbij we de eerste rij en kolom van (A(1) , B(1) ) We herhalen nu het proc´ed´e voor de matrix (A, onberoerd laten. Als we dit n − 1 maal doen, dan krijgen we de gewenste matrix (A(n−1) , B(n−1) ) = (U,C). Merk op dat stap 1 er in feite op neer komt (A, B) links te vermenigvuldigen met de permutatiema-

48

trix P1r . Bij stap 2 vermenigvuldigen we links met de Frobeniusmatrix   1 0 ··· 0    −a21 /a11 1 · · · 0      −a /a F1 =  31 11 0 · · · 0   .. .. ..   . . .   −an1 /a11 0 · · · 1 We kunnen dus schrijven (A(1) , B(1) ) = F1 P1r (A, B) en de uiteindelijke matrix is dus van de vorm (U,C) = (A(n−1) , B(n−1) ) = Fn−1 P(n−1)rn−1 Fn−2 P(n−2)rn−2 · · · F1 P1r1 (A, B) waarbij ri ≥ i. In het bijzonder is Fn−1 P(n−1)rn−1 Fn−2 P(n−2)rn−2 · · · F1 P1r1 A = U een bovendriehoeksmatrix. Stel nu P = P(n−1)rn−1 P(n−2)rn−2 · · · P1r1 P is een permutatiematrix. Door herhaaldelijk gebruik te maken van (4.8) vinden we 0 Fn−1 Fn−2 · · · F20 F10 PA = U

waarbij de Fi0 allen Frobeniusmatrices zijn. Tenslotte is 0

0

−1 PA = F1−1 F2−1 · · · Fn−1 U = LU

waarbij

0

0

−1 L = F1−1 F2−1 · · · Fn−1

een benedendriehoeksmatrix is, die gemakkelijk uit te rekenen is door (4.9). We hebben met andere woorden het volgend resultaat bewezen. Stelling 4.3.1 (LU decompositie) Voor elke vierkante reguliere matrix A bestaat er een permutatiematrix P, een benedendriehoeksmatrix L en een bovendriehoeksmatrix U zodat PA = LU

(4.10)

Wat is het belang van zulk een decompositie? Als we van de matrix A de LU-decompositie kennen, dan kunnen we het stelsel AX = B (4.11) gemakkelijk oplossen, ongeacht het rechterlid B. Immers, (4.11) is equivalent met PAX = LUX = PB 49

Los eerst het stelsel LY = PB op. Dit kan gemakkelijk, omdat L een benedendriehoeksmatrix is. Los daarna het stelsel UX = Y op. Ook dit is gemakkelijk, want U is een bovendriehoeksmatrix. Hoe kunnen we L, U en P in de praktijk bepalen? Uiteindelijk hoeven we enkel de matrix A door elementaire rijoperaties in bovendriehoeksvorm te brengen. We vatten het algoritme samen als volgt. (s)

Algoritme 4.3.2 Stel T (0) = A. T (s) wordt uit T (s−1) berekend als volgt. We schrijven ti j voor het element op plaats (i, j) in de matrix T (s) . stap 1 Kies de index r zo dat (s−1)

|trs

(s−1)

| = max{|tis

| : i ≥ s}

(s−1)

Als trs = 0, dan is de matrix A singulier. stap 2 Verwissel de r-de en de s-de rij. We noemen de nieuwe matrix T . stap 3 Stel nu tis

(s)

= lis = t is /t ss voor i > s

(s) tik (s) tik

= t ik − list sk voor i, k > s = t ik anders

Stel T = T (n−1) . Dan is 1

0

··· 0

t11

t12

· · · t1n

  t21 L=  ... 

1 .. .

  0 ··· 0 en U =  ..    ... . 

t21 .. .

 · · · t2n  ..  .  

tn1

tn2

··· 1

0

· · · tnn







0



Opmerkingen 4.3.3 1) Uit (4.9) blijkt dat de elementen op de j-de kolom van L al bekend zijn bij de j-de iteratie (op eventuele verwisselingen van de rijen achteraf na). We plaatsen deze onder de diagonaal, op plaatsen die anders nul zijn. Dit spaart geheugenplaatsen uit. 2) Ook de j-de rij van U is bekend vanaf de j-de iteratie. 3) In ons algoritme passen we de zogenaamde parti¨ele pivotering toe. Bij volledige pivotering vervangt men de eerste en tweede stap in het algoritme door: Kies r en l zodat (s−1) (s−1) |trl | = max{|ti j | : i, j ≥ s} Verwissel dan de r-de en de s-de rij, en de l-de en s-de kolom. Men moet dan niet alleen links, maar ook rechts met een permutatiematrix vermenigvuldigen, zodat ook de derde stap ingewikkelder wordt. We gaan hier niet verder in op de details. Volledige 50

pivotering wordt toegepast als de matrix niet ge¨ekwilibreerd is, dit wil zeggen dat de matrix elementen van verschillende grootteorden bevat. 4) We kunnen ook het aantal bewerkingen berekenen dat nodig is om de LU-decompositie uit te voeren (en dus om het lineaire stelsel op te lossen. Het aantal bewerkingen bij iteratiestap s is (n − s) + 2(n − s)2 Gebruik makende van de formules n

∑i

i=1 n 2

∑i

i=1

=

n(n + 1) 2

=

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

vinden we voor het aantal bewerkingen n−1

n−1

s=1

i=1

∑ ((n − s) + 2(n − s)2) =

∑ (i + 2i2)

n(n − 1) (3 + 4n + 2) 6 n(n − 1)(4n + 5) = 6 =

en dit stijgt dus volgens 2n3 /3. Het aantal bewerkingen stijgt dus snel met het aantal vergelijkingen. Anders gezegd: als we het aantal vergelijkingen en onbekenden verdubbelen, dan wordt het aantal nodige bewerkingen met 8 vermenigvuldigd. 5) Onderstel dat er bij de Gauss eliminatie geen pivoteringen nodig zijn. Om de LU-decompositie te bepalen kan men dan ook als volgt te werk gaan. In L staan er (n − 1)n/2 onbekende elementen, namelijk degenen die onder de diagonaal staan. In U staan er n(n + 1)/2 onbekende elementen: de elementen op en boven de diagonaal. In het totaal zijn er dus n2 onbekenden te bepalen, en hiervoor heeft men n2 vergelijkingen, namelijk LU = A Dit stelsel met n2 vergelijkingen en n2 onbekenden blijkt gemakkelijk op te lossen te zijn. Alnaargelang de volgorde waarin men de vergelijkingen oplost spreekt men van de methodes van Crout of Banachiewiz. Deze blijken numeriek en theoretisch equivalent te zijn met de eliminatiemethode van Gauss. Als men hierin ook de pivotering wil inbouwen, dan worden de algoritmes ingewikkelder. Voorbeeld 4.3.4 Beschouw de matrix 



1

1

1

 A = T (0) =  2

1

 2

0

1

3

51

eerste iteratie: s = 1 stap 1: r = 2 stap 2: verwissel de eerste en de tweede rij. Dit komt eropneer A links te vermenigvuldigen met P12 .   2 1 2   T = 1 1 1 0

1

3

stap 3 

2

1

 T (1) =  1/2

2

1/2

0



 0

1

3

tweede iteratie: s = 2 stap 1: r = 3 stap 2: verwissel de tweede en de derde rij. Dit komt eropneer T (1) links te vermenigvuldigen met P23 .   2 1 2   T = 0 1 3 1/2

1/2

0

2

1

2

 T (2) =  0

1

 3 

stap 3 

1/2 We vinden dus



1/2



−3/2 





1

0

0

2

1

2

 L= 0

1

  0  en U =  0

1

 3 

1

0

1/2

1/2

0

−3/2

Verder is P = P23 P12 . A linksvermenigvuldigen met P komt erop neer de cyclische permutatie ( 1 3 2 ) op de rijen van A toe te passen, met andere woorden   2 1 2   PA =  0 1 3  1

1

1

Met kan onmiddellijk verifi¨eren dat LU = PA, zoals het moet.

52

4.4

Tridiagonale stelsels

Onderstel dat de matrix A van het stelsel AX volgende vorm is.  a1 c1 0 b a c 2 2  2   0 b3 a3  A= . .. ..  .. . .   0 0 0 In dit geval kunnen L en U  1 0 0 0 β  2 1 0 0   0 β3 1 0  L= . .. .. ..  .. . . .   0 0 0 0 0

0

0

= B tridiagonaal is, dit wil zeggen dat A van de 0

···

0

0

···

0

c3 .. .

···

0 .. .

0

· · · an−1

0



0    0   ..  .    cn−1 

0 0 0 0 · · · bn an in bidiagonale vorm gekozen worden, met andere woorden    ··· 0 0 α1 c1 0 0 · · · 0 0 0 α ··· 0 0 c2 0 · · · 0 0  2          ··· 0 0 0 0 α c · · · 0 0 3 3   en U =  .  . . . . . .. ..  . . . . . .  . . . . . .  . .       0 0 0 0 · · · αn−1 cn−1  ··· 1 0

0 · · · βn

0

1

0

0

0

···

0

αn

Om dit aan te tonen gebruiken we de methode van Crout: Als we de elementen op de eerste rij in de identiteit LU = A aan elkaar gelijkstellen dan vinden we a1 = α1 en c1 = c1 De tweede rij levert

b2 = β2 α1 , a2 = β2 c1 + α2 en c2 = c2

In het algemeen levert de i-de rij bi = βi αi−1 , ai = βi ci−1 + αi en ci = ci We krijgen de volgende recursieve formule, die toelaat om αi en βi te bepalen: α1 = a1 , en voor i = 2, . . . , n is bi αi−1 = ai − βi ci−1

βi =

(4.12)

αi

(4.13)

De methode werkt niet als een van de αi nul is. Dit doet zich voor als det(A) = 0, of als pivoteringen nodig zijn vooraleer men de matrix in LU-vorm kan ontbinden. Dit laatste geval doet zich bijvoorbeeld voor bij de matrix ! 0 1 A= 1 0 53

Opmerkingen 4.4.1 1) Men kan bovenstaande opmerking veralgemenen tot matrices met bandbreedte m. Dit zijn matrices waarvan alleen de elementen op de hoofddiagonaal en de m − 1 diagonalen boven en onder de hoofddiagonaal verschillend van nul zijn. Voor L kan men dan een matrix nemen met enkel elementen verschillend van nul op de hoofddiagonaal en de m − 1 diagonalen eronder. Analoog voor U. 2) De methode is ook van toepassing in de situatie waarbij A1

C1

0

0

···

0

B  2  0  A= .  ..   0

A2

C2

0

···

0

B3 .. .

A3 C3 .. .. . .

···

0 .. .

0

0

0

· · · An−1

0

0

0

0

···



Bn

0



0    0   ..  .    Cn−1  An

waarbij Ai , Bi en Ci vierkante matrices van een “kleine”dimensie zijn. In het algoritme (4.10)-(4.11) moet men dan vermenigvuldigen met de inverse van de matrix αi .

4.5

Singuliere waarden ontbinding

Neem twee Euclidische ruimten E en F, van dimensies respectievelijk N en m, en een lineaire afbeelding f : E → F. Het is onze bedoeling om orthonormale basissen voor E en F te vinden tenopzichte waarvan de matrix van f zo eenvoudig mogelijk wordt. We herhalen eerst dat de toegevoegde afbeelding f† : F → E gekarakteriseerd wordt door de formule (zie “Algebra en Meetkunde”, §6.3). h~x, f † (~y)i = h f (~x),~yi voor elke ~x ∈ E en ~y ∈ F. We bekijken nu de samenstelling f†◦ f : E → E De afbeelding f † ◦ f is zelftoegevoegd, en bijgevolg bestaat er een orthonormale basis (zie “Algebra en Meetkunde”, §7.3) U = {~u1 , · · · ,~un } van E zodat de matrix [ f † ◦ f ]U ,U een diagonaalmatrix is. Bovendien zijn alle eigenwaarden nietnegatief: als u een (genormaliseerde) eigenvector is met eigenwaarde λ, dan is ( f † ◦ f )(~u) = λ~u

54

en dus λ = = = =

h~u, λ~ui h~u, ( f † ◦ f )(~u)i h f (~u), f (~u)i k f (~u)k2 ≥ 0

We kunnen de volgorde in de basis U dus zo herschrijven dat  2 σi ~ui als i ≤ r † ( f ◦ f )(~ui ) = ~0 als i > r Hierbij is r = rg( f † ◦ f ). Voor i = 1, · · · , r stellen we ~vi =

f (~ui ) σi

Dan is f (~ui ) f (~u j ) , i σi σj 1 = h~ui , ( f † ◦ f )(~u j )i σi σ j 1 = h~ui , σ j~u j i σi σ j

h~vi ,~v j i = h

σ2j

=

σi σ j = δi j

δi j

Vul {~v1 , · · · ,~vr } aan tot een orthonormale basis V = {~v1 , · · · ,~vm } van F. Schrijf   σ1 0 · · · 0    0 σ2 · · · 0   D1 =  . .. ..  . . .   .  0 Voor i = 1, · · · , r geldt

0

· · · σr

f (~ui ) = σi~vi

en voor i = r + 1, · · · , m hebben we ( f † ◦ f )(~ui ) = ~0, zodat h f (~ui ), f (~ui )i = h~ui , ( f † ◦ f )(~ui )i = 0 zodat f (~ui ) = ~0 55

We vinden daarom [ f ]U ,V =

D1

0

0

0

!

Hieruit volgt ook dat rg( f ) = r = rg( f † ◦ f ) en we hebben bewezen: Stelling 4.5.1 Zij f : E → F een lineaire afbeelding tussen twee eindigdimensionale Euclidische ruimten. Dan is r = rg( f ) = rg( f † ◦ f ), en er bestaan orthonormale basissen U en V van respectievelijk E en F zodat ! D1 0 [ f ]U ,V = 0 0 waarbij D1 een r × r-diagonaalmatrix is, met positieve elementen σi op de diagonaal. We kunnen dit resultaat herinterpreteren in termen van matrices. Zij A een re¨ele m × n-matrix, en bekijk de lineaire afbeelding f : Rn → Rm . We rusten Rn en Rm uit met het standaard inwendig product. Uit stelling 4.5.1 volgt het bestaan van orthonormale basissen U en V van respectievelijk Rn en Rm zodat ! D1 0 [ f ]U ,V = 0 0 Zij M en N de overgangsmatrices van de standaardbasissen voor Rn en Rm naar U en V . Dan zijn M en N orthogonale matrices, en modat A de matrix van f is tenopzichte van de standaardbasissen volgt ! D 0 1 N t AM = 0 0 Stelling 4.5.2 Onderstel dat A een m × n-matrix is. Er bestaat een m × m-orthogonale matrix N en een n × n-orthogonale matrix M zodat ! D1 0 t ∈ Mmn (R) (4.14) N AM = 0 0 Men noemt de getallen σ1 , σ2 , · · · , σr de singuliere waarden van de matrix A. (4.14) noemt men de singuliere waardenontbinding van A. Meetkundige interpretatie We werken met de hierboven beschreven basissen U van E en V van F. Beschouw de eenheidsbol B in E. Deze heeft als vergelijking x12 + x22 + · · · + xn2 ≤ 1 56

Dit betekent dat we B onder de volgende vorm kunnen schrijven. ( ) n

n

i=1

i=1

∑ xi~ui | ∑ xi2 ≤ 1

B=

Het beeld van B onder de lineaire afbeelding f is nu ( n

n

i=1

i=1

)

∑ xi f (~ui) | ∑ xi2 ≤ 1

f (B) = (

r

n

i=1

i=1

= (

∑ xiσi~vi | ∑ xi2 ≤ 1 r

r

y2 ∑ yi~vi | ∑ σi2 ≤ 1 i=1 i=1 i

=

)

)

Het beeld van de eenheidsbol is dus de r − 1-dimensionale ellipsoide met vergelijking r

y2i ∑ σ2 ≤ 1 i=1 i yr+1 = · · · = ym = 0 De singuliere waarden zijn dus de assen van het beeld van de eenheidsbol. Het conditiegetal van een matrix Onderstel dat A een vierkante reguliere matrix is (dus rg(A) = r = n = m). Het beeld van de eenheidsbol is dan de ellipsoide met kleinste as σmin en grootste as σmax . Voor elke X ∈ Rn geldt dan σmin kXk ≤ kAXk ≤ σmax kXk Beschouw nu het lineair stelsel AX = B Een fout ∆B op het rechterlid veroorzaakt een fout ∆X op de oplossing X: A(X + ∆X) = B + ∆B of

A∆X = ∆B

Omdat σmin k∆Xk ≤ k∆Bk en kBk ≤ σmax kXk vinden we voor de relatieve fout op X de afschatting k∆Xk σmax k∆Bk ≤ kXk σmin kBk Per definitie noemen we k=

σmax σmin

57

het conditiegetal van de matrix A. k is een betere maat voor de conditie van de matrix A dan de determinant van A. Beschouw bijvoorbeeld de 100 × 100-matrix met 0.1 op de hoofddiagonaal en nullen elders. De determinant is heel klein: det(A) = 10−100 maar het conditiegetal is 1. A zet de eenheidsbol om in een kleinere bol.

4.6

Veralgemeende inversen van matrices

Onderstel dat een probleem gesteld is, waarbij gevraagd wordt om vijf onbekende parameters te bepalen. Informatie over de parameters kan men verkrijgen door metingen uit te voeren. Men voert vijftig metingen uit. Men beschikt dan eigenlijk over vijftig betrekkingen tussen de vijf onbekende parameters. Dit zijn natuurlijk meer vergelijkingen dan nodig. Door de meetfouten zal het stelsel van vijftig vergelijkingen met vijf onbekenden een strijdig stelsel zijn. De zaak is nu om de meest plausibele waarde voor de vijf onbekende parameters naar voor te schuiven. We zullen dit hier doen in het geval dat de gegeven vergelijkingen lineair zijn. Desnoods linearizeren we ons stelsel vergelijkingen: we benaderen onze vergelijkingen door lineaire. Onze taktiek is de volgende: beschouw een lineair stelsel AX = B van m vergelijkingen met n onbekenden. Het stelsel heeft niet noodzakelijk een oplossing, en daarom stellen we ons tevreden met de vectoren X waarvoor kAX − Bk minimaal is. We werken hierbij met de gewone Euclidische norm. Als A de matrix van een niet-injectieve lineaire afbeelding is, dan zal kAX − Bk minimaal zijn voor meer dan een waarde van X. We kiezen de oplossing X van minimale lengte. In het vervolg zullen we zien dat de afbeelding die B afbeeldt op X zoals hierboven beschreven, lineair is. We noemen deze de veralgemeende inverse of Moore-Penrose inverse van A. Zoals in de voorgaande paragraaf beschouwen we twee Euclidische ruimten E en F, van respectievelijk dimensie n en m, en een lineaire afbeelding f : E−→F We noteren verder K = Ker ( f ) ⊂ E en L = Im ( f ) ⊂ F We kunnen nu schrijven E = K ⊕ K ⊥ en F = L ⊕ L⊥ We gebruiken volgende notaties voor de orthogonale projecties pK : E → K

pK ⊥ : E → K ⊥

pL : F → L

pL⊥ : E → L⊥ 58

De lineaire afbeelding f beperkt zich tot een bijectie h = f|K ⊥ : K ⊥ −→L Per definitie stellen we f + = h−1 ◦ pL

(4.15)

Neem ~y ∈ F, en onderstel dat f + (~y) =~x. Dan is f (~x) = pL (~y) en ~x ∈ K ⊥ Herhaal dat pL (~y) de vector in L is die het dichtst bij ~y gelegen is. Dus met andere woorden k~z −~yk ≥ kpL (~y) −~yk voor elke~z ∈ L. Voor elke ~x1 ∈ E is f (~x1 ) ∈ L, en dus is k f (~x1 ) −~yk ≥ kpL (~y) −~yk = k f (~x) −~yk en dus is k f (~x) −~yk minimaal voor ~x = f + (~y). Als nu~x1 een andere vector is waarvoor k f (~x1 )−~yk minimaal is, dan is f (~x) = f (~x1 ), en~x−~x1 ∈ K. Maar dan is ~x1 =~x + (~x1 −~x) ∈ K ⊥ ⊕ K en k~x1 k2 = k~xk2 + k~x1 −~xk2 ≥ k~xk2 zodat we bewezen hebben dat ~x van minimale lengte is. Dit is de procedure die we hierboven beschreven. We zullen nu de eigenschappen van f + bestuderen. Merk om te beginnen op dat f + = f −1 als f een isomorfisme is.

Stelling 4.6.1 f + ◦ f = pK ⊥ en f ◦ f + = pL

(4.16)

f + ◦ f en f ◦ f + zijn zelftoegevoegd. Dit wil zeggen dat f † ◦ f +† = f + ◦ f f +† ◦ f † = f ◦ f +

(4.17) (4.18)

f ◦ f+◦ f = f f+◦ f ◦ f+ = f+

(4.19) (4.20)

Bovendien geldt

59

Bewijs. f ◦ f + = pL volgt onmiddellijk uit de definitie van f + . Neem ~x ∈ E willekeurig, en stel f + ( f (~x)) = ~x1 . Dan is ~x1 ∈ K ⊥ , en f (~x1 ) = f (~x). Dus is f (~x1 − ~x) = ~0, en ~x1 −~x ∈ K. Dit wil juist zeggen dat ~x1 de orthogonale projectie van ~x op K ⊥ is. (4.17) en (4.18) volgen nu onmiddellijk uit het feit dat orthogonale projecties zelftoegevoegd zijn. Tenslotte hebben we, voor alle ~x ∈ E en ~y ∈ F ( f ◦ f + ◦ f )(~x) = pL ( f (~x)) = f (~x) ( f + ◦ f ◦ f + )(~y) = f + (pK ⊥ (~y)) = (h−1 ◦ pK ⊥ ◦ pK ⊥ )(~y) = f + (~y)  We zullen nu zien dat (4.17)-(4.20) A+ volledig karakteriseren. Stelling 4.6.2 Onderstel dat g : F → E voldoet aan volgende eigenschappen: g◦ f f ◦g f ◦g◦ f g◦ f ◦g

= = = =

f † ◦ g† g† ◦ f † f g

(4.21) (4.22) (4.23) (4.24)

Dan is g = f + . Bewijs. g = (4.19) = (4.17), (4.18), (4.21), (4.22) = = (4.23) = (4.17), (4.18) = (4.20) =

g◦ f ◦g g◦ f ◦ f+◦ f ◦ f+◦ f ◦ f+◦ f ◦g f † ◦ g† ◦ f † ◦ f +† ◦ f + ◦ f +† ◦ f † ◦ g† ◦ f † ( f ◦ g ◦ f )† ◦ f +† ◦ f + ◦ f +† ◦ ( f ◦ g ◦ f )† f † ◦ f +† ◦ f + ◦ f +† ◦ f † f+◦ f ◦ f+◦ f ◦ f+ f+ 

Gevolg 4.6.3 Voor een lineaire afbeelding f : E → F geldt dat †

f ++ = f en f + = ( f † )+ Bewijs. Vervang in stelling 4.6.2 f door f + en g door f . Aan (4.21)-(4.24) is voldaan, zodat f = f ++ . De tweede gelijkheid wordt op dezelfde manier bewezen.  Kies nu orthonormale basissen E en F voor E en F. Als [ f ]E ,F = A, dan is de veralgemeende inverse van de matrix A per definitie A+ = [ f + ]F ,E 60

Uit de twee bovenstaande stellingen volgt onmiddellijk dat A+ gekenmerkt wordt door volgende vergelijkingen. A+ A AA+ AA+ A A+ AA+

= = = =

At A+t A+t At A A+

(4.25) (4.26) (4.27) (4.28)

Hoe kunnen we de veralgemeende inverse van een matrix A berekenen? Als we de singuliere waarden ontbinding van A kennen, dan kennen we ook de veralgemeende inverse. We gebruiken dezelfde notatie als in stelling 4.5.1. We hebben ! D1 0 [ f ]U ,V = ∈ Mmn (R) 0 0 Van deze matrix kunnen we de veralgemeende inverse gemakkelijk uitrekenen: ! −1 D 0 1 [ f + ]V ,U = ∈ Mnm (R) 0 0 en dus is A+ = [ f + ]F ,E = N[ f + ]V ,U Mt = N

61

D−1 1

0

0

0

! Mt

(4.29)

Hoofdstuk 5 Eigenwaarden en eigenvectoren 5.1

Algemene eigenschappen

We herhalen definitie en eigenschappen van eigenwaarden en eigenvectoren, uit de cursus “Algebra en Meetkunde”. Zij A een (re¨ele) n × n-matrix. X is een eigenvector van A als X 6= 0, en AX = λX

(5.1)

We noemen λ een eigenwaarde van A. De eigenwaarden zijn de wortels van de karakteristieke vergelijking a11 − λ a12 ··· a1n a21 a22 − λ · · · a2n PA (λ) = det(A − λIn ) = . (5.2) .. .. = 0 . . .. a an2 · · · ann − λ n1 PA is een veelterm van graad n in λ, en wordt de karakteristieke veelterm genoemd. Men kan ook eigenwaarden en eigenvectoren van een willekeurige lineaire afbeelding defini¨eren (zie de cursus “Algebra en Meetkunde”). Als A en B matrices zijn van eenzelfde lineaire afbeelding, maar tenopzichte van verschillende basissen, dan is PA = PB , en bijgevolg hebben A en B dezelfde eigenwaarden. Als A n lineair onafhankelijke eigenvectoren heeft, dan vormen deze een basis van eigenvectoren. De lineaire afbeelding bepaald door A heeft een diagonale matrix tenopzichte van deze basis van eigenvectoren, en we zeggen dat A diagonalizeerbaar is. Herhaal verder de volgende eigenschappen. Stelling 5.1.1 Voor elke re¨ele symmetrische matrix A bestaat er een orthonormale basis van eigenvectoren. Bijgevolgd bestaat er een orthogonale matrix M zodat Mt AM = D diagonaal is. 62

Stelling 5.1.2 (formule van Cayley-Hamilton) Een vierkante matrix A voldoet steeds aan zijn karakteristieke vergelijking: PA (A) = 0 Om de eigenwaarden en eigenvectoren van een matrix A te bepalen kunnen we als volgt tewerk gaan. We schrijven de karakteristieke vergelijking (5.2) op, en we bepalen de wortels λi . Bij elke λi bepalen we de bijhorende eigenvectoren, door (5.1) op te lossen. Voor grote waarden van n is dit niet erg praktisch, en daarom zoeken we naar andere methoden. De cirkels van Gerschgorin Onderstel dat A een matrix is waarvan de waarden op de diagonaal veel groter zijn dan die buiten de diagonaal. Men kan er dan van uitgaan dat de waarden op de diagonaal dicht tegen de eigenwaarden liggen. Dit wordt bewezen in de volgende eigenschap. Stelling 5.1.3 Schrijf

n

ri =

∑

|ai j |

j=1, j6=i

Alle eigenwaarden van A liggen dan binnen de cirkels met middelpunt aii en straal ri . Bewijs. Onderstel dat AX = λX. Dan is n

∑ ai j x j = λxi

j=1

voor i = 1, 2, · · · , n. Kies nu de index i zodanig dat |xi | ≥ |x j | voor elke andere index j. Dan is n

(λ − aii )xi =

∑

ai j x j

j=1, j6=i

en

n

n

xj |λ − aii | ≤ ∑ |ai j | ≤ xi j=1, j6=i

∑

|ai j | = ri

j=1, j6=i

 De methode van Krylov Dit is een manier om de karakteristieke vergelijking op te stellen zonder een determinant te berekenen. De karakteristieke vergelijking is van de vorm n−1

λn + ∑ bi λi = 0 i=0

63

(5.3)

Uit de formule van Cayley-Hamilton weten we dat n−1

An + ∑ bi Ai = 0 i=0

We vermenigvuldigen dit met een willekeurige kolomvector Y . n−1

AnY + ∑ bi AiY = 0 i=0

Deze vectori¨ele betrekking kunnen we beschouwen als een lineair stelsel in b0 , b1 , · · · , bn−1 . Oplossen van dit stelsel (bijvoorbeeld met behulp van Gauss eliminatie) levert de co¨effici¨enten b0 , b1 , · · · , bn−1 en dus de karakteristieke vergelijking.

5.2

De methode van de machten

We onderstellen dat A diagonalizeerbaar is, dit wil zeggen dat er een basis van eigenvectoren bestaat. We schrijven {X1 , X2 , · · · , Xn } voor deze basis van eigenvectoren, en we noteren λ1 , λ2 , · · · , λn voor de bijhorende eigenwaarden. Onderstel nu dat de eigenwaarde λ1 in absolute waarde groter is dan alle andere, m.a.w. |λ1 | > |λi | voor i = 2, 3, · · · , n. Neem nu een willekeurige vector Y . We kunnen schrijven: n

Y = ∑ αi Xi i=1

Als we beide leden linksvermenigvuldigen met Am , dan vinden we n

n

i=1

i=1

AmY = ∑ αi Am Xi = ∑ αi λm i Xi Deel nu beide leden door λm 1.

en

n λm AmY = α X + α i 1 1 ∑ λim Xi λm 1 1 i=2

AmY = α1 X1 m→∞ λm 1 lim

(5.4)

We kunnen besluiten dat de richting van AmY in de limiet voor m → ∞ de richting aanneemt van de eigenvector X1 behorend bij de dominante eigenwaarde. Om λ1 te bepalen redeneren we als volgt: vermenigvuldig (5.4) met een willekeurige rijvector E t (bijvoorbeeld een basisvector): E t AmY lim = α1 E t X1 m→∞ λm 1 64

We hebben ook

E t Am+1Y = α1 E t X1 m→∞ λm+1 1 lim

Als we beide betrekkingen door mekaar delen, dan vinden we E t Am+1Y m→∞ E t AmY

λ1 = lim

(5.5)

Om X1 te bepalen nemen we dus een startvector Y , en berekenen achtereenvolgens AY , A2Y , A3Y , enz. Omdat enkel de richting van X1 belangrijk is, normaliseren we onze benadering voor de eigenvector bij elke iteratie: we zorgen ervoor dat de grootste component in absolute waarde 1 wordt. Opeenvolgende benaderingen van λ1 vinden we met behulp van (5.5). (0)

Algoritme 5.2.1 Kies een startvector Y = X1 die als grootste component in absolute waarde 1 heeft. (s) Onderstel dat een benadering X1 van X1 gekend is, en dat de grootste component in absolute (s) waarde van X1 gelijk is aan 1. (s) (s) (s) Bepaal AX1 . λ1 is de component van AX1 die het grootste is in absolute waarde, en (s)

(s+1) X1

=

X1

(s)

λ1

In bovenstaand algoritme benaderen we λ1 met behulp van (5.5): we nemen voor Y een van de basisvectoren. Welk is de beste keuze voor Y . Of met andere woorden, als we een benadering X hebben voor de eigenvector X1 , welk is dan de beste schatting voor de bijhorende eigenwaarde. Dit is het getal µ waarvoor kAX − µXk minimaal is. Lemma 5.2.2 kAX − µXk is minimaal voor µ=

X t AX kXk2

(5.6)

Bewijs. kAX − µXk2 = (AX − µX)t (AX − µX) = kXk2 µ2 − 2X t AXµ + kAXk2 We stellen de afgeleide naar µ gelijk aan nul: 2kXk2 µ − 2X t AX = 0 Het maximum wordt dus bereikt voor µ=

X t AX kXk2 

65

Opmerking 5.2.3 Als λ1 een meervoudige eigenwaarde is, dan werkt de methode nog altijd. Onderstel dat met λ1 twee lineair onafhankelijke vectoren X1 en X2 overeenstemmen. Als we redeneren zoals voor (5.4), dan vinden we dat AmY = α1 X1 + α2 X2 m→∞ λm 1 lim

Als we vertrekken met een andere startvector Y 0 die geen veelvoud is van Y , dan vinden we AmY 0 = α01 X1 + α02 X2 m→∞ λm 1 lim

en op deze manier vinden we twee lineair onafhankelijke eigenvectoren die behoren bij λ1 . Opmerking 5.2.4 Wat gebeurt er als toevallig α1 = 0? In dit geval is AmY =0 m→∞ λm 1 lim

In de praktijk zal dit zich niet voordoen: door afrondingsfouten krijgt AmY een component langs X1 die niet helemaal nul is. Initieel zal AmY de richting van de eigenvector X2 behorend bij de in absolute waarde tweedegrootste eigenwaarde aannemen. Nadien convergeert de rij toch naar de richting van X1 . Opmerking 5.2.5 Als |λ2 /λ1 | ≈ 1, dan zal de methode der machten traag convergeren. Men kan de convergentie versnellen door een verschuiving uit te voeren: men vervangt A door A − sI. De eigenvectoren blijven dan dezelfden, en de bijhorende eigenwaarden zijn λ0i = λi − s. Door een geschikte keuze van s wordt de verhouding |λ2 /λ1 | kleiner, zodat de convergentie versnelt. Als bijvoorbeeld λ1 = 8, λ2 = −7.9 en λ3 = 2, dan is |λ2 /λ1 | ≈ 1, zodat de convergentie traag verloopt. Kiest men s = −4, dan worden de nieuwe eigenwaarden λ01 = 12, λ02 = −3.9 en λ03 = 6. Nu is |λ02 /λ01 | = 1/2 zodat de convergentie verbetert. Het nadeel aan deze methode is dat we s op goed geluk af moeten kiezen. Opmerking 5.2.6 Met de techniek uit de vorige opmerking kunnen we nog een tweede eigenwaarde en bijhorende eigenvector bepalen. Neem weer het voorbeeld uit voorgaande opmerking: λ1 = 8, λ2 = −7.9 en λ3 = 2, en stel nu s = 4. De nieuwe eigenwaarden worden nu λ01 = 4, λ02 = −11.9 en λ03 = −2. λ02 is dus dominant geworden. Bepalen van de overige eigenwaarden Zoals blijkt uit de vorige opmerking kunnen we met behulp van de methode der machten steeds de grootste en de kleinste eigenwaarde, en de bijhorende eigenwaarden bepalen. Onderstel nu dat A re¨eel en symmetrisch is. We kunnen dan nog verder gaan. Onderstel dat de dominante eigenwaarde λ1 bepaald is, en normaliseer de bijhorende eigenvector X1 in Euclidische norm: onderstel dat kX1 k = 1. Stel nu B = A − λ1 X1 X1t 66

Merk op dat X1 X1t het product is van een 3 × 1 en een 1 × 3 matrix is, en dus een 3 × 3-matrix. A is re¨eel en symmetrisch, en dit betekent dat de (genormaliseerde) eigenvectoren een orthonormale basis van Rn vormen (zie stelling 5.1.1). Voor i 6= 1 vinden we daarom dat BXi = AXi − λ1 X1 X1t Xi = λXi en, voor i = 1:

BX1 = AX1 − λ1 X1 X1t X1 = 0

De eigenvectoren van B zijn dus X1 , X2 , . . . , Xn , met bijhorende eigenwaarden 0, λ2 , . . . , λn . Toepassen van de methode der machten levert een eigenwaarde verschillend van λ1 . In theorie kan men op deze manier alle eigenwaarden en eigenvectoren bepalen. Voor grotere waarden van n stapelen de fouten zich op, en dit brengt mee dat men met de methode slechts enkele eigenwaarden kan bepalen: degenen die het grootst zijn in absolute waarde. De inverse machtenmethode Onderstel dat λ 6= 0 een eigenwaarde is van A, met bijhorende eigenvector X: AX = λX Als A regulier is, dan vinden we

A−1 X = λ−1 X

en dus is λ−1 een eigenwaarde van A−1 , en de bijhorende eigenvector is X. Onderstel nu dat µ geen eigenwaarde is van A. Dan is 0 geen eigenwaarde van A − µIn , en A − µIn is dus regulier. De eigenwaarden van (A − µIn )−1 zijn dus (λi − µ)−1 (i = 1, . . . , n). Als µ een benadering is voor e´ e´ n van de eigenwaarden, bijvoorbeeld λk , dan is (λk − µ)−1 de dominante eigenwaarde. Als µ zeer dicht bij λk , dan is deze eigenwaarde v´ee´ l groter dan alle anderen, zodat de methode der machten (toegepast op (A − µIn )−1 ) zeer snel convergeert. Als we de eigenwaarden ongeveer weten liggen (bijvoorbeeld met behulp van de cirkels van Gershgorin), dan kunnen we deze allen bepalen met deze methode. Nadeel is dat we bij elke iteratieslag moeten vermenigvuldigen met de inverse van A − µIn . We moeten dus een lineair stelsel oplossen. Een verfijning van deze methode is de volgende. Rayleigh quotient iteratie Men vertrekt van een benadering µ0 van een van de eigenwaarden, en van een benadering X0 van de bijhorende eigenvector. Men voert e´ e´ n iteratie van de inverse machtenmethode uit. Men vindt een benadering X1 voor de eigenvector, en bepaalt hieruit een nieuwe benadering µ1 van de eigenwaarde, met behulp van lemma 5.2.2: men stelt µ1 gelijk aan het Rayleigh quotient µ1 =

X1t AX1 kX1 k2

Men herhaalt dan het proc´ed´e. We verkrijgen volgend algoritme. 67

Algoritme 5.2.7 Vertrek van benaderingen µ0 en X0 van de eigenwaarde λ en de bijhorende eigenvector X. µk+1 en Xk+1 worden uit µk en Xk berekend door oplossing van het stelsel (A − µk In )Xk+1 = Xk Men vervangt Xk+1 door Xk+1 /kXk=1 k, en men stelt t µk+1 = Xk+1 AXk+1

Dit algoritme convergeert zeer snel: als A re¨eel en symmetrisch is, dan kan men bewijzen dat het van orde drie is.

5.3

De methode van Jacobi

Onderstel dat A re¨eel en symmetrisch is. Dan weten we dat er een orthogonale matrix Q bestaat zodanig dat Qt AQ = D een diagonaalmatrix is. De elementen op de diagonaal zijn de eigenwaarden van A, en de kolommen van Q zijn de (genormaliseerde) eigenvectoren van A. Als we Q kunnen bepalen is het probleem opgelost. Een typevoorbeeld van een orthogonale transformatie is een rotatie. Onze taktiek bestaat er nu in om opeenvolgend rotaties uit te voeren zodanig dat de getransformeerde matrix A dichter en dichter bij een diagonaalmatrix komt. Onderstel dat we een rotatie over een hoek α uitvoeren in het vlak bepaald door de r-de en de s-de eenheidsvector. De orthogonale matrix N die deze rotatie bepaalt is r s   1 0 ··· 0 ··· 0 ··· 0  0 1 ··· 0 ··· 0 ··· 0 . . . . ..   .. .. .. .. .    r  0 0 · · · cos α · · · − sin α · · · 0  N =  .. .. .. .. ..  . . . . .   s  0 0 · · · sin α · · · cos α · · · 0  . . .. .. ..   .. .. . . . 0 0 ··· 0 ··· 0 ··· 1 Als we deze rotatie uitvoeren, dan wordt de nieuw matrix A = N t AN We gaan deze nu berekenen. Om de berekeningen te vereenvoudigen stellen we r = 1 en s = 2. In het algemeen zijn de berekeningen volledig analoog, maar veel omslachtiger om neer te schrijven. 68

Met notaties als in § 4.2 hebben we    a11 a12 A13 cos α    A =  a21 a22 A23  en N =  sin α A31

A32

A33

− sin α

0

cos α

0

0



 0  In−2

zodat t A = N AN cos α  =  − sin α



0 cos α

 =  − sin α 0

sin α

0

a11

a12

A13

cos α

 0   a21

a22

 A23   sin α

0 sin α cos α 0





cos α

In−2 A31 A32 A33  0 a11 cos α + a12 sin α  0   a12 cos α + a22 sin α

A31 cos α + A32 sin α

In−2

en we krijgen dus voor A volgende matrix :  a11 cos2 α + 2a12 cos α sin α + a22 sin2 α   cos α sin α(a22 − a11 ) + (cos2 α − sin2 α)a12

− sin α cos α

0



 0 

0 0 In−2 −a11 sin α + a12 cos α

A13

−a12 sin α + a22 cos α

 A23 



−A31 sin α + A32 cos α A33

cos α sin α(a22 − a11 ) + (cos2 α − sin2 α)a12 a22 cos2 α − 2a12 cos α sin α + a11 sin2 α

A31 cos α + A32 sin α

−A31 sin α + A32 cos α

cos αA13 + sin αA23



 − sin αA13 + cos αA23  A33 De matrix A heeft volgende eigenschappen: 1) A is re¨eel en symmetrisch. Immers t

A = (N t AN)t = N t AN = A 2) A en A hebben dezelfde eigenwaarden. Immers, A en A stellen dezelfde lineaire afbeelding voor, maar ten opzichte van verschillende (orthonormale) basissen. 3) ai j = ai j als i en j beiden verschillend zijn van r en s. Dit kunnen we zien aan de expliciete vorm van de matrix A. 4) n

∑

n

n

∑ a2i j = ∑

n

∑ a2i j

(5.7)

i=1 j=1

i=1 j=1

Om (5.7) te bewijzen gaan we als volgt te werk. Onderstel dat N een orthogonale matrix is, en X een kolomvector. Stel Y = NX. Dan is n

n

i=1

i=1

∑ y2i = kY k2 = kNXk2 = kXk2 = ∑ xi2 69

Immers, N is orthogonaal, en behoudt dus de lengte. Hieruit volgt dat, voor een vierkante matrix A, en B = NA: n

n

∑∑

b2i j

i=1 j=1

n

=∑

n

∑ a2i j

(5.8)

i=1 j=1

Als nu X een rijvector is, en N een orthogonale matrix, dan hebben we, met Y = XN n

n

i=1

i=1

∑ y2i = kY t k2 = kNt X t k2 = kX t k2 = ∑ xi2

en we kunnen hieruit concluderen dat (5.8) ook nog geldt als B = AN. Combineren van de twee eigenschappen geeft (5.7). Het basisidee van de methode van Jacobi bestaat er nu in om α zo te kiezen dat ars = asr = 0. Uit de expliciete formule voor A zien we dat ars = asr = sin α cos α(ass − arr ) + ars (cos2 α − sin2 α)

(5.9)

en het volstaat dus om α zo te kiezen dat tg 2α =

2ars arr − ass

(5.10)

Uit (5.7), toegepast op het geval n = 2, volgt a2rr + a2rs + a2sr + a2ss = a2rr + a2rs + a2sr + a2ss Als ars = asr = 0, dan hebben we dus dat a2rr + a2ss = a2rr + 2a2sr + a2ss en

n

n

i=1

i=1

∑ a2ii = ∑ a2ii + 2a2sr

De “massa¨op de hoofddiagonaal is dus toegenomen. Stel A1 = A. We passen dezelfde procedure nogmaals toe, en maken een ander niet-diagonaal element nul. We verkrijgen een nieuwe matrix A2 . Jammer genoeg wordt dan het element ars dat we in de eerste stap nul maakten, opnieuw verschillend van nul. Belangrijk evenwel is dat de “massa¨op de hoofddiagonaal opnieuw toeneemt. Men verkrijgt aldus een rij matrices A0 , A1 , · · ·, met stijgende massa op de hoofddiagonaal, en, vanwege (5.7), dalende massa buiten de hoofddiagonaal. Men hoopt dat lim Ak = D

k→∞

een diagonaalmatrix is. In dat geval is ∞

Q = lim N0 N1 · · · Nk = ∏ Nk k→∞

k=0

de gevraagde orthogonale matrix. De kolommen van Q zijn de genormalizeerde eigenvectoren. 70

Praktische werkwijze Er zijn drie methodes. 1) De klassieke werkwijze Dit is de meest voor de hand liggende methode. Het is de bedoeling om de massa buiten de hoofddiagonaal zo snel mogelijk nul te maken. Bij elke iteratie neemt de massa op de hoofddiagonaal toe met tweemaal het kwadraat van de nul gemaakte term. Men zoekt dus telkens de grootste |ars | (met r 6= s), en maakt dit element nul. Dit optimaliseert het aantal iteraties. Het nadeel is dat men telkens het grootste element moet zoeken. 2) De reeksmethode Men doorloopt de niet-diagonale elementen in een vast bepaalde volgorde, en maakt deze telkens nul. Het aantal nodige iteraties is dan uiteraard veel groter dan bij de klassieke werkwijze. 3) De drempelmethode Men doorloopt de niet-diagonale elementen in een vast bepaalde volgorde, maar men voert de iteratie enkel uit als het element |ars | in kwestie boven een gegeven drempelwaarde ligt. Als alle elementen buiten de diagonaal kleiner zijn geworden dan de drempel, dan herbegint men met een lagere drempel. Men kan aantonen dat de methode van Jacobi convergeert in de drie gevallen. Berekening van cos α en sin α (5.10) levert tg 2α. Om de matrix N neer te schrijven hebben we echter cos α en sin α nodig. Uit (5.10) volgt dat 2ars R arr − ass cos 2α = R sin 2α =

met R=

q

(arr − ass )2 + 4a2rs

Hieruit halen we dat 1 arr − ass + (5.11) 2 2R 1 arr − ass sin2 α = − (5.12) 2 2R ars sin α cos α = (5.13) R Bij het gebruik van deze formules kan een numerieke instabiliteit optreden. Indien het algoritme convergeert worden de elementen ars buiten de diagonaal na verloop van tijd erg klein, zodat cos2 α =

R2 ≈ (arr − ass )2 en

arr − ass 1 ≈± 2R 2 71

(5.14)

Als het rechterlid van (5.14) −1/2 is, dan is (5.11) numeriek instabiel. Als het rechterlid 1/2 is, dan is (5.12) numeriek instabiel. De numerieke instabiliteit wordt omzeild door het algoritme als volgt aan te passen. Als arr > ass , dan gebruikt men (5.11) en (5.13). Als arr < ass , dan gebruikt men (5.11) en (5.12).

72

Hoofdstuk 6 Interpolatie De oplossingen van de problemen in de voorgaande hoofdstukken waren steeds getallen (oplossing van vergelijkingen) of vectoren (oplossing van stelsels). Vanaf nu zullen we te maken hebben met problemen waarbij zowel de gegevens als de oplossingen mogelijk functies zijn. We zullen dus te doen hebben met functies die gegeven zijn in een numerieke vorm, dit wil zeggen dat we de functies slechts kennen in een aantal spilpunten x0 , x1 , . . . , xn . Anders gezegd, de functie y = f (x) is gegeven door een tabel x0 x1 .. .

f0 = f (x0 ) f1 = f (x1 )

xn

fn = f (xn )

Onderstel dat f in deze vorm gegeven is. Een vraagstuk dat zich nu stelt is het volgende: gegeven een waarde x die niet gelijk is aan een van de spilpunten. Wat is een plausibele waarde voor f (x)? In dit hoofdstuk behandelen we het vraagstuk van de interpolatie: we bepalen een veelterm pn van graad n zodanig dat pn (xi ) = fi (6.1) voor i = 0, . . . , n. We stellen tevens een formule op voor de restterm rn (x) = f (x) − pn (x) Onze formules zullen gebaseerd zijn op de stelling van Rolle (zie de cursus “Analyse”, volume 1): Als f : [a, b] → R een functie is die continu is over [a, b] en afleidbaar over (a, b), en f (a) = f (b), dan bestaat er een ξ ∈ (a, b) zodat f 0 (ξ) = 0

73

6.1

De formule van Lagrange

De elementaire veeltermen van Lagrange Om het probleem van Lagrange op te lossen kijken we eerst naar het speciaal geval waarin alle fi ’s nul zijn, op f j na. Onderstel f j = 1, met andere woorden fi = δi j voor i = 0, 1, . . . , n. De oplossing van het probleem van Lagrange noemen we L j . L j moet dus voldoen aan L j (xi ) = δi j (6.2) Stel L j (x) = k(x − x0 )(x − x1 ) · · · (x − x j−1 )(x − x j+1 ) · · · (x − xn ) Dan is voldaan aan (6.2) voor i 6= j. Om ook i = j in orde te brengen stellen we L j (x) =

n (x − x0 )(x − x1 ) · · · (x − x j−1 )(x − x j+1 ) · · · (x − xn ) x − xi = ∏ (x j − x0 )(x j − x1 ) · · · (x j − x j−1 )(x j − x j+1 ) · · · (x j − xn ) i=0, i6= j x j − xi

(6.3)

De interpolatieveelterm van Lagrange We keren terug naar het algemeen probleem (6.1). Stel pn (x) = f0 L0 (x) + f1 L1 (x) + · · · + fn Ln (x)

(6.4)

pn voldoet aan (6.1) en wordt de interpolerende veelterm van Lagrange genoemd. In de volgende stelling bepalen we de restterm. Stelling 6.1.1 Onderstel dat I een open interval is dat de spilpunten x0 , x1 , . . . , xn bevat. Onderstel dat f (n+1) bestaat en eindig is over I. Dan bestaat er een ξ ∈ I zodat f (x) = pn (x) +

(x − x0 )(x − x1 ) · · · (x − xn ) (n+1) f (ξ) (n + 1)!

(6.5)

Bewijs. De formule is juist als x = xi een der spilpunten is. Dan is namelijk de voorgestelde restterm nul, en f (xi ) = pn (xi ). Onderstel dus dat x 6= xi , voor i = 0, 1, . . . , n. Neem x vast, en beschouw de hulpfunctie g(t) = f (t) − pn (t) − ( f (x) − pn (x))

(t − x0 )(t − x1 ) · · · (t − xn ) (x − x0 )(x − x1 ) · · · (x − xn )

De nulpunten van g zijn x0 , x1 , . . . , xn en x. In het totaal zijn er dus n + 2 verschillende nulpunten. Vanwege de stelling van Rolle heeft g0 een nulpunt tussen elke koppel nulpunten van g. g0 heeft dus zeker n + 1 verschillende nulpunten. Op dezelfde manier heeft g00 zeker n verschillende nulpunten. 74

Als we deze redenering herhalen, dan vinden we dat g(n+1) minstens een nulpunt ξ ∈ I bevat. Nu is f (x) − pn (x) g(n+1) (t) = f (n+1) (t) − (n + 1)! (x − x0 )(x − x1 ) · · · (x − xn ) Stel t = ξ, en los op naar f (x). We vinden (6.5).



Voorbeeld 6.1.2 We beschouwen de functie 2 f (x) = √ π

Z x

2

e−t dt

0

Gegeven is f (1.0) = 0.84270 f (1.1) = 0.88021 f (1.2) = 0.91031 Gevraagd wordt f (1.15). We benaderen dit door p2 (1.15): (1.15 − 1.1)(1.15 − 1.2) (1.0 − 1.1)(1.0 − 1.2) (1.15 − 1.0)(1.15 − 1.1) (1.15 − 1.0)(1.15 − 1.2) + 0.91031 +0.88021 (1.1 − 1.0)(1.1 − 1.2) (1.2 − 1.1)(1.2 − 1.1) = 0.89619

f (1.15) ≈ p2 (1.15) = 0.84270

In dit geval kunnen we ook een bovengrens vinden voor de fout r(x). 2 2 f 0 (x) = √ e−x π 2 4x f 00 (x) = − √ e−x π 2 2 f 000 (x) = √ (4x2 − 2)e−x π

en

2 2 2 | f 000 (x)| ≤ max √ (4x2 − 2)e−x ≤ √ (4(1.2)2 − 2)e−1 ≈ 1.6 1≤x≤1.2 π π

zodat |r(1.15)| < |(1.15 − 1.0)(1.15 − 1.1)(1.15 − 1.2)|

1.6 ≈ 0.000098 6

Opmerking 6.1.3 In ons voorbeeld werd het aantal spilpunten opgelegd. In de praktijk wordt f dikwijls gegeven door een tabel, en dan kan men het aantal spilpunten zo groot kiezen als men wil. Men zal dan achtereenvolgens p1 (x), p2 (x), p3 (x), . . . berekenen, waarbij men de spilpunten kiest die het dichtst bij x liggen, en men stopt als de berekende waarden dicht bij elkaar liggen. In dit geval is het niet effici¨ent om (6.4) te gebruiken. Het is beter om een recursief algoritme op te stellen, en dit is wat we hierna zullen doen. 75

6.2

Het algoritme van Neville

We noteren pn1 ,n2 ,...,nk (x) voor de interpolatieveelterm van Lagrange die de spilpunten xn1 , xn2 , . . . , xnk gebruikt. Dit is een veelterm van graad k − 1. Zo is p3,6,9 de veelterm van graad 2 die de spilpunten x3 , x6 en x9 gebruikt. Lemma 6.2.1 x − xn 1 k pn1 ,n2 ,...,nk ,nk+1 (x) = xnk+1 − xnk x − xnk+1

pn1 ,n2 ,···,nk−1 ,nk (x) pn1 ,n2 ,···,nk−1 ,nk+1 (x)

(6.6)

Bewijs. Beide leden zijn veeltermen van graad k. Het volstaat om te bewijzen dat deze gelijk zijn in de k + 1 punten xn1 , xn2 , . . . , xnk+1 . Voor i = 1, 2, . . . , k − 1 vinden we xn − xn f (x ) n 1 i i k f (xni ) = xnk+1 − xnk xni − xnk+1 f (xni ) Bovendien geldt 0 1 f (xnk ) = xnk+1 − xnk xnk − xnk+1 en 1 f (xnk+1 ) = xnk+1 − xnk

xn − xn k+1 k 0

f (xnk ) ? f (xnk+1 ) ?



en dit bewijst het Lemma.

Lemma 6.2.1 wordt systematisch gebruikt bij het opstellen van de volgende tabel. Dit is het algoritme van Neville. De gezochte waarden staan op de diagonaal. x1 x2 x3 x4

x − x1 x − x2 x − x3 x − x4

p1 (x) p2 (x) p1,2 (x) p3 (x) p1,3 (x) p1,2,3 (x) p4 (x) p1,4 (x) p1,2,4 (x) p1,2,3,4 (x)

Hierbij is pi (x) = fi . Voorbeeld 6.2.2 Beschouw de situatie uit voorbeeld 6.1.2. Het algoritme van Neville geeft volgende tabel, voor x = 1.15. xi 1.1 1.2 1.0

x − xi 0.05 −0.05 0.15

fi 0.88201 0.91031 0.89526 0.84270 0.89897 0.89619

76

6.3

De methode der gedeelde verschillen

Het nadeel bij het algoritme van Neville is dat we pn (x) slechts in een punt berekenen. We zouden een formule willen opstellen die toelaat om pn te bepalen, en die zo is dat we niet helemaal opnieuw moeten beginnen met rekenen als we een spilpunt toevoegen. Om dit te verwezenlijken schrijven we pn (x) onder de vorm pn (x) = a0 + a1 (x − x0 ) + a2 (x − x0 )(x − x1 ) + · · · + an (x − x0 )(x − x1 ) · · · (x − xn−1 )

(6.7)

Als we hierin x = x0 , x1 , . . . , xn invullen, dan vinden we   f 0 = a0      f1 = a0 + a1 (x1 − x0 ) f2 = a0 + a1 (x2 − x0 ) + a2 (x2 − x0 )(x2 − x1 )  ..   .    f = a + a (x − x ) + a (x − x )(x − x ) + · · · + a (x − x )(x − x ) · · · (x − x ) n n n n n n 0 1 n 0 2 n 0 1 0 1 n−1 (6.8) Dit is een driehoekig lineair stelsel dat gemakkelijk op te lossen is. Omdat de punten xi onderling verschillen is de oplossing enig. Bovendien veranderen de berekende co¨effici¨enten a0 , . . . , an niet als men spilpunten toevoegt. Immers, dit komt erop neer aan het stelsel (6.8) een vergelijking toe te voegen, met nieuwe onbekende an+1 . Hieruit volgt ook dat a j enkel afhangt van de functiewaarden f0 , f1 , . . . , f j . Daarom voeren we de volgende notatie in a j = f [x0 , x1 , . . . , x j ] (6.9) Als we de co¨effici¨enten van xn vergelijken in (6.7) en n

n

n

i=0

i=0

pn (x) = ∑ fi Li (x) = ∑ fi

∏

(x − x j )

∏

(xi − x j )

j=0, j6=i n j=0, j6=i

dan vinden we

n

an = f [x0 , x1 , . . . , xn ] = ∑

i=0

fi n

∏

(6.10)

(xi − x j )

j=0, j6=i

Hieruit vinden we bijvoorbeeld dat f [x0 , x1 ] =

f (x1 ) − f (x0 ) x1 − x0

(6.11)

Uit (6.10) volgt eveneens dat an een lineaire combinatie is van f0 , f1 , . . . , fn , en dat an een symmetrische functie is van x0 , x1 , . . . , xn . Dit laatste betekent dat an niet verandert als we de volgorde van de xi ’s veranderen. We beweren nu dat an = f [x0 , x1 , . . . , xn ] =

f [x1 , x2 , . . . , xn ] − f [x0 , x1 , . . . , xn−1 ] xn − x0 77

(6.12)

Het grote voordeel aan (6.12) is dat we de co¨effici¨enten an nu recursief kunnen berekenen. We bewijzen (6.12) met behulp van (6.10). Noem het rechter- en linkerlid van (6.12) resp. RL en LL. Dan geldt RL =

n 1 ∑ xn − x0 i=1

fi

−

n

∏

(xi − x j )

1 n−1 ∑ xn − x0 i=0

j=1, j6=i

=

fn

∏

(xi − x j )

j=0, j6=i

f0

+

n−1

fi n−1

n

(x0 − x j ) ∏ (xn − x j ) ∏ j=1 j=0

+

1 n−1 ∑ xn − x0 i=1

fi



n−1

∏

(xi − x j )

1  1 − xi − xn xi − x0

j=1, j6=i

Nu is

1  1 1  1 = − xn − x0 xi − xn xi − x0 (xi − xn )(xi − x0 )

zodat

n

RL = ∑

i=0

fi n

∏

= LL

(xi − x j )

j=0, j6=i

De berekening van de co¨effici¨enten a j kan nu uitgevoerd worden met behulp van de volgende tabel x0 x1 x2 x3

f [x0 ] f [x1 ] f [x0 , x1 ] f [x2 ] f [x1 , x2 ] f [x0 , x1 , x2 ] f [x3 ] f [x2 , x3 ] f [x1 , x2 , x3 ] f [x0 , x1 , x2 , x3 ]

(6.13)

Als men een spilpunt toevoegt, dan verandert de tabel niet, maar wordt slechts uitgebreid met een nieuwe rij. Om de interpolerende veelterm in een punt uit te rekenen schrijft men deze als volgt op. 
 p3 (x) = a0 + (x − x0 ) a1 + (x − x1 ) a2 + (x − x2 )a3

6.4

Differentierekening

Tot nu toe waren de spilpunten x0 , x1 , . . . , xn willekeurig. In vele gevallen liggen de spilpunten op een gelijke afstand h van elkaar. We kunnen dan schrijven xk = x0 + kh 78

waarbij k = 0, 1, . . . , n. De overeenstemmende functiewaarden noteren we nog steeds door f0 , f1 , · · · , fn . Tabel (6.13) neemt nu de volgende vorm aan x0

f0

x1

f1

x2

f2

x3

f3

f1 − f0 h f2 − f1 h f3 − f2 h

f2 − 2 f1 + f0 2h2 f3 − 2 f2 + f1 2h2

(6.14) f3 − 3 f2 + 3 f1 − f0 6h3

We zullen alles nu in een meer elegante vorm herschrijven. We noteren RN voor de vectorruimte waarvan de elementen numerieke rijen ( f0 , f1 , f2 , · · ·) = ( fn ) zijn. We bekijken ook de operator (in feite de lineaire afbeelding) ∆ : RN −→RN gedefinieerd als volgt (∆ fn ) = ( f1 − f0 , f2 − f1 , f3 − f2 , . . .) = (∆ f1 , ∆ f2 , ∆ f3 , . . .) = ( fn+1 − fn ) Op een rij ( fn ) kunnen we de operator ∆ verschillende malen toepassen. We vinden bijvoorbeeld (∆(∆ fn )) = (∆2 fn ) = ( f2 − 2 f1 + f0 , f3 − 2 f2 + f1 , . . .) = ( fn+2 − 2 fn+1 + fn ) De tabel van de gedeelde verschillen (6.14) neemt nu de volgende vorm aan x0

f0

x1

f1

x2

f2

x3

f3

∆ f0 h ∆ f1 h ∆ f2 h

∆2 f 0 2h2 ∆2 f 1 2h2

(6.15) ∆3 f0 6h3

In het algemeen zien we dat

∆n f 0 (6.16) n!hn We voeren nu een verandering van veranderlijke uit, door x = x0 + sh te stellen. De algemene term van de interpolatieveelterm neemt volgende vorm aan. f [x0 , x1 , · · · , xn ] =

∆ j f0 sh(s − 1)h · · · (s − j + 1)h j!h j s(s − 1)(s − 2) · · · (s − j + 1) j = ∆ f0 j!   s j = ∆ f0 j

a j (x − x0 )(x − x1 ) · · · (x − x j−1 ) =

79

en we kunnen de interpolerende veelterm in de volgende elegante vorm herschrijven (de interpolatieveelterm van Newton). n   s j pn (x) = pn (x0 + sh) = ∑ ∆ f0 (6.17) j=0 j Behalve ∆ kunnen we nog andere lineaire transformaties van RN bekijken. We noteren I voor de identieke afbeelding. E is de shift operator: E : RN −→RN wordt gegeven door (E fn ) = ( f1 , f2 , f3 , · · ·) = ( fn+1 ) We zien onmiddellijk dat

∆ = E − I en E = ∆ + I

Doordat E en I commuteren (dit wil zeggen dat E ◦I = I ◦E), mogen we het binomium van Newton toepassen, en we vinden     k k−1 k k−2 k k k ∆ = (E − I) = E − E + E − · · · + (−1)k I (6.18) 1 2 Als we (6.18) toepassen op de rij ( fn ), dan vinden we     k k k fn+k−2 − · · · + (−1)k fn ∆ fn = fn+k − fn+k−1 + 2 1

(6.19)

Op analoge wijze vinden we       k k 2 k E = (I + ∆) = I + ∆+ ∆ +···+ ∆k−1 + ∆k 1 2 k−1

(6.20)

      k k 2 k fn+k = fn + ∆ fn + ∆ fn + · · · + ∆k−1 fn + ∆k fn 1 2 k−1

(6.21)

k

en

k

Als toepassing van (6.21) geven we volgende eigenschap van machtreeksen. Stelling 6.4.1 Beschouw een rij ( fn ). Als x ligt in het inwendige van het convergentieinterval van de machtreeksen ∞ ∞ fn n ∆n f 0 n x en x ∑ ∑ n=0 n! n=0 n! dan hebben we dat

∞

∞ fn n ∆n f0 n x x = e x ∑ ∑ n=0 n! n=0 n!

80

Bewijs. Machtreeksen mogen binnen hun convergentiegebeid vermenigvuldigd worden als gewone veeltermen. De co¨effici¨ent van xn in het product  ∞ xm  ∞ ∆s f  0 ∑ m! ∑ s! xs m=0 s=0 is

  1 ∆s f 0 1 n n s fn ∑ (n − s)! s! = n! ∑ s ∆ f0 = n! s=0 s=0 n



De laatste gelijkheid volgt uit (6.21). Voorbeeld 6.4.2 We stellen fn = In dit geval is

1 n+1

∞

∞ ex − 1 fn n xn x = = ∑ ∑ x n=0 n! n=0 (n + 1)!

zodat e−x

∞

∞ fn n 1 − e−x xn n x = = (−1) ∑ ∑ x (n + 1)! n=0 n! n=0

en uit stelling 6.4.1 volgt dat ∆n f0 =

(−1)n n+1

Als toepassing zullen we een - soms bruikbare - formule opstellen voor de n-de parti¨ele som van een machtreeks. We merken eerst op dat n−1

∑ fk = (I + E + E 2 + · · · + E n−1) f0

k=1

Nu is n−1

n−1

∆ ◦ ( ∑ E k ) = (E − I) ◦ ( ∑ E k ) k=0

k=0

= En − I = (∆ + I)n − I n   n i = ∑ ∆ i=1 i  n−1  n = ∆◦(∑ ∆k ) k + 1 k=0 Als ∆ : RN → RN bijectief zou zijn, dan zouden we kunnen besluiten dat  n−1 n−1  n k ∑ E = ∑ k + 1 ∆k k=0 k=0 81

(6.22)

∆ is geen bijectie, maar toch is (6.22) geldig. We bewijzen (6.22) per inductie op n. Voor n = 1 is de formule duidelijk geldig. Onderstel dat ze waar is voor n. Dan geldt  n n−1  n k ∑ E = ∑ k + 1 ∆k + (∆ + I)n k=0 k=0  n−1  n   n n k k = ∑ ∆ +∑ ∆ k=0 k + 1 k=0 k    n−1  n+1 k n+1 n = ∑ ∆ + ∆ n+1 k=0 k + 1  n  n+1 k = ∑ ∆ k=0 k + 1 en dit bewijst (6.22). Uit (6.22) volgt ook dat  n n   n−1 i−1 ∑ fk = ∑ i ∆ f0 i=1 k=0     n n 2 = n f0 + ∆ f0 + ∆ f0 + · · · + ∆n−1 f0 2 3

(6.23)

n

Als toepassing berekenen we

∑ k2. Stel in (6.23) fk = (k + 1)2. Dan is

k=0

f0 = 1, ∆ f0 = 3, ∆2 f0 = 2 en ∆k f0 = 0 voor k > 2. We vinden dat n−1

∑ (k + 1)2 =

k=0

n

∑ k2

= n+

k=1

n(n − 1)(n − 2) n(n − 1) 3+ 2 2 6

1 n(n + 1)(2n + 1) 6 Dezelfde methode kan gebruikt worden om de som van de r-de machten van de eerste n natuurlijke getallen op te stellen. We zien dat dit een veelterm van graad r + 1 in n is. =

6.5

De minimaxbenadering

De rest in de interpolatieformule van Lagrange wordt gegeven door de formule rn (x) = (x − x0 )(x − x1 ) · · · (x − xn )

f (n+1) (ξ) (n + 1)!

Op f (n+1) hebben we uiteraard geen enkele vat. In sommige gevallen kan men de spilpunten xi vrij kiezen, en dan gaat men die zo bepalen dat max{|(x − x0 )(x − x1 ) · · · (x − xn )| : x ∈ [a, b]} 82

minimaal is. In deze paragraaf bespreken we hoe we deze spilpunten kunnen bepalen. Door een lineaire transformatie kunnen we het interval [a, b] herleiden tot het interval [−1, 1]. Onze methode maakt gebruik van de veeltermen van Chebyshev. We zullen deze eerst behandelen. De veeltermen van Chebyshev De n-de veelterm van Chebyshev is per definitie Tn (x) = cos(nbgcos x)

(6.24)

(6.24) definieert Tn (x) enkel voor x ∈ [−1, 1]. We zullen eerst aantonen dat Tn (x) een veelterm is, die gedefinieerd is voor elke waarde van x. Stel x = cos θ, zodat Tn (x) = cos nθ. Uit de formule van Simpson cos(n + 1)θ + cos(n − 1)θ = 2 cos θ cos nθ volgt onmiddellijk dat Tn+1 (x) = 2xTn (x) − Tn−1 (x)

(6.25)

(6.25) laat toe om de veeltermen van Chebyshev recursief te berekenen. T0 (x) T1 (x) T2 (x) T3 (x) T4 (x) T5 (x) T6 (x)

= = = = = = =

1 x 2x2 − 1 4x3 − 3x 8x4 − 8x2 + 1 16x5 − 20x3 + 5x 32x6 − 48x4 + 18x2 − 1

We zien dus dat Tn een veelterm van graad n is. Stelling 6.5.1 De veeltermen van Chebyshev vormen een orthogonaal stel over [−1, 1] ten opzichte van de gewichtsfunctie 1 w(x) = √ 1 − x2 Bewijs. Neem n 6= m. We moeten aantonen dat Z 1

I=

cos(nbgcos x) cos(mbgcos x) √

−1

dx 1 − x2

Om dit te doen voeren we volgende substitutie uit. θ = bgcos x ; dθ = − √ 83

dx 1 − x2

=0

en vinden I = − = −

Z 0 π

1 2

cos nθ cos mθdθ

Z 0 π


cos(n + m)θ + cos(n − m)θ dθ

= 0  Uit (6.25) volgt dat de hoogste graadsterm van Tn (x) gelijk is aan 2n−1 . We zullen nu de grafiek van Tn (x) schetsen. Tn (x) = 0 als cos nθ = 0, of als nθ = (2k + 1)

π 2

met k ∈ Z. De verzameling nulpunten is dus    2k + 1  π : k = 0, 1, 2, . . . , n − 1 cos 2n Er zijn dus n enkelvoudige nulpunten, allen gelegen in het interval [−1, 1], zoals het hoort voor orthogonale veeltermen. We bepalen vervolgens de extreme waarden. d dθ sin nθ d Tn (x) = (cos nθ) = n dx dθ dx sin θ Als Tn0 (x) = 0, dan is dus sin nθ = 0, en dus is θ = kπ/n, waarbij k ∈ Z. k = 0 bepaalt geen nulpunt, en dus vinden we n − 1 nulpunten     k cos π : k = 1, 2, . . . , n − 1 n Deze liggen telkens tussen 2 nulpunten van Tn in, zodat in elk van deze punten een extremum bereikt wordt. De waarde die in dit extremum bereikt wordt is   k  cos nbgcos cos π = cos kπ = (−1)k n Laten we ook de waarden in de eindpunten bepalen. Tn (−1) = cos nπ = (−1)n en Tn (1) = cos 0 = 1 We kunnen besluiten dat Tn het interval [−1, 1] afbeeldt op [−1, 1]. De functiewaarden 1 en −1 worden bereikt in n + 1 punten, namelijk in de n − 1 nulpunten van de afgeleide, en ook in de punten 1 en −1. We kunnen dus de grafiek van Tn schetsen. We doen dit hieronder voor kleine waarden van n 84

y=T_3(x)

1.5

y=T_4(x)

1.5

1

1

0.5

y-as

y-as

0.5 0

0 -0.5 -0.5

-1

-1.5 -1.5

-1

-0.5

0

0.5

1

-1 -1.5

1.5

-1

-0.5

x-as y=T_5(x)

1.5

0

0.5

1

1.5

0.5

1

1.5

x-as y=T_6(x)

1.5

1

1

0.5

y-as

y-as

0.5 0

0 -0.5 -0.5

-1

-1.5 -1.5

-1

-0.5

0

0.5

1

-1 -1.5

1.5

-1

-0.5

x-as

0 x-as

Figuur 6.1: De veeltermen van Chebychef Oplossing van het minimaxprobleem Stelling 6.5.2 max{|(x − x0 )(x − x1 ) · · · (x − xn )| : x ∈ [−1, 1]} bereikt haar minimale waarde als men voor x0 , x1 , . . . , xn de n + 1 nulpunten van de Chebyshev veelterm Tn+1 kiest. Bewijs. De hoogste graadsterm van Tn+1 (x) is 2n . We moeten dus aantonen dat het minimum bereikt wordt door de veelterm (x − x0 )(x − x1 ) · · · (x − xn ) =

1 Tn+1 (x) 2n

Onderstel dat er een monische veelterm qn+1 (x) bestaat, verschillend van Tn+1 (x)/2n , en van graad n + 1, zodat   Tn+1 (x) 1 max{|qn+1 (x)| : x ∈ [−1, 1]} < max : x ∈ [−1, 1] = n (6.26) n 2 2 85

We stellen rn (x) = qn+1 (x) −

Tn+1 (x) 2n

rn is dan een veelterm van graad ten hoogste n. 1 In een punt waar Tn+1 2n zijn maximale waarde 2n bereikt, geldt dat rn (x) < 0, vanwege (6.26). −1 In een punt waar Tn+1 2n zijn minimale waarde 2n bereikt, geldt dat rn (x) > 0, eveneens vanwege (6.26). In de n + 2 punten waar Tn+1 2n extreem is, neemt rn dus afwisselend een positief en negatief teken aan. Dit betekent dat rn n + 1 nulpunten heeft. Aangezien rn van graad ten hoogste n is, betekent dit dat rn = 0. Dus is Tn+1 (x) qn+1 (x) = 2n en dit is strijdig met onze onderstelling, zodat onze stelling bewezen is. 

86

Hoofdstuk 7 Benadering volgens de methode van de kleinste kwadraten 7.1

Algemene methode

Een functie f is numeriek gegeven: de functiewaarde is bekend in de spilpunten x1 , x2 , . . . , xn . In het vorige hoofdstuk hebben we een veelterm van graad n − 1 geconstrueerd die in de spilpunten dezelfde waarde bereikte als f . Soms heeft dit niet veel zin. In vele gevallen is het zo dat de de functiewaarden in de spilpunten slechts in benadering bekend zijn. In dat geval heeft het geen zin om een veelterm te nemen die precies door de gegeven punten gaat. In andere gevallen heeft men bijkomende informatie over de functie f , bijvoorbeeld dat deze lineair is. We beschikken nu in feite over te veel gegevens: we hebben n gegevens, en er zijn slechts twee co¨effici¨enten te bepalen. In het algemeen gaan we als volgt te werk: we kiezen m + 1 functies φ0 (x), φ1 (x), . . . , φm (x) Gevraagd wordt om co¨effici¨enten a0 , a1 , . . . , am zo te bepalen dat f (x) zo goed mogelijk benaderd wordt door m

∑ a j φ j (x)

j=0

De meest gebruikte methode bestaat erin om te eisen dat 2 n  m f − a φ (x ) ∑ i ∑ j j i i=1

j=0

minimaal is. Nog iets algemener is de situatie waarbij niet alle benaderde functiewaarden fi even betrouwbaar zijn. In dat geval hecht men gewichten wi > 0 aan elke metingen, en men vraagt om de co¨effici¨enten a j zo te bepalen dat  2 n m H(a0 , a1 , . . . , am ) = ∑ wi fi − ∑ a j φ j (xi ) (7.1) i=1

87

j=0

minimaal is. Een nodige voorwaarde om een extremum te hebben is dat de parti¨ele afgeleiden van H naar elk van de ak nul zijn, dit wil zeggen dat n

∑ wi



i=1

 m fi − ∑ a j φ j (xi ) φk (xi ) = 0

(7.2)

j=0

(7.2) is een lineair stelsel van m + 1 vergelijkingen met m + 1 onbekenden a0 , a1 , . . . , am . Als de determinant van dit stelsel niet nul is, dan heeft het stelsel een unieke oplossing, en men kan aantonen dat hiermee een minimum overeenstemt. Men noemt (7.2) het stelsels normaalvergelijkingen. Benadering door veeltermen In principe kan men voor φ j (x) gelijk welke functie nemen. Een voor de hand liggende keuze is φ j (x) = x j Dit betekent dat men f tracht te benaderen door een veelterm van graad m. De normaalvergelijkingen worden nu   n m j k w f − a x ∑ i i ∑ j i xi = 0 i=1

of

m

∑



j=0

n

∑ wixi

j+k

j=0 i=1



n

a j = ∑ wi fi xik

(7.3)

i=1

voor k = 0, 1, . . . , m. Alles gaat goed zolang m niet te groot wordt. Voor m ≥ 6 hebben we niet alleen het probleem dat we een lineair stelsel van relatief grote omvang moeten oplossen, het blijkt bovendien dat dit stelsel meestal slecht geconditionneerd is. Een belangrijk geval is het geval n = 1: men benadert de grafiek van f door een rechte. Men spreekt dan van lineaire regressieanalyse. In veel praktische gevallen weet men ook dat f (0) = 0. Het stelsel (7.3) herleidt zich dan tot 1 vergelijking met 1 onbekende, namelijk n

n

i=1

i=1

a1 ∑ wi xi2 = ∑ wi fi xi en

7.2

∑ni=1 wi fi xi a1 = n ∑i=1 wi xi2

Benadering met Forsythe veeltermen

Zoals al gezegd is het lineaire stelsel (7.3) slecht geconditionneerd. Als we ervoor kunnen zorgen dat het stelsel (7.3) een diagonaal stelsel wordt, dan is er uiteraard geen probleem. Om dit te

88

verwezenlijken zullen we niet langer werken met φ j (x) = x j , maar we zullen voor φ j (x) = p j (x) een veelterm van graad j nemen. Voor twee functies f en g defini¨eren we n

h f , gi = ∑ wi fi gi

(7.4)

i=1

Merk eerst op dat h•, •i alle eigenschappen van een inwendig product heeft: h•, •i is bilineair, symmetrisch en niet-negatief. Alleen de positiviteit hebben we niet: het is mogelijk dat h f , f i = 0, terwijl de functie f toch niet identiek nul is: het volstaat om een functie f te nemen die 0 is in elk van de spilpunten x1 , x2 , . . . , xn . Het lineair stelsel (7.2) kunnen we nu herschrijven onder de vorm m

∑ hp j , pk ib j = h f , pk i

(7.5)

j=0

voor j = 0, 1, . . . , m. Het stelsel (7.5) is een diagonaalstelsel indien hp j , pk i = 0 voor j 6= k. Vandaar de volgende definitie. Definitie 7.2.1 Zij p0 , p1 , . . . , pm een stel veeltermen van graad 0, 1, 2, . . . , m. We noemen dit stel orthogonaal over de spilpunten x1 , x2 , . . . , xn ten opzichte van de gewichten wi indien hp j , pk i = 0 voor j 6= k. De veeltermen p j worden ook Forsythe veeltermen genoemd. Constructie van de Forsythe veeltermen Een eerste manier bestaat erin om het orthogonalizatieproc´ed´e van Gram-Schmidt toe te passen op 1, x, x2 , . . . , xm . In de hiernavolgende stelling geven we een meer praktische formule. Stelling 7.2.2 Stel p−1 (x) = 0 en p0 (x) = 1. Er bestaan getallen α j+1 en β j ( j = 0, 1, . . . , n) zodat p j+1 (x) = (x − α j+1 )p j (x) − β j p j−1 (x) orthogonaal zijn over de spilpunten x1 , x2 , . . . , xn ten opzichte van de gewichten wi . Bewijs. We werken per inductie op j. Voor j = 0 wordt (7.6) p1 (x) = (x − α1 )p0 Dan is

hp1 , p0 i = hx − α1 , 1i = 0

als

hx, 1i ∑ni=1 wi xi = n h1, 1i ∑i=1 wi Onderstel nu dat p0 , p1 , . . . , pk gekend zijn en orthogonaal zijn. We stellen, volgens (7.6) α1 =

pk+1 (x) = (x − αk+1 )pk (x) − βk pk−1 (x) 89

(7.6)

We moeten αk+1 en βk zo bepalen dat hpk+1 , p j i = hxpk , p j i − αk+1 hpk , p j i − βk hpk−1 , p j i = 0

(7.7)

voor j = 0, 1, . . . , k. Omdat hxpk , p j i = hpk , xp j i is aan (7.7) altijd voldaan voor j = 0, 1, . . . , k − 2. Voor j = k − 1 herleidt (7.7) zich tot hpk , xpk−1 i − βk hpk−1 , pk−1 i = 0 en hieraan is voldaan als βk =

hpk , xpk−1 i ∑n wi xi pk (xi )pk−1 (xi ) = i=1 n hpk−1 , pk−1 i ∑i=1 wi pk−1 (xi )2

(7.8)

We kunnen de formule voor βk nog iets eenvoudiger herschrijven. Stel in (7.6) j = k − 1 en neem het inwendig product met pk . Dan volgt hpk , pk i = hpk , xpk−1 i − αk hpk , pk−1 i − βk−1 hpk , pk−2 i = hpk , xpk−1 i en

hpk , pk i ∑ni=1 wi pk (xi )2 βk = = hpk−1 , pk−1 i ∑ni=1 wi pk−1 (xi )2

(7.9)

Voor j = k vinden we voor (7.7) hpk , xpk i − αk+1 hpk , pk i = 0 en hieraan is voldaan als αk+1 =

hpk , xpk i ∑ni=1 wi xi pk (xi )2 = n hpk , pk i ∑i=1 wi pk (xi )2

(7.10) 

Opmerkingen 7.2.3 1) De stelling geeft niet alleen het bestaan van, maar ook een formule voor α j+1 en β j : (7.9) en (7.10) laten toe om α j+1 en β j te berekenen. 2) De orthogonale veeltermen zijn slechts bepaald op een constante na. In het bewijs van stelling 7.2.2 worden ze zo geconstrueerd dat de co¨effici¨ent van de hoogste graadsterm 1 is. 3) Uiteraard mag de noemer n

∑ wi pk (xi)2

i=1

van (7.9) en (7.10) niet nul zijn. Zolang k < n is dit geen probleem, want een veelterm van graad k < n heeft ten hoogste n − 1 nulpunten, zodat zeker e´ e´ n van de pk (xi ) 6= 0, en dan is de noemer verschillend van nul. Vanaf k = n is er wel een probleem. Om dit in te zien merken we op dat we in feite in de ndimensionale Euclidische ruimte aan het werken zijn: voor ons inwendig product zijn twee functies f en g aan elkaar gelijk van zodra ze in de n spilpunten gelijk zijn. In een n-dimensionale Euclidische ruimte is het maximaal aantal ortogonale, van nul verschillende vectoren n. Als we dus p0 , . . . , pn−1 geconstrueerd hebben, dan zijn we dus aan het einde van onze mogelijkheden: 90

elke vector die loodrecht staat op p0 , . . . , pn−1 is nul. De veelterm pn (x) die men uiteindelijk vindt door (7.6) herhaaldelijk toe te passen is dus noodzakelijk nul in elk van de spilpunten. Aangezien gr(pn ) = n geldt dus noodzakelijk dat pn (x) = (x − x1 )(x − x2 ) . . . (x − xn ) Bij de volgende iteratie is de noemer van (7.10) nul, en het algoritme is niet langer van toepassing. Voor ons heeft dit echter geen belang, aangezien we enkel p0 , p1 , . . . , pn−1 nodig hebben. Voorbeeld 7.2.4 We stellen wi = 1, en werken met de spilpunten 0, ±1, ±2. We vinden achtereenvolgens p0 (x) = 1 α1 =

1 5 ∑ xi = 0 5 i=1

p1 (x) = x ∑5i=1 p1 (xi )2 =2 ∑5i=1 p0 (xi )2 ∑5i=1 xi3 = =0 ∑5i=1 xi2

β1 = α2

p2 (x) = xp1 (x) − 2p0 (x) = x2 − 2 ∑5i=1 (xi2 − 2)2 7 β2 = = 5 ∑5i=1 xi2 ∑5i=1 xi p2 (xi )2 α3 = =0 ∑5i=1 p2 (xi )2 7 1 p3 (x) = xp2 (x) − p1 (x) = (5x3 − 17x) 5 5 5 2 36 ∑i=1 p3 (xi ) β3 = = 5 2 35 ∑i=1 p2 (xi ) ∑5i=1 xi p3 (xi )2 =0 ∑5i=1 p3 (xi )2 36 1 p4 (x) = xp3 (x) − p2 (x) = (35x4 − 155x2 + 72) 35 35 5 2 4 ∑i=1 p4 (xi ) β4 = = 5 2 7 ∑i=1 p3 (xi ) α4 =

∑5i=1 xi p4 (xi )2 =0 ∑5i=1 p4 (xi )2 4 p5 (x) = xp4 (x) − p3 (x) = x5 − 5x3 + 4x = x(x − 1)(x + 1)(x − 2)(x + 2) 7 Voorbeeld 7.2.5 We stellen wi = 1, en nemen als spilpunten de nulpunten van de n-de Chebyshev veelterm Tn (x). Als orthogonale veeltermen vinden we dan precies T0 (x), T1 (x), . . . , Tn (x). α5 =

91

Voor we dit bewijzen formuleren we eerst enkele eenvoudige Lemmas. Lemma 7.2.6 Onderstel dat de getallen eiα0 , eiα1 , . . . , eiαn−1 in het complexe vlak de hoekpunten zijn van een regelmatige n-hoek. Dan is n−1

∑ cos α j = 0

j=0

Bewijs. Stel α0 = α. Dan zijn de hoeken die optreden {α +

2kπ | k = 0, . . . , n − 1} n

De eiα j zijn dan de nulpunten van de veelterm P(z) = (e−iα z)n − 1 De som van de nulpunten is min de co¨effici¨ent van zn−1 , en is dus 0. Daarom is n−1

∑ eiα j = 0

j=0

en het resultaat volgt na het nemen van het re¨eel gedeelte van beide leden.



Lemma 7.2.7 Neem k ∈ {1, 2, . . . , 2n − 1}. Dan is n−1

∑ cos(

j=0

2j+1 kπ) = 0 2n

Bewijs. Omdat cos x = cos(2π − x) is n−1

∑ cos(

j=0

2j+1 1 2n−1 2j+1 kπ) = ∑ cos( kπ) 2n 2 j=0 2n

De punten {ei

2 j+1 2n π

| j = 0, 1, . . . , 2n − 1}

zijn de hoekpunten van een regelmatige 2n-hoek. Ook de punten {ei

2 j+1 2n kπ

| j = 0, 1, . . . , 2n − 1}

zijn de hoekpunten van een regelmatige veelhoek. Het aantal hoekpunten is een deler van 2n. Het aantal hoekpunten is ook groter dan 1, omdat k < 2n. Vanwege lemma 7.2.6 is dus 2n−1

∑

j=0

cos(

2j+1 kπ) = 0 2n 92



en dit bewijst ons Lemma. Bewijs. van voorbeeld 7.2.5 Herhaal dat Tn (x) = cos(nbgcos x). De n-nulpunten x0 , x1 , . . . , xn−1 worden gegeven door de formule xi = cos

2i + 1 π 2n

en we vinden, voor 0 ≤ l < k < n, n−1

hTk , Tl i =

∑ Tk (xi)Tl (xi)

i=0 n−1

=

∑ cos

i=0 n−1

(2i + 1)l (2i + 1)k π cos π 2n 2n

1 (2i + 1) 1 n−1 (2i + 1) cos (k + l)π + cos (k − l)π ∑ ∑ 2 i=0 2n 2 i=0 2n = 0 =

Hierbij gebruikten we de formules van Simpson, en lemma 7.2.7 (in acht genomen dat 0 < k + l < 2n en 0 < k − l < 2n).  Zonder bewijs formuleren we ook nog volgend voorbeeld. Voorbeeld 7.2.8 Stel wi = 1, en neem als spilpunten 0, 1, 2, . . . , n. Dan worden de pk (x) gegeven door de formule   pk (x) = ∆k x(x − 1) · · · (x − k + 1)(x − n − 1)(x − n − 2) · · · (x − n − k)

93

Hoofdstuk 8 Numeriek Afleiden 8.1

Afleiden van de interpolatieformule van Lagrange

Gegeven is een functie y = f (x). We zullen onderstellen dat f gekend is in n + 1 spilpunten x0 , x1 , . . . , xn . We zullen bovendien onderstellen dat de spilpunten ekwidistant zijn: xk = x0 + kh Deze onderstelling is strikt genomen niet nodig, maar vereenvoudigt wel de formules. In de praktijk gebeurt het dikwijls dat een functie f numeriek gegeven wordt door een tabel. Gevraagd wordt om f 0 te berekenen in elk van de spilpunten. Men kan daarna f 0 benaderen in een willekeurig punt door interpolatie. We onderstellen dat de spilpunten x0 , x1 , . . . , xn in een open interval I gelegen zijn, en dat f (n+1) (x) eindig is over I. We kunnen dan de interpolatieformule van Lagrange (zie stelling 6.1.1) opschrijven: voor elke x ∈ I hebben we n

n

i=0

j=0

f (x) = ∑ fi Li (x) + ∏ (x − x j )

f (n+1) (ξ(x)) (n + 1)!

(8.1)

We leiden (8.1) af en vinden n

f (x) = ∑ 0

fi Li0 (x)

+

i=0

 f (n+1) (ξ(x)) d  n (x − x j ) (n + 1)! dx ∏ j=0 n

+

∏ (x − x j ) j=0

d f (n+1) (ξ(x)) dx (n + 1)!

De laatste term is moeilijk te manipuleren, omdat we niets weten over de functie ξ(x). Als we onderstellen dat de afgeleide van f (n+1) ◦ ξ eindig is over I, dan is de laatste term nul in elk der spilpunten. Dit is de reden waarom we de formule enkel gebruiken in de spilpunten. In een spilpunt

94

xk vinden we  d n (x − x ) = j dx ∏ x=xk j=0

n

∏

(xk − x j )

j=0, j6=k n

= h k(k − 1)(k − 2) · · · 2.1.(−1)(−2) · · · (−(n − k)) = hn (−1)n−k k!(n − k)! We krijgen dus volgende formule voor de afgeleide in de spilpunten. n

f (xk ) = ∑ fi Li0 (xk ) + (−1)n−k hn 0

i=0

8.2

k!(n − k)! (n+1) f (ξk ) (n + 1)!

(8.2)

Praktische formules

Na enig rekenwerk kunnen we (8.2) expliciet uitschrijven in de gevallen n = 2, 3, . . .. We doen dit voor n = 2. In dit geval hebben we (x − x1 )(x − x2 ) 1 = 2 (x − x0 − h)(x − x0 − 2h) (x0 − x1 )(x0 − x2 ) 2h 1 (x − x0 )(x − x2 ) = − 2 (x − x0 )(x − x0 − 2h) L1 (x) = (x1 − x0 )(x1 − x2 ) h 1 (x − x0 )(x − x1 ) = 2 (x − x0 )(x − x0 − h) L2 (x) = (x2 − x0 )(x2 − x1 ) 2h

L0 (x) =

Afleiden geeft 1 (2x − 2x0 − 3h) 2h2 1 L10 (x) = − 2 (2x − 2x0 − 2h) h 1 L20 (x) = (2x − 2x0 − h) 2h2 L00 (x) =

Substitutie in (8.2) geeft, na enig gereken 1 h2 (−3 f0 + 4 f1 − f2 ) + f (3) (ξ0 ) 2h 3 2 1 h f 0 (x1 ) = (− f0 + f2 ) − f (3) (ξ1 ) 2h 6 1 h2 f 0 (x2 ) = ( f0 − 4 f1 + 3 f2 ) + f (3) (ξ2 ) 2h 3

f 0 (x0 ) =

(8.3) (8.4) (8.5)

Opmerkingen 8.2.1 1) In de praktijk beschikt men gewoonlijk over veel meer dan drie spilpunten, zodat men kan kiezen tussen de drie bovenstaande formules. Gewoonlijk kiest men dan (8.4), 95

omdat de fout tweemaal kleiner is, en omdat de formule iets eenvoudiger is dan de twee anderen. We merken ook op dat bovenstaande formules numeriek instabiel zijn: in (8.4) worden twee getallen die ongeveer aan mekaar gelijk zijn van mekaar afgetrokken. Dit is het geval voor alle numerieke formules om afgeleiden te berekenen. 2) Als we n = 1 stellen in (8.2), dan vinden we de volgende formule: h 1 f 0 (x0 ) = ( f1 − f0 ) − f 00 (ξ0 ) h 2 Vergelijk deze met de definitie van de afgeleide: 1 f 0 (x0 ) = lim ( f (x0 + h) − f (x0 )) h→0 h Als we (8.4) uitschrijven voor n = 2, 3, 4, dan vinden we de hiernavolgende formules. Deze zijn uitgeschreven in “moleculaire vorm”. Dit betekent bijvoorbeeld dat de eerste formule in het geval n = 3 kan herschreven worden als f 0 (x0 ) = -3

h3 1 (−11 f0 + 18 f1 − 9 f2 + 2 f3 ) − f (4) (ξ0 ) 6h 4

4 -1 —— ——

h2 3 3D

2hD

N

rest

2hD

-1 0 1 N

—— ——

rest − h6 D3

2hD

1 -4 3 N

—— ——

rest

6hD

-11 18 -9 2 N —— —— ——

rest − h4 D4

6hD

-2 -3 6 -1 N

—— —— ——

rest

6hD

1 -6 3 2 N

—— —— ——

h rest − 12 D4

6hD

-2 9 -18 11 N

—— —— ——

rest

h3 4 4D

12hD

-25 48 -36 16 -3 N —— —— —— ——

rest

h4 5 5D

2

h2 3 3D

3

h3 4 12 D

3

96

12hD

-3 -10 18 -6 1 N

—— —— —— ——

h rest − 20 D5

12hD

1 -8 0 8 -1 N

—— —— —— ——

h rest − 30 D5

12hD

-1 6 -18 10 3 N

—— —— —— ——

h rest − 20 D5

12hD

3 -16 36 -48 25 N

—— —— —— ——

rest

97

4

4

4

h4 5 5D

Hoofdstuk 9 Numerieke Integratie 9.1

De Newton-Cotes formules

Ons doel is om algoritmen te ontwikkelen die toelaten om Z b

f (x)dx a

numeriek te berekenen. Onderstel dat a, b, x0 , . . . , xn gelegen zijn in een gesloten interval I waarover f tenminste n + 1 maal continu differentieerbaar is. We onderstellen verder dat f gekend is in de spilpunten x0 , x1 , . . . , xn . Net als in het voorgaande hoofdstuk vertrekken we van de interpolatieformule van Lagrange. n

f (x) = ∑ fi Li (x) + rn (x)

(9.1)

i=0

met

f (n+1) (ξ(x)) n (x − x j ) rn (x) = (n + 1)! ∏ j=0

(9.1) integreren geeft Z b a

We stellen ci =

Rb a

n

Z b

i=0

a

f (x)dx = ∑ fi

Z b

Li (x)dx +

a

rn (x)dx

Li (x)dx, en vinden als integratieformule Z b a

n

f (x)dx = ∑ ci fi + E( f )

(9.2)

i=0

waarbij de gemaakte fout E( f ) gegeven wordt door de formule Z b

E( f ) = a

98

rn (x)dx

(9.3)

Opmerkingen 9.1.1 1) Als f een veelterm van graad ten hoogste n is, dan is de fout rn (x) in de interpolatieformule 0, en dus is ook E( f ) = 0. We zeggen dat de interpolatieformule juist is. 2) E is linear als functie van f : E(α f + βg) = αE( f ) + βE(g) We maken nu welbepaalde keuzes voor de spilpunten. Stel h = (b − a)/n, a = x0 , xk = x0 + kh en b = xn . a en b zijn dus spilpunten, en de spilpunten zijn ekwidistant. De integratieformules die h

a = x0

x1

b x2

x n–1

xn

Figuur 9.1: De gesloten Newton-Cotes formules men aldus verkrijgt noemt men de gesloten Newton-Cotes formules. Bij de open Newton-Cotes formules gebruikt men a en b niet als spilpunten. Men stelt h = (b − a)/(n + 2), x−1 = a, xk = x−1 + (k + 1)h, en dus b = xn+1 . In beide gevallen vinden we voor de a

x1

x0

xn

b

Figuur 9.2: De open Newton-Cotes formules co¨effici¨enten ci de volgende formule Z b

ci =

a

n

x−xj (−1)n−i dx = ∏ hn i!(n − i)! j=0, j6=i xi − x j

Z b

n

∏

a j=0, j6=i

(x − x j )dx

(9.4)

We zullen nu de co¨effici¨enten ci berekenen in een aantal eenvoudige gevallen. Voorbeeld 9.1.2 De open formule met 1 spilpunt In dit geval is x0 = (a + b)/2, en L0 = 1. (9.4) geeft c0 = b − a. De integratieformule neemt volgende vorm aan Z b

f (x)dx ≈ (b − a) f ((a + b)/2)

a

Dit is natuurlijk een erg ruwe benadering! In de praktijk gaat men gewoonlijk als volgt tewerk: men deelt het interval [a, b] op in een aantal deelintervallen, en past op elk van die deelintervallen de Newton-Cotes formule toe. Dit komt er in dit geval op neer de integraal te benaderen door een Riemann som. Voorbeeld 9.1.3 De trapeziumregel We schrijven de gesloten formule met twee spilpunten neer. Gebruik makend van (9.4) vinden we 1 c0 = − h

Z x1 x0

1 (x − x1 )dx = − h 99

Z 0

tdt = −h

h 2

en c1 = zodat

Z x1 x0

1 h

Z x1 x0

(x − x0 )dx =

h 2

h f (x)dx ≈ ( f0 + f1 ) 2

(9.5)

We berekenen dus in feite de oppervlakte van een trapezium. Als we het interval [a, b] opdelen in n

y y = f(x)

x0

x1 x Figuur 9.3: De trapeziumregel gelijke deelintervallen [xi−1 , xi ] en op elk van die deelintervallen de trapeziumregel toepassen, dan vinden we Z xn h f (x)dx ≈ ( f0 + 2 f1 + 2 f2 + . . . + 2 fn−1 + fn ) (9.6) 2 x0 De trapeziumregel wordt meestal in deze vorm gebruikt. Voorbeeld 9.1.4 De formule van Simpson We schrijven de gesloten formule met drie spilpunten neer. Om de berekeningen te vereenvoudigen onderstellen we eerst dat de spilpunten 0, 1 en 2 zijn. Uit (9.4) volgt ci =

(−1)i i!(2 − i)!

Z 2

2

∏

(x − j)dx

0 j=0, j6=i

en dus c0

1 = 2

zodat

Z 2

(x − 1)(x − 2)dx =

0

Z 2

1 3

4 3 0 Z 2 1 1 x(x − 1)dx = = 2 0 3

c1 = − c2

Z 2

x(x − 2)dx =

 1 g(0) + 4g(1) + g(2) 3 0 Om de algemene formule te vinden nemen we als substitutie g(t)dt ≈

x = x0 + th en dx = hdt 100

We krijgen Z x2

Z 2

f (x)dx = h x0

0

≈ of

Z x2 x0

f (x0 + ht)dt

 h f (x0 ) + 4 f (x0 + h) + f (x0 + 2h) 3 h f (x)dx ≈ ( f0 + 4 f1 + f2 ) 3

(9.7)

Dit is de formule van Simpson, en het is een van de veelgebruikte formules om bepaalde integralen te benaderen. Als we het interval [a, b] opdelen in n intervallen van gelijke lengte, en op elk van die intervallen de formule van Simpson toepassen, dan vinden we Z x2n x0

h f (x)dx ≈ ( f0 + 4 f1 + 2 f2 + 4 f3 + 2 f4 + · · · + 2 f2n−2 + 4 f2n−1 + f2n ) 3

(9.8)

Op deze manier kan men natuurlijk zoveel integratieformules opstellen als men wil. We geven enkele voorbeelden van eenvoudige formules. We vermelden ook de restterm, die we zullen bespreken in de volgende paragraaf. Gesloten formules Z x1

h3 h ( f0 + f1 ) − f 00 (ξ) 2 12 h5 h ( f0 + 4 f1 + f2 ) − f (4) (ξ) f (x)dx = 3 90 3h 3h5 (4) f (x)dx = ( f0 + 3 f1 + 3 f2 + f3 ) − f (ξ) 8 80 2h 8h7 (6) f (x)dx ≈ (7 f0 + 32 f1 + 12 f2 + 32 f3 + 7 f4 ) − f (ξ) 45 945 f (x)dx =

x0 Z x2 x0 Z x3 x0 Z x4 x0

(9.9) (9.10) (9.11) (9.12)

Open formules Z x1 x−1 Z x2 x−1 Z x3 x−1 Z x4 x−1

h3 00 f (ξ) 3 3h 3h3 00 f (x)dx = ( f0 + f1 ) + f (ξ) 2 4 4h 28h5 (4) f (x)dx = (2 f0 − f1 + 2 f2 ) + f (ξ) 3 90 5h 95h5 (4) f (x)dx = (11 f0 + f1 + f2 + 11 f3 ) + f (ξ) 24 144 f (x)dx = 2h f0 +

(9.15) staat bekend als de formule van Milne. 101

(9.13) (9.14) (9.15) (9.16)

Voorbeeld 9.1.5 We berekenen het getal π met de formule π=

Z 2p

4 − x2 dx

0

We gebruiken de formule van Simpson met h = 0.2, en vinden

1, 6

A B 0

C 1, 2

2

Figuur 9.4: Berekening van π

π≈

0.2 ( f0 + 4 f1 + 2 f2 + 4 f3 + 2 f4 + · · · + 4 f9 + f10 ) = 3.127008 3

Deze benadering voor π is een beetje teleurstellend. De reden voor de minder goede benadering is de verticale raaklijn in het punt x = 2. In een omgeving van x = 2 zijn de afgeleiden van orde 1, 2, . . . willekeurig groot, zodat ook de fout in de interpolatieformule van Lagrange groot wordt. Het is niet moeilijk om dit ook kwalitatief in te zien: bij de interpolatieformule van Lagrange benaderen we een functie f door een veelterm. Als deze functie f ergens een verticale raaklijn heeft, dan kan deze benadering nooit goed zijn, omdat veeltermen nooit een verticale raaklijn hebben. Het is mogelijk deze moeilijkheden te vermijden door de kwartcirkel te verdelen in drie delen A, B en C, zoals op bovenstaande figuur. De formule van Simpson geeft A+B ≈

0.2 ( f0 + 4 f1 + 2 f2 + 4 f3 + 2 f4 + 4 f5 + f6 ) = 2.246991 3

en B +C =

Z 1.6 p

4 − y2 dy ≈

0

0.2 ( f0 + 4 f1 + 2 f2 + 4 f3 + 2 f4 + 4 f5 + 2 f6 + 4 f7 + f8 ) = 2.814534 3

Tenslotte is π = (A + B) + (B +C) − B ≈ 2.246991 + 2.814534 − (1.2)(1.6) = 3.141525 en dit is een veel betere benadering dan de voorgaande.

102

9.2

De rest in de formule van Newton-Cotes

We beschouwen een Newton-Cotes formule met n + 1 spilpunten. We weten dat E( f ) = 0 als f een veelterm is van graad ten hoogste n. Het kan zijn dat ook E(xn+1 ) = 0. In dat geval is E( f ) = 0 voor veeltermen van graad ten hoogste n + 1. We defini¨eren het getal k als het maximale getal waarvoor E(xi ) = 0, voor elke i ≤ k. Het getal k kan eenvoudig bepaald worden: men bepaalt gewoon achtereenvolgens E(xn+1 ), E(xn+2 ), enz, en men stopt wanneer men een resultaat vindt dat niet gelijk is aan nul. Hoe groter k, hoe beter onze interpolatieformule! Onderstel dat de spilpunten 0, 1, . . . , n zijn. Als we ons inspireren op de rest bij de formule van Lagrange, dan verwachten we dat E( f ) van de volgende vorm is: E( f ) = c

f (k+1) (ξ) (k + 1)!

Voor een willekeurig interval vinden we dan, via de subsitutie x = a + th, dat Z b

Z n

f (x)dx = h a

f (a + th)dt 0

n

=

∑ ci fi + hchk+1

i=0

f (k+1) (ξ) (k + 1)!

en vandaar de volgende definitie: Definitie 9.2.1 Een integratieformule wordt een simplexformule genoemd indien men de rest kan schrijven onder de vorm (k+1) (ξ) k+2 f E( f ) = ch (9.17) (k + 1)! Zonder bewijs geven we volgende stelling Stelling 9.2.2 De Newton-Cotesformules zijn simplexformules. Voor een gegeven Newton-Cotes formule kunnen we nu gemakkelijk h en c bepalen. We geven twee voorbeelden. De rest in de formule van het trapezium We bepalen eerst het getal k. We weten dat k ≥ 1. Onderstel eerst dat het integratieinterval [0, 1] is. Z 1 1 1 1 1 E(x2 ) = x2 dx − (02 + 12 ) = − = − 6= 0 2 3 2 6 0 Daarom is k = 1. Aangezien 1 c d2 2 E(x2 ) = − = (x ) = c 6 2! dx2 103

vinden we dat c = −1/6. De formule voor de rest is E( f ) = −

1 00 f (ξ) 12

Voor een willekeurig integratieinterval krijgen we de volgende formule voor de rest E( f ) = −

h3 00 f (ξ) 12

(9.18)

Men ziet dit in door de substitutie x = a + ht toe te passen. De rest in de formule van Simpson We bepalen eerst het getal k. We weten dat k ≥ 2. Onderstel eerst dat het integratieinterval [0, 2] is. Z 2

1 x3 dx − (03 + 4 · 13 + 23 ) = 0 3 0 Z 2 1 4 E(x4 ) = x4 dx − (03 + 4 · 14 + 24 ) = − 3 15 0 3

E(x ) =

We concluderen dat k = 3 en c = −4/15. De formule voor de rest is E( f ) = −

h5 (4) f (ξ) 90

(9.19)

We hebben dus een voorbeeld van een integratieformule waarvoor k > n. Dit verklaart ook het belang van de formule van Simpson. Enerzijds is deze tamelijk eenvoudig, en langs de andere kant is deze veel nauwkeuriger dan bijvoorbeeld de trapeziumregel.

9.3

Rombergintegratie

We vertrekken van de trapeziumregel met rest Z x1 x0

h h3 f (x)dx = ( f0 + f1 ) − f 00 (ξ) 2 12

Onderstel dat f 00 (x) in benadering constant is, en stel c = − f 00 (x)/12. We passen de regel van het trapezium toe op de intervallen [x0 , x2 ], [x0 , x1 ] en [x1 , x2 ]. We krijgen enerzijds Z x2

I=

f (x)dx = x0

2h ( f0 + f2 ) + 8h3 c 2

en anderzijds Z x2

Z x1

f (x)dx = x0

Z x2

f (x)dx + x0

x1

h f (x)dx = ( f0 + 2 f1 + f2 ) + h3 c + h3 c 2 104

Stel

h A0,0 = h( f0 + f2 ) en A1,0 = ( f0 + 2 f1 + f2 ) 2

dan hebben we in feite Z x2 x0

f (x)dx = A0,0 + 8h3 c = A1,0 + 2h3 c

Hieruit kunnen we c elimineren, en we vinden Z x2 x0

h 1 f (x)dx = (4A1,0 − A0,0 ) = ( f0 + 4 f1 + f2 ) 3 3

Dit is niets anders dan de regel van Simpson, en is dus, zoals verwacht, een betere benadering voor de integraal. Meer algemeen kunnen we de integraal berekenen met behulp van de trapeziumformule door het integratieinterval op te delen in achtereenvolgens 1, 2, 4, 8, 16, . . . deelintervallen. De opeenvolgende benaderingen die we verkrijgen voor de integraal noemen we A0,0 , A1,0 , A2,0 , . . . , An,0 , . . .. Uit An−1,0 en An,0 vinden we dan een nog betere benadering, namelijk 1 An,1 = (4An,0 − An−1,0 ) 3 Op deze manier vinden we een nieuwe rij benaderingen. Uitgaande van de rij (An,1 ) kan men alweer een nieuwe rij nog betere benaderingen (An,2 ) opstellen. Om dit te doen heeft men een meer verfijnde vorm van de rest van de trapeziumformule nodig. In de onderstelling dat f voldoende afleidbaar is, kan men bewijzen dat er constanten bestaan zodat Z x1 x0

h f (x)dx = ( f0 + f1 ) + c1 h3 + c2 h5 + · · · + cN h2N+1 f (2N) (ξ) 2

(9.20)

Wat we hierboven gedaan hebben is het elimineren van de constante c3 . De fout wordt dan van de orde h5 , en dit is inderdaad wat we voor de formule van Simpson gevonden hebben. De volgende stap bestaat erin de term in h5 te elimineren, zodat de fout van de orde h7 wordt. Laten we dit expliciet uitvoeren. Als we in (9.20) alle termen vanaf c3 h7 verwaarlozen, dan vinden we I ≈ A1,1 + c02 (2h)5 Als we het interval in twee delen, en de formule op elk interval toepassen, dan krijgen we I ≈ A2,1 + 2c02 h5 We elimineren c02 , en stellen het resultaat gelijk aan A2,2 . Dit levert A2,2 =

1 (16A2,1 − A1,1 ) 15

Meer algemeen krijgen we Am,2 =

1 (16Am,1 − Am−1,1 ) 15 105

We kunnen deze redenering herhalen, en vinden de volgende algemene formule Am,n+1 =

4n+1 Am,n − Am−1,n 4n+1 − 1

(9.21)

Dit leidt tot de volgende tabel A0,0 A1,0

A1,1

A2,0

A2,1

A2,2

A3,0 .. .

A3,1 .. .

A3,2 .. .

(9.22) A3,3 .. .

We verwachten dat de kolommen sneller en sneller naar de gezochte integraal I zullen convergeren. We hebben inderdaad het volgende resultaat (zonder bewijs). Stelling 9.3.1 Als f analytisch over het interval [a, b], dan convergeren de kolommen van (9.22) allemaal naar I. Bovendien convergeert elke kolom sneller dan de vorige, dit wil zeggen dat I − Am,n+1 lim =0 m→∞ I − Am,n Opmerkingen 9.3.2 1) De methode heeft slechts zin als de functiewaarden in de spilpunten voldoende nauwkeurig bekend zijn. Zoniet maken de fouten op de fi alle verbeteringen zinloos. 2) Bij overgang van Am−1,0 naar Am,0 hoeft men f (x) alleen te berekenen in de extra toegevoegde spilpunten. Voorbeeld 9.3.3 We passen de methode van Romberg toe op Z 1 dx

= ln 2 = 0.69314 71805 59945 . . . 1+x Rombergintegratie levert ons de volgende tabel 0

.75000 00000 00 .70833 33333 34

.69444 44444 45

.69702 38095 25

.69325 39682 53

.69317 46031 73

.69412 18503 72

.69315 45306 53

.69314 79014 73

.69314 74776 37

.69339 12022 06

.69314 76528 17

.69314 71942 93

.69314 71830 68

.69314 71819 1

.69320 82082 69

.69314 72102 9

.69314 71807 87

.69314 71805 73

.69314 71805 65

9.4

.69314 71805 64

De integratieformules van Gauss

Tot nu toe kozen we onze spilpunten ekwidistant. Men kan ze ook willekeurig verspreid over [a, b] nemen. We kiezen n spilpunten x1 , x2 , . . . , xn , en we integreren de formule van Lagrange: Z b a

n

f (x)dx = ∑ ci fi + E( f ) i=1

106

Merk op dat we onze spilpunten anders nummeren dan in de voorgaande paragrafen. De integratieformule is zeker juist voor f een veelterm van graad ten hoogste n − 1. Het idee van Gauss bestaat er nu in om de spilpunten xi zo te kiezen dat E(xn ) = E(xn+1 ) = · · · = E(x2n−1 ) = 0 Het getal k is dan 2n−1, zodat we een zeer nauwkeurige integratieformule hebben. Om de spilpunten x1 , x2 , . . . , xn te bepalen moeten we in feite een stelsel van n vergelijkingen met n onbekenden oplossen. We gaan echter anders te werk. Vooreerst voeren we de lineaire substitutie x=

a+b b−a + t 2 2

uit, zodat het integratieinterval [−1, 1] wordt. Beschouw nu de veeltermen van Legendre Pn (x). Deze werden bestudeerd in volumes 3 en 4 van de cursussen “Analyse¨en “Differentiaal- en Integraalrekening”. We herhalen de volgende eigenschappen: 1) Pn is een veelterm van graad n; 2) Pn heeft n enkelvoudige nulpunten x1 , x2 , . . . , xn die allemaal in het interval [−1, 1] gelegen zijn; 3) Als f een veelterm van graad strikt kleiner dan n is, dan is Z 1 −1

f (x)Pn (x)dx = 0

We kiezen nu voor de spilpunten in onze integratieformule de n nulpunten x1 , x2 , . . . , xn van Pn (x). Zij f een veelterm van graad ten hoogste 2n − 1. We delen f (x) door Pn (x), en we krijgen f (x) = q(x)Pn (x) + r(x) waarbij q(x) en r(x) veeltermen zijn van graad kleiner dan of gelijk aan n − 1. Daarom is E(r(x)) = 0 en Z 1

E(q(x)Pn (x)) =

−1

n

q(x)Pn (x)dx − ∑ ci q(xi )Pn (xi ) i=1

= 0−0 = 0 We concluderen dat E( f ) = 0, zodat we het idee van Gauss verwezenlijkt hebben. We zullen de methode van Gauss nog veralgemenen, maar eerst hebben we een interpolatieformule nodig. De interpolatieformule van Hermite Beschouw n spilpunten x1 , x2 , . . . , xn . Gevraagd wordt om een veelterm P(x) te vinden zodanig dat P(xi ) = f (xi ) = fi P0 (xi ) = f 0 (xi ) = fi0 107

(9.23) (9.24)

voor i = 1, . . . , n, waarbij f een gegeven functie is. We gaan als volgt tewerk: ons inspirerend op de interpolatieveelterm van Lagrange proberen we een veelterm P van de vorm n

n

i=1

i=1

P(x) = ∑ fi hi (x) + ∑ fi0 hi (x) waarbij hi (x) en hi (x) veeltermen van graad 2n − 1 zijn. Als we ervoor kunnen zorgen dat hi (x j ) = δi j h0i (x j ) = 0

hi (x j ) = 0

(9.25)

0 hi (x j ) = δi j

(9.26)

dan is aan (9.23) en (9.24) voldaan en ons probleem is opgelost. We nemen de elementaire veelterm van Lagrange n

Li (x) =

x−xj j=1, j6=i xi − x j

∏

De veelterm Li2 (x) is dan van graad 2n − 2. Bovendien is Li2 (x j ) = δi j d 2 (L )(x j ) = 0 dx i

(9.27) (9.28)

voor i 6= j. We stellen hi (x) = ri (x)Li2 (x) en hi (x) = si (x)Li2 (x) waarbij r(x) en s(x) veeltermen van graad 1 zijn. Aan (9.25) en (9.26) wordt dan voldaan voor i 6= j. We bepalen tenslotte ri (x) en si (x) door te eisen dat aan (9.25) en (9.26) wordt voldaan voor i = j. Dit geeft ri (xi ) = 1 si (xi ) = 0

ri0 (xi ) + 2Li0 (xi ) = 0 s0i (xi ) = 1

zodat ri (x) = 1 − 2Li0 (xi )(x − xi ) en si (x) = x − xi

(9.29)

en we hebben de veelterm P(x) gevonden. We noemen P(x) de interpolatieveelterm van Hermite We willen vervolgens een formule vinden voor de fout r(x) = f (x) − P(x). We onderstellen dat f (2n) eindig is over een interval I dat de spilpunten x1 , x2 , . . . , xn bevat, en beschouwen de hulpfunctie n

F(t) = f (t) − P(t) − k ∏(t − xi )2 i=1

We zien onmiddellijk dat x1 , x2 , . . . , xn nulpunten zijn van zowel F(t) als F 0 (t). Kies x ∈ I vast, en verschillend van elk van de spilpunten x1 , x2 , . . . , xn . We stellen k=

f (x) − P(x) ∏ni=1 (x − xi )2 108

zodat ook x een nulpunt is van F(t). Gebruik makend van de stelling van Rolle vinden we dat F 0 (t) n nulpunten x10 , x20 , . . . , xn0 heeft die tussen de punten x1 , x2 , . . . , xn en x liggen. Deze nieuwe nulpunten zijn verschillend van x1 , x2 , . . . , xn , en dus heeft F 0 2n verschillende nulpunten x1 , x2 , . . . , xn , x10 , x20 , . . . , xn0 . Herhaaldelijk gebruik makend van de stelling van Rolle vinden we dat F 00 2n − 1 nulpunten heeft, F (3) 2n − 2, enzovoort. F (2n) heeft tenminste 1 nulpunt ξ ∈ I. F (2n) (ξ) = f (2n) (ξ) − k(2n)! = 0 We concluderen dat

f (2n) (ξ) (2n)!

k= en f (x) = P(x) +

f (2n) (ξ(x)) n (x − xi )2 (2n)! ∏ i=1

We vatten onze resultaten samen in de volgende stelling. Stelling 9.4.1 Onderstel dat f : I → R een eindige afgeleide van orde 2n heeft, en neem spilpunten x1 , x2 , . . . , xn ∈ I. Dan is f (x) = P(x) + r(x) waarbij n

P(x) =

n

∑ fihi(x) + ∑ fi0hi(x)

i=1

hi (x) =

i=1

0


1 − 2L (xi )(x − xi ) Li2 (x)

hi (x) = (x − xi )Li2 (x) r(x) =

f (2n) (ξ(x)) n (x − xi )2 (2n)! ∏ i=1

De integratieformule van Gauss We onderstellen dat f : I = [a, b] → R voldoet aan de voorwaarden van stelling 9.4.1. Zij w : I = [a, b] → R een gewichtsfunctie (m.a.w. w(x) ≥ 0). We vermenigvuldigen de interpolatieformule van Hermite met w(x) en integreren van a tot b. Z b a

n

n

i=1

i=1

w(x) f (x)dx = ∑ Hi fi + ∑ H i fi0 + E( f )

met Z b

Hi =

a

w(x)hi (x)dx

Z b

Hi =

a

w(x)hi (x)dx

1 E( f ) = (2n)!

Z b a

n

f (2n) (ξ(x))w(x) ∏(x − xi )2 dx i=1

109

(9.30)

(9.30) is een integratieformule die juist is voor veeltermen tot en met graad 2n − 1. Het nadeel van de formule is dat men de afgeleiden van f in de spilpunten moet kennen. De truc bestaat er nu in om de spilpunten zo te kiezen dat de co¨effici¨enten H i wegvallen. Dit kunnen we doen door gebruik te maken van orthogonale veeltermen. Zij pn (x) een rij orthogonale veeltermen over [a, b] ten opzichte van de gewichtsfunctie w(x). Uit de resultaten van volume 4 van de cursussen “Analyse¨en “Differentiaal- en integraalrekening”weten we dat zulk een rij veeltermen bestaat. Bovendien heeft pn (x) n enkelvoudige nulpunten die allen in het interval [a, b] liggen. We kiezen als spilpunten nu de nulpunten van pn (x). Merk op dat pn (x) = A(x − x1 ) · · · (x − xi ) · · · (x − xn ) = ai (x − xi )Li (x) zodat

1 b w(x)pn (x)Li (x)dx = 0 ai a a aangezien gr(Li (x)) = n − 1 < n. We vinden dus een integratieformule Z b

Hi =

Z

w(x)(x − xi )Li2 (x)dx =

Z b a

n

w(x) f (x)dx = ∑ Hi fi + E( f )

(9.31)

i=1

die juist is voor veeltermen tot en met graad 2n − 1. Voor de co¨effici¨enten Hi vinden we Z b

Hi =

w(x)hi (x)dx

a

Z b

= a

Z b

= a

w(x)Li2 (x)dx − 2Li0 (xi )

Z b a

w(x)(x − xi )Li2 (x)dx

w(x)Li2 (x)dx

(9.32)

We kunnen voor Hi een meer eenvoudige formule opstellen. Neem f (x) = L j (x) in (9.31). E(L j ) = 0, omdat L j van graad n − 1 ≤ 2n − 1 is, en Z b a

n

w(x)L j (x)dx = ∑ Hi L j (xi ) = H j i=1

en we vinden

Z b

Hi =

a

w(x)Li (x)dx

(9.33)

De integratieformule (9.31) is dus de formule die men verkrijgt door te vertrekken van de interpolatieformule van Lagrange met als spilpunten de nulpunten van pn (x), en die te integreren na vermenigvuldiging met w(x). We zullen nu een formule voor de rest opstellen. Onderstel dat f (2n) (ξ(x)) een continue functie van x is. Met behulp van de stelling van het gemiddelde vinden we 1 E( f ) = (2n)! =

Z b a

n

f (2n) (ξ(x))w(x) ∏(x − xi )2 dx i=1

Z f (2n) (ξ) b

A2 (2n)!

a

w(x)pn (x)2 dx

De rest is dus van de vorm k f (2n) (ξ). 110

De formules van Legendre-Gauss Als het integratieinterval eindig is, dan kan men het door een lineaire substitutie herleiden tot [−1, 1]. Kiest men w(x) = 1, dan vindt men als orthogonale veeltermen de veeltermen van Legendre Pn . 1 dn (1 − x2 )n Pn (x) = (−1)n 2n dxn De integratieformule neemt de volgende vorm aan Z 1 −1

n

f (x)dx = ∑ Hi f (xi ) + i=1

22n+1 (n!)4 f (2n) (ξ) 3 (2n + 1)((2n)!)

(9.34)

In de hiernavolgende tabel geven we de spilpunten xi , en de co¨effici¨enten Hi . Spilpunten xi

Co¨effici¨enten Hi

n=2

±0.577 350 27

1

n=3

0

0.888 888 89

±0.774 596 67

0.555 555 56

±0.339 981 04

0.652 145 15

±0.861 136 31

0.347 854 85

0

0.568 888 89

±0.538 469 31

0.478 628 67

±0.906 179 85

0.236 926 89

n=4 n=5

De formules van Laguerre-Gauss We werken over het interval [0, +∞), met gewichtsfunctie w(x) = e−x . De orthogonale veeltermen die we aldus verkrijgen zijn de veeltermen van Laguerre Ln (x) = ex

dn n −x (x e ) dxn

meer bepaald L0 (x) L1 (x) L2 (x) L3 (x) L4 (x)

= = = = =

1 1−x 2 − 4x + x2 6 − 18x + 9x2 − x3 24 − 96x + 72x2 − 16x3 + x4

111

√ Voor√n = 2 zijn de spilpunten de nulpunten 2 ± 2 van L2 (x). De bijhorende gewichten zijn (2 ∓ 2)/4, en we krijgen de volgende integratieformule. Z +∞

e−x f (x)dx =

0

√ √ √  f (4) (ξ) √ 1 (2 + 2) f (2 − 2) + (2 − 2) f (2 + 2) + 4 6

Na een lineaire transformatie vinden we een formule die toelaat om de Laplacegetranformeerde te benaderen. √ √ Z +∞ √ √ 1 2− 2 2 + 2  f (4) (ξ) −xp L ( f )(p) = e f (x)dx = (2 + 2) f ( ) + (2 − 2) f ( ) + 4p p p 6p5 0 De algemene integratieformule luidt als volgt Z +∞ 0

n

e−x f (x)dx = ∑ Hi f (xi ) + i=1

(n!)2 (2n) (ξ) f (2n)!

(9.35)

In de hiernavolgende tabel geven we de spilpunten xi , en de co¨effici¨enten Hi .

n=2 n=3

n=4

n=5

Spilpunten xi

Co¨effici¨enten Hi

0.585 786 44

0.853 553 39

3.414 213 56

0.146 446 61

0.415 774 56

0.711 093 01

2.294 280 36

0.278 517 73

6.289 945 08

0.010 389 26

0.322 547 69

0.603 154 10

1.745 761 10

0.357 418 69

4.536 620 30

0.038 887 91

9.395 070 91

0.000 539 29

0.263 560 32

0.521 755 61

1.413 403 06

0.398 666 81

3.596 425 77

0.075 942 45

7.085 810 01

0.003 611 76

12.640 800 84

0.000 023 37

De formules van Hermite-Gauss 2

We werken over het interval (−∞, +∞), met gewichtsfunctie w(x) = e−x . De orthogonale veeltermen die we aldus verkrijgen zijn de veeltermen van Hermite. Hn (x) = (−1)n ex 112

2

dn −x2 (e ) dxn

De algemene integratieformule luidt als volgt Z +∞

−x2

e −∞

√ n! π (2n) f (x)dx = ∑ Hi f (xi ) + n f (ξ) 2 (2n)! i=1 n

(9.36)

In de hiernavolgende tabel geven we de spilpunten xi , en de co¨effici¨enten Hi . Spilpunten xi

Co¨effici¨enten Hi

n=2

±0.707 106 78

0.886 226 93

n=3

0

1.181 635 90

±1.224 744 87

0.295 408 98

±0.524 647 62

0.804 914 09

±1.650 680 12

0.081 312 84

0

0.945 308 72

±0.958 572 46

0.393 619 32

±2.020 182 87

0.019 953 24

n=4 n=5

Voor een overzicht van integratieformules opgesteld met behulp van nog andere orthogonale veeltermen, en voor meer uitgebreide tabellen verwijzen we naar M. Abramovitz en I. Stegun, Handbook of Mathematical functions, Dover Publ. Inc., New York, Hoofdstuk 25.

113

Hoofdstuk 10 Differentiaalvergelijkingen In de eerste vier paragrafen bespreken we differentiaalvergelijkingen van orde 1 y0 = f (x, y)

(10.1)

y(x0 ) = y0

(10.2)

met beginvoorwaarde De algoritmes die we zullen ontwikkelen zijn ook van toepassing op stelsels differentiaalvergelijkingen d~y = f (x,~y) dx met beginvoorwaarde ~y(x0 ) =~y0 . Doordat een differentiaalvergelijking van orde n steeds te herleiden is tot een stelsel differentiaalvergelijkingen van n vergelijkingen en n onbekenden van orde 1, vinden we in principe ook algoritmen voor differentiaalvergelijkingen van hogere orde met beginvoorwaarden. Een numeriek veel ingewikkelder probleem is het oplossen van een differentiaalvergelijking van orde 2 met randvoorwaarden. Een korte bespreking geven we in paragraaf 5. In een laatste paragraaf worden parti¨ele differentiaalvergelijkingen zeer kort behandeld.

10.1

De methode van Euler-Heun

We zullen (10.1) integreren met behulp van de methode van de discretisatie. Dit wil zeggen dat we de onbekende functie y zullen trachten te berekenen in de spilpunten xn = x0 + nh

n = 0, 1, 2, . . .

Het getal h noemen we de pas. De benaderde waarden voor y(xn ) noteren we yn . (10.1) en (10.2) zijn equivalent met Z x

y(x) = y0 +

f (x, y)dx x0

114

(10.3)

We benaderen de integraal in het rechterlid met een van de methodes uit het voorgaande hoofdstuk. Stel x = x1 , en bereken de integraal in het rechterlid met de trapeziumregel. Dit levert h y1 ≈ y0 + ( f (x0 , y0 ) + f (x1 , y1 )) 2

(10.4)

(10.4) is een vergelijking in y1 , en is van de vorm y = g(y). We lossen deze op via iteratie. Een eerste benadering voor y1 is (10.5) y∗1 = y0 + h f (x0 , y0 ) Dit betekent dat we de kromme y = y(x) vervangen door de raaklijn in het punt (x0 , y0 ). (10.5) noemt men de voorspellingsformule Men vervangt y∗1 in het rechterlid van (10.4) en voert een y1* y1

y0 x0 0

x1

1

Figuur 10.1: De methode van Euler-Heun eerste iteratie uit

h (10.6) y1 = y0 + ( f (x0 , y0 ) + f (x1 , y∗1 )) 2 (10.6) noemt men de verbeteringsformule. Meestal stopt men na e´ e´ n iteratie. Het heeft geen zin om de wortels van (10.4) nauwkeuriger te bepalen, aangezien (10.4) zelf reeds een benadering is, omdat men de rest in de trapeziumformule niet kan meenemen. Als we ons tevreden stellen met de voorspellingsformule (10.5), dan noemen we onze methode de methode van Euler. (10.5) gecombineerd met (10.6) noemen we de methode van Heun. Nadat y1 berekend is, gaan we op dezelfde manier tewerk om achtereenvolgens y2 , y3 , . . . te berekenen. Merk op dat de afrondingsfouten en methodefouten groter en groter worden naarmate we verder en verder van x0 weggaan. Als men een kleinere pas h neemt, dan worden de methodefouten kleiner. Men heeft dan evenwel meer iteraties nodig om bij een gegeven waarde van x te geraken, zodat de afrondingsfouten dan weer groter worden. In de praktijk heeft het dus niet altijd zin om h te klein te nemen.

10.2

De methode van Runge-Kutta

Er bestaan een groot aantal meer ingewikkelde varianten van de methode van Euler-Heun. Een veel gebruikte variante is die van Runge-Kutta. Hierbij gebruikt men de methode van Simpson in plaats van de trapeziumregel. Stel x = x2 in (10.3), en benader de integraal met de regel van Simpson:
h y2 ≈ y0 + f (x0 , y0 ) + 4 f (x1 , y1 ) + f (x2 , y2 ) (10.7) 3 115

y0 x0

y2

y1 x1

x2

Figuur 10.2: De methode van Runge-Kutta In het rechterlid staan er nu twee grootheden die we niet kennen: y1 en y2 . We gaan als volgt tewerk. We benaderen y1 met behulp van de methode van Heun, via (10.5) en (10.6). Als we een benadering y1 voor y(x1 ) hebben, dan hebben we ook een benadering voor y0 (x1 ), namelijk y01 = f (x1 , y1 ). We hebben dus y∗1 = y0 + h f (x0 , y0 ) h y1 = y0 + ( f (x0 , y0 ) + f (x1 , y∗1 )) 2 0 y1 = f (x1 , y1 ) We schrijven dan de veelterm p3 van Hermite neer voor de spilpunten x0 en x1 . p3 is van graad 3. We benaderen dan y2 in het rechterlid van (10.7) door p3 (x2 ). Gebruik makend van (9.29) vinden we  (x − x )2  2 1 h0 (x) = 1 + (x − x0 ) h h2  (x − x )2  2 0 h1 (x) = 1 − (x − x1 ) 2 h h (x − x1 )2 h0 (x) = (x − x0 ) h2 (x − x0 )2 h1 (x) = (x − x1 ) h2 en h1 (x2 ) = −4 h1 (x2 ) = 4h

h0 (x2 ) = 5 h0 (x2 ) = 2h Als voorspellingsformule voor y2 krijgen we dus

y∗2 = p3 (x2 ) = 5y0 − 4y1 + 2hy00 + 4hy01

(10.8)

De verbeteringsformule is (10.7). We vatten alles samen als volgt. We schrijven fi = f (xi , yi ) en fi∗ = f (xi , y∗i ). y∗1 = y0 + h f0 116

h y1 = y0 + ( f0 + f1∗ ) 2 y∗2 = 5y0 − 4y1 + 2h f0 + 4h f1 = y0 − 2h f1∗ + 4h f1 h y2 = y0 + ( f0 + 4 f1 + f2∗ ) 3

(10.9)

In deze formules worden dikwijls kleine getallen bij grotere opgeteld. Het is numeriek stabieler om eerst de kleine getallen bij elkaar op te tellen, en daarom herschrijven we ons algoritme als volgt. 1 y2 = y0 + (k1 + 4k3 + k4 ) (10.10) 6 met k1 = 2h f0 = 2h f (x0 , y0 ) k1 k2 + ) 4 4 k 1 = 2h f1∗ = 2h f (x0 + h, y0 + ) 2 = 2h f2∗ = 2h f (x0 + 2h, y0 − k2 + 2k3 )

k3 = 2h f1 = 2h f (x0 + h, y0 + k2 k4

Hiermee is tevens het spilpunt x1 uit ons algoritme verdwenen. We hernummeren daarom de spilpunten: x2 wordt vervangen door x1 , en h door h/2. Het oorspronkelijk interval [x0 , x2 ] (na hernummering [x0 , x1 ]) krijgt dus lengte h. We verkrijgen nu het algoritme Runge-Kutta I. k1 = h f (x0 , y0 ) k1
h k2 = h f x0 + , y0 + 2 2 h k1 k2
k3 = h f x0 + , y0 + + 2 4 4
k4 = h f x0 + h, y0 − k2 + 2k3 1 y1 = y0 + (k1 + 4k3 + k4 ) 6

(10.11)

Een lichte variante van dit algoritme is het algoritme van Runge-Kutta II. Dit is iets eenvoudiger, en levert in de praktijk even goede resultaten op. k1 = h f (x0 , y0 ) h k1
k2 = h f x0 + , y0 + 2 2 h k2
k3 = h f x0 + , y0 + 2 2
k4 = h f x0 + h, y0 + k3 1 y1 = y0 + (k1 + 2k2 + 2k3 + k4 ) 6 117

(10.12)

10.3

De orde van een eenstapsmethode

De methoden uit de vorige twee paragrafen zijn van de vorm y1 = y0 + hφ(x0 , y0 , h)

(10.13)

We spreken van een eenstapsmethode. Onderstel nu dat z de exacte oplossing is van de differentiaalvergelijking (10.1) met beginvoorwaarde (10.2). Het differentiequoti¨ent van de exacte oplossing wordt gegeven door de formule   z(x0 + h) − y0 als h 6= 0 (10.14) ∆(x0 , y0 , h) = h  f (x0 , y0 ) als h = 0 Voor de benaderde oplossing is het differentiequoti¨ent y1 − y0 = φ(x0 , y0 , h) h Hoe kleiner het verschil tussen deze twee differentiequoti¨enten, hoe beter y1 de werkelijke functiewaarde z(x1 ) benadert. Daarom stellen we τ(x0 , y0 , h) = ∆(x0 , y0 , h) − φ(x0 , y0 , h)

(10.15)

Voor een “redelijk¨algoritme geldt steeds dat lim τ(x0 , y0 , h) = 0

h→0

Als we aannemen dat ∆ continu is in h = 0 volgt dus: lim φ(x0 , y0 , h) = f (x0 , y0 )

h→0

Als dit het geval is, dan noemen we onze methode consistent. We noemen onze methode van orde p indien τ(x0 , y0 , h) = O(h p ) Dit betekent dat

τ(x0 , y0 , h) ≤M hp

voor h voldoende klein. Hieraan is voldaan indien de eerste niet triviale term in de Taylorontwikkeling van τ(x0 , y0 , h) als functie van h in h p staat. We ontwikkelen z in Taylor: z(x0 + h) = z(x0 ) + hz0 (x0 ) + en

h2 00 z (x0 ) + · · · 2

h ∆(x0 , y0 , h) = z0 (x0 ) + z00 (x0 ) + · · · 2 118

Nu is z0 (x0 ) = f (x0 , y0 ) z0 (t) = f (t, z(t)) ∂f ∂f z00 (t) = (t, z(t)) + (t, z(t))z0 (t) ∂x ∂y ∂f ∂f z00 (x0 ) = (x0 , y0 ) + (x0 , y0 ) f (x0 , y0 ) ∂x ∂y zodat

h ∆(x0 , y0 , h) = f (x0 , y0 ) + 2



 ∂f ∂f (x0 , y0 ) + (x0 , y0 ) f (x0 , y0 ) + · · · ∂x ∂y

(10.16)

Voorbeeld 10.3.1 De methode van Euler In dit geval is φ(x0 , y0 , h) = f (x0 , y0 ), en   ∂f h ∂f (x0 , y0 ) + (x0 , y0 ) f (x0 , y0 ) τ(x0 , y0 , h) = 2 ∂x ∂y De methode van Euler is dus van orde 1. Voorbeeld 10.3.2 De methode van Euler-Heun Nu is
 1 φ(x0 , y0 , h) = f (x0 , y0 ) + f x0 + h, y0 + h f (x0 , y0 ) 2   ∂f h ∂f (x0 , y0 ) + (x0 , y0 ) f (x0 , y0 ) + · · · = f (x0 , y0 ) + 2 ∂x ∂y Tot op orde 1 vallen de Taylor ontwikkelingen van φ en ∆ samen. Hieruit volgt dat de methode van Euler-Heun van orde twee is. Voorbeeld 10.3.3 De gemodifieerde methode van Euler Hierbij is
h h φ(x0 , y0 , h) = f x0 + , y0 + f (x0 , y0 ) 2 2 en dus
h h y1 = y0 + h f x0 + , y0 + f (x0 , y0 ) 2 2 Toon zelf aan dat de gemodifieerde methode van Euler van orde 2 is.

(10.17)

(10.18)

Voorbeeld 10.3.4 De methode van Runge-Kutta Een berekening analoog aan die van hierboven (maar ingewikkelder) toont aan dat de algoritmes van Runge-Kutta (I en II) van orde 4 zijn.

119

10.4

Meerstapsmethodes

Bij de voorgaande methodes werd enkel yn gebruikt om yn+1 te berekenen. Gebruikt men yn−k , yn−k+1 , . . . , yn , dan heeft men een meerstapsmethode. Een probleem is natuurlijk dat men in het begin slechts 1 waarde van y kent, namelijk y0 . Men lost dit op door y1 , . . . , yk te berekenen met een eenstapsmethode, bijvoorbeeld een algoritme van Runge-Kutta. Men vertrekt van Z xn+1 yn+1 = yn + f (x, y)dy xn

Om de integraal te berekenen vervangt men f (x, y) door een interpolatieveelterm waarbij men xn−k , xn−k+1 , . . . , xn en eventueel xn+1 als spilpunten gebruikt. Als men enkel xn−k , xn−k+1 , . . . , xn gebruikt, dan bekomt men een expliciete formule (de AdamsBashfort formule). Voorbeeld 10.4.1 We interpoleren in de punten xn−1 en xn . We schrijven fn = f (xn , yn ). De interpolatieveelterm is x − xn−1 x − xn p(x) = fn−1 + fn xn−1 − xn xn − xn−1 en we vinden fn−1 xn+1 fn = yn − (x − xn )dx + h xn h h = yn + (3 fn − fn−1 ) 2 Z

yn+1

Z xn+1 xn

(x − xn−1 )dx (10.19)

Op analoge wijze vinden we, gebruik makend van de spilpunten xn−2 , xn−1 en xn : yn+1 = yn +

h (23 fn − 16 fn−1 + 5 fn−2 ) 12

(10.20)

en, gebruik makend van de spilpunten xn−3 , xn−2 , xn−1 en xn : yn+1 = yn +

h (55 fn − 59 fn−1 + 37 fn−2 − 9 fn−3 ) 24

(10.21)

Als men naast xn−k , xn−k+1 , . . . , xn ook xn+1 als spilpunt gebruikt, dan bekomt men een impliciete formule (de Adams-Moulton formule). Gebruik makend van spilpunten xn−1 , xn en xn+1 vinden we yn+1 = yn +

h (5 fn+1 + 8 fn − fn−1 ) 12

(10.22)

Met spilpunten xn−2 , xn−1 , xn en xn+1 vinden we yn+1 = yn +

h (9 fn+1 + 19 fn − 5 fn−1 + fn−2 ) 24 120

(10.23)

De impliciete formules kunnen als verbeteringsformules gebruikt worden. Deze moeten dan gecombineerd worden met een expliciete formule als voorspellingsformule. Er zijn allerlei combinaties mogelijk. Een veelgebruikt algoritmes is het volgende. 1) Bereken y1 , y2 , y3 met Runge-Kutta; 2) gebruik (10.21) als voorspellingsformule; 3) gebruik (10.23) als verbeteringsformule.

10.5

Differentiaalvergelijkingen met randvoorwaarden

We bekijken nu de differentiaalvergelijking y00 = f (x, y, y0 )

(10.24)

y(a) = y0 en y(b) = y1

(10.25)

met randvoorwaarden

Eerste methode: de shooting method We trachten het vraagstuk op te lossen door y0 (a) te raden: we nemen als beginvoorwaarden y(a) = y0 en y0 (a) = λ en integreren de vergelijking met een van de methodes uit de vorige paragrafen. De verkregen oplossing noemen we yλ . Meestal is yλ (b) 6= y1 . We krijgen een algebraische vergelijking

•

y0

.

•

y1

a

b

x

Figuur 10.3: De shooting method yλ = y1 met als onbekende de parameter λ. We kunnen deze oplossen met behulp van de methode van de koorde. Het nadeel is wel dat we bij elke iteratie een numerieke integratie moeten uitvoeren. Voorbeeld 10.5.1 We bekijken de lineaire differentiaalvergelijkingen y00 + p(x)y0 + q(x)y = r(x) 121

met randvoorwaarden (10.25). We bepalen eerst de oplossing yI (x) die voldoet aan de beginvoorwaarden y(a) = y0 en y0 (a) = 1 Daarna bepalen we de oplossing yII (x) van de homogene vergelijking y00 + p(x)y0 + q(x)y = 0 met beginvoorwaarden y(a) = 0 en y0 (a) = 1 De functie y(x) = yI (x) + αyII (x) is dan een oplossing van de volledige vergelijking die voldoet aan de eerste randvoorwaarde y(a) = y0 . Het volstaat nu om α zodanig te kiezen dat y(b) = y1 (b) + αy2 (b) = y1 Discretisatie Men stelt h = (b − a)/n en xk = a + kh. Het interval [a, b] wordt dus in n gelijke delen verdeeld. De bedoeling is om de oplossing y(x) in de spilpunten te berekenen. In de differentiaalvergelijking (10.24) vervangt men de afgeleiden door een formule voor numerieke afleiding: y(k) wordt vervangen door een uitdrukking van de vorm n

∑ αiyi

i=0

Dit levert in elk spilpunt een (lineaire) vergelijking met y0 , y1 , . . . , yn als onbekenden. Voorbeeld 10.5.2 We bekijken het Sturm-Liouville probleem y00 + λy = 0 met randvoorwaarden y(0) = y(π) = 0 De eigenwaarden van het Sturm-Liouville probleem zijn die waarden van λ waarvoor er van nul verschillende oplossingen zijn (zie Volume 4 van de cursus “Analyse”). In dit geval zijn de eigenwaarden λ = m2 , met m = 1, 2, 3, . . .. De bijhorende eigenfuncties zijn ym = sin mx We proberen deze eigenfuncties numeriek terug te vinden. We delen het interval (0, π) op in n gelijke delen: xk = kh, met h = π/n. We hebben een numerieke formule voor de tweede afgeleide nodig. Uit de formule van Taylor volgt h2 00 h3 000 h4 (4) y + y + y (ξ1 ) 2 k 6 k 24 h2 h3 h4 (4) = yk − hy0k + y00k − y000 + y (ξ2 ) 2 6 k 24

yk+1 = yk + hy0k + yk−1

122

Optelling geeft y00 (xk ) =

yk+1 − 2yk + yk−1 + O(h2 ) h2

We verkrijgen het lineair stelsel   y0 = 0 y + (λh2 − 2)yk + yk−1 = 0 k = 1, 2 . . . , n − 1  k−1 yn = 0 Dit is een lineair homogeen stelsel. Om niet-triviale oplossingen te hebben moet de determinant nul zijn. Dit levert een vergelijking in λ, waaruit een benadering van de eigenwaarde kan gevonden worden. Met die waarden van λ kunnen we dan het lineaire stelsel oplossen, waaruit de eigenfuncties volgen. Voor n = 2 wordt het stelsel   y0 = 0 y + (λh2 − 2)y1 + y2 = 0  0 y2 = 0 Er is een niet-triviale oplossing als λ=

2 8 = 2 ≈ 0.81 2 h π

We vinden dus in benadering de eerste eigenwaarde λ = 1 terug. Voor n = 3 wordt het stelsel, na vereenvoudiging  (λh2 − 2)y1 + y2 = 0 y1 + (λh2 − 2)y2 = 0 Als we de determinant van het stelsel gelijk aan nul stellen vinden we (λh2 − 2)2 = 1 Een van de wortels is λ1 =

9 ≈ 0.91 π2

λ2 =

27 ≈ 2.73 π2

De andere wortel is

Collocatiemethodes We bekijken de lineaire differentiaalvergelijking y00 + a(x)y0 + b(x)y = f (x)

(10.26)

y(x0 ) = y0 en y(x1 ) = y1

(10.27)

met randvoorwaarden 123

Het basisidee bestaat erin om een stel functies ϕ0 (x), . . . , ϕn (x) te kiezen zodat ϕ0 (x0 ) = y0 en ϕ0 (x1 ) = y1 en

ϕi (x0 ) = ϕi (x1 ) = 0

Een lineaire combinatie van de vorm n

ϕ(x) = ϕ0 (x) + ∑ αk ϕk (x) k=1

voldoet dan steeds aan de randvoorwaarden (10.25). Men tracht dan α1 , . . . , αn zo te bepalen dat ϕ(x) zo goed mogelijk voldoet aan de differentiaalvergelijking: ϕ00 +a(x)ϕ0 +b(x)ϕ moet “zo dicht mogelijk bij” f liggen. Dit kan op verschillende manieren gebeuren. Methode der spilpunten Men kiest n spilpunten ξ1 , ξ2 , . . . , ξn , en men bepaalt de αi zodanig dat ϕ00 (ξk ) + a(ξk )ϕ0 (ξk ) + b(ξk )ϕ(ξk ) = f (ξk ) in elk van de spilpunten. Dit geeft een lineair stelsel van n vergelijkingen met n onbekenden α1 , . . . , αn . Methode van de functieruimte Men bepaalt de co¨effici¨enten α1 , . . . , αn zodanig dat de afstand tussen de functies f en ϕ00 + a(x)ϕ0 + b(x)ϕ minimaal is. Hierbij heeft men verschillende mogelijke keuzes voor de afstand tussen twee functies, bijvoorbeeld rZ x 1 d( f , g) = ( f (x) − g(x))2 dx x0

10.6

Parti¨ele differentiaalvergelijkingen

Numeriek integreren van parti¨ele differentiaalvergelijkingen is een vak op zichzelf! We beperken ons hier tot een enkel eenvoudig voorbeeld. Gegeven is een rechthoek ABCD waarvan de zijden evenwijdig zijn met de co¨ordinaatassen, en een functie f (x, y) gedefinieerd op deze rechthoek. Gevraagd wordt om een functie z(x, y) te bepalen zodat ∂2 z ∂2 z + = f (x, y) ∂x2 ∂y2 binnen de rechthoek, met randvoorwaarden z=0 ∂z =0 ∂x

op de intervallen [A, B] en [C, D] op de intervallen [A, D] en [C, B] 124

D

C 1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

k

A

B

h

Figuur 10.4: Parti¨ele differentiaalvergelijkingen We discretiseren onze rechthoek: we leggen een rooster over de beschouwde rechthoek, en zoeken z in de roosterpunten. Als we dit doen zoals op de tekening, dan vinden we 21 roosterpunten P1 , P2 , . . . , P21 waarin we de functiewaarden z1 , z2 , . . . , z21 trachten te benaderen. De roosterpunten op de intervallen [A, B] en [C, D] geven we geen nummer, want de functiewaarde is er toch nul. We vervangen de parti¨ele afgeleiden door een benaderende formule. In P10 hebben we bijvoorbeeld z9 − 2z10 + z11 ∂2 z ≈ ∂x2 h2 2 ∂ z z3 − 2z10 + z17 ≈ 2 ∂y k2 en de differentiaalvergelijking geeft ons in dit punt z9 − 2z10 + z11 z3 − 2z10 + z17 + = f10 h2 k2 We vinden een dergelijke vergelijking in elk punt. Dit geeft een lineair stelsel van 21 vergelijkingen met 21 onbekenden. Oplossen geeft een benadering in elk van de roosterpunten. Merk op dat de matrix van het lineair stelsel ijl is: in elk van de vergelijkingen treden slechts 5 co¨effici¨enten op die van nul verschillen. Merk ook op dat we in een punt op de rand moeten rekening houden met de randvoorwaarden. In P1 bijvoorbeeld, wordt de vergelijking ∂z −3z1 + 4z2 − z3 =0≈ ∂x 2h of −3z1 + 4z2 − z3 = 0 Als het gebied waarover we werken een willekeurige vorm heeft, dan is het uiteraard minder eenvoudig om met de randvoorwaarden rekening te houden. 125

Index Adams-Bashfort formule 120 Adams-Moulton formule 120 afrondingsfout 7 algoritme van Neville 76 Cauchy rij 20 conditiegetal 11 contractie 20 contractiestelling 20 dekpuntstelling 14 dekpuntstelling 20 discretisatie 114 drijvende komma representatie 8 eenstapsmethode 118 eigenvector 62 eigenwaarde 62 exponent 8 formule van Milne 101 formule van Newton 14 formule van Simpson 100 Forsythe veeltermen 89 gesloten Newton-Cotes formule 99 interpolatieveelterm van Hermite 108 interpolatieveelterm van Newton 80 kolomcriterium 34 machinegetallen 7 machinenauwkeurigheid 9 mantisse 8 methode van Aitken 27 methode van de vaste richting 15 methode van Euler 115 methode van Gauss-Seidel 45 methode van Heun 115 methode van Runge-Kutta 115 metrische ruimte 20 metrische ruimte 31 modelfouten 7 Moore-Penrose inversematrix 58

Morrey 30 Newton-Raphson 30 normaalvergelijkingen 88 numeriek stabiel algoritme 10 open Newton-Cotes formule 99 overflow 8 permutatiematrix 46 regula falsi 16 rijcriterium 34 Runge-Kutta I 117 Runge-Kutta II 117 samentrekking 20 simplexformule 103 singuliere waarden 56 totale rekenfout 7 trapeziumregel 99 tridiagonale matrix 53 underflow 8 veeltermen van Chebyshev 83 veeltermen van Hermite 112 veeltermen van Laguerre 111 veralgemeende inverse matrix 58 verbeteringsformule 115 voorspellingsformule 115 voorwaarde (V ) 34

126

Bibliografie [1] G.H. Golub, C.F. Van Loan, Matrix Computations, The John Hopkins Universtity Press, Baltimore, 1988. [2] I. Jacques, C. Judd, Numerical Analysis, Chapman and Hill, London 1987. [3] A. Ralston, A first course in Numerical Analysis, Mc Graw-Hill, New York, 1965. [4] J. Stoer, R. Bulirsch, Introduction to Numerical Analysis, Springer Verlag, Berlin, 1980.

127