Wiskunde: Voortgezette Analyse
S. Caenepeel
Syllabus 136 bij IR-WISK 11453 “Wiskunde: Voortgezette Analyse” Tweede Bachelor Ingenieurswetenschappen: Architectuur (SD-ID 004015)
2012
Inhoudsopgave I
Reeksen
4
1
Rijen en reeksen
5
1.1
Numerieke rijen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5
1.2
Numerieke reeksen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7
2
II 3
4
Reeksen van functies
9
2.1
Uniforme convergentie van een rij functies . . . . . . . . . . . . . . . . . . . . . .
9
2.2
Uniforme convergentie van een reeks functies . . . . . . . . . . . . . . . . . . . . 14
2.3
Machtreeksen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
2.4
Taylorreeksen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
2.5
Goniometrische reeksen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
2.6
De Fourierintegraal . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
Gewone en parti¨ele differentiaalvergelijkingen Differentiaalvergelijkingen
48 49
3.1
Definitie en voorbeelden . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
3.2
Beginvoorwaarden en randvoorwaarden . . . . . . . . . . . . . . . . . . . . . . . 53
Differentiaalvergelijkingen van eerste orde
55
4.1
Juiste differentiaalvergelijkingen . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
4.2
Homogene differentiaalvergelijkingen . . . . . . . . . . . . . . . . . . . . . . . . 60
4.3
Lineaire differentiaalvergelijkingen van eerste orde . . . . . . . . . . . . . . . . . 63
1
5
6
7
8
9
III
Normale differentiaalstelsels en vergelijkingen
67
5.1
De methode van afleiding en eliminatie . . . . . . . . . . . . . . . . . . . . . . . 67
5.2
Lineaire differentiaalstelsels . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
5.3
Lineaire differentiaalvergelijkingen van orde n . . . . . . . . . . . . . . . . . . . . 73
5.4
Constante co¨effici¨enten . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76
5.5
De methode der eigenwaarden en eigenvectoren . . . . . . . . . . . . . . . . . . . 86
Oplossing van differentiaalvergelijkingen door reeksontwikkeling
87
6.1
Inleiding . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87
6.2
Oplossing in een omgeving van een gewoon punt . . . . . . . . . . . . . . . . . . 88
6.3
Oplossing in een omgeving van een regelmatig singulier punt . . . . . . . . . . . . 93
Eerste integralen
99
7.1
Eerste integralen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99
7.2
Het bepalen van eerste integralen . . . . . . . . . . . . . . . . . . . . . . . . . . . 101
Lineaire parti¨ele differentiaalvergelijkingen van orde 1
104
8.1
Inleiding . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104
8.2
Homogene lineaire parti¨ele differentiaalvergelijking van orde e´ e´ n . . . . . . . . . . 106
8.3
De volledige lineaire vergelijking van orde e´ e´ n . . . . . . . . . . . . . . . . . . . 108
8.4
Het vraagstuk van Cauchy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110
Variatierekening
113
9.1
De vergelijkingen van Euler . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113
9.2
Het vraagstuk met nevenvoorwaarden . . . . . . . . . . . . . . . . . . . . . . . . 121
Numerieke Analyse
125
10 Het oplossen van vergelijkingen en stelsels
126
10.1 Algoritmen voor het oplossen van vergelijkingen . . . . . . . . . . . . . . . . . . 126 10.2 De orde van een iteratieve methode . . . . . . . . . . . . . . . . . . . . . . . . . . 131 10.3 Het oplossen van stelsels . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133
2
10.4 Stelsels lineaire vergelijkingen . . . . . . . . . . . . . . . . . . . . . . . . . . . . 137 11 Afleiden en integreren
140
11.1 De interpolatieformule van Lagrange . . . . . . . . . . . . . . . . . . . . . . . . . 140 11.2 Numeriek Afleiden . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142 11.3 Numeriek integreren . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 145 11.4 Differentiaalvergelijkingen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 152
3
Deel I Reeksen
4
Hoofdstuk 1 Rijen en reeksen 1.1
Numerieke rijen
In deze paragraaf herhalen we kort enkele belangrijke eigenschappen van limieten van rijen. Definitie 1.1.1 lim un = l ⇐⇒ ∀ε > 0, ∃N : n > N ⇒ |un − l| < ε
n→∞
We noemen l de limiet van de rij (un ), en we zeggen dat de rij convergent is. Een rij die niet convergent is, noemen we divergent. Op dezelfde manier defini¨eren we lim un = +∞ ⇐⇒ ∀α ∈ R, ∃N : n > N ⇒ un > α
n→∞
en lim un = −∞ ⇐⇒ ∀α ∈ R, ∃N : n > N ⇒ un < α
n→∞
We zeggen dan dat (un ) divergeert naar +∞, resp. −∞. Voorbeelden 1.1.2 1 =0 n→∞ 2n lim n2 = +∞ n→∞ √ lim − n = −∞ lim
n→∞
lim (−1)n bestaat niet
n→∞
Indien een rij divergent is, en ook niet divergeert naar ±∞, dan noemen we de rij schommelend of oscillerend. Een rij heeft niet altijd een limiet (zie het voorbeeld hierboven). Wel is het zo dat de limiet, indien hij bestaat, uniek is. Stelling 1.1.3 Een convergente rij is begrensd.
5
Bewijs. Kies ε = 1 in de definitie. Dan geldt voor n > N : l − 1 < un < l + 1 zodat {uN+1 , uN+2 , uN+3 , · · ·} begrensd is. Aangezien {u1 , u2 , u3 , · · · uN } eindig en dus begrensd is, volgt dat de rij begrensd is. De omgekeerde eigenschap geldt niet: een begrensde rij kan divergent zijn: zie voorbeeld 4 hierboven. Wel geldt dat een begrensde monotone rij convergent is (zie stelling 1.1.5). Eerst hebben we een definitie nodig. Definitie 1.1.4 Neem een rij (un ). (un ) heet strikt dalend als ∀m, n ∈ N0 : m > n =⇒ um < un (un ) heet niet stijgend als ∀m, n ∈ N0 : m > n =⇒ um ≤ un (un ) heet strikt stijgend als ∀m, n ∈ N0 : m > n =⇒ um > un (un ) heet niet dalend als ∀m, n ∈ N0 : m > n =⇒ um ≥ un In elk van deze vier gevallen noemen we un een monotone rij. Stelling 1.1.5 Een niet dalende (niet stijgende) rij, die naar boven (naar onder) begrensd is, convergeert naar haar bovenste (onderste) grens. Bewijs. Onderstel un niet dalend, en noteer m = supVun . We hebben de volgende eigenschap voor het supremum gezien: ∀ε > 0, ∃N : uN > m − ε Aangezien un niet dalend is, geldt voor elke n > N : m − ε < uN ≤ un ≤ m < m + ε of |un − m| < ε zodat lim un = m
n→∞
Het bewijs voor een niet stijgende rij verloopt analoog. Volledig analoog kunnen we aantonen:
6
Stelling 1.1.6 Een niet dalende (niet stijgende) rij die niet begrensd is, divergeert naar +∞ (−∞). Bewijs. We bewijzen alleen het geval waarin de rij (un ) niet dalend en niet begrensd is. Dan bestaat voor elke α ∈ R een index N zodat uN > α. Voor elke n > N geldt dan un ≥ uN > α. Opmerking 1.1.7 De stellingen 1.1.5 en 1.1.6 gelden ook voor rijen die slechts vanaf een zekere index N monotoon zijn. Het is soms belangrijk te kunnen bewijzen dat een rij convergeert, zonder daarom de limiet te hoeven kennen. Zulk een criterium zullen we nu opstellen. Definitie 1.1.8 Een rij (un ) is een Cauchyrij indien ∀ε > 0, ∃N > 0 : n, m > N =⇒ |un − um | < ε
(1.1)
Stelling 1.1.9 Een rij is convergent dan en slechts dan als ze een Cauchyrij is. Bewijs. Onderstel eerst dat (un ) convergent is, en stel lim un = l, m.a.w. n→∞
∀ε > 0, ∃N > 0 : n > N =⇒ |un − l| <
ε 2
Neem nu n, m > N. Dan geldt |un − um | ≤ |un − l| + |l − um | <
ε ε + =ε 2 2
Dit bewijst een implicatie. Het bewijs van de omgekeerde implicatie is moeilijker, en we laten dat hier achterwege.
1.2
Numerieke reeksen
Definitie 1.2.1 Gegeven is een numerieke rij (un ). Een uitdrukking van de vorm u1 + u2 + u3 + · · · + un + · · · of, korter, ∞
∑ un
n=1
noemt men een numerieke reeks. Stel n
sn = u1 + u2 + · · · + un = ∑ ui i=1
7
∞
Men noemt de numerieke rij (sn ) de rij der parti¨ele sommen van de reeks
∑ un. Indien de rij der
n=1
∞
parti¨ele sommen convergeert naar een getal s ∈ R, dan zegt men dat de reeks
∑ un convergent is
n=1
met som s, en met noteert dit als volgt ∞
∑ un = s
n=1
Dit betekent dus dat ∀ε > 0, ∃N ∈ N : ∀n > N : |sn − s| < ε Voor de parti¨ele sommen van een reeks hebben we volgende formule: sn+k − sn = un+1 + un+2 + · · · + un+k . Het convergentiecriterium van Cauchy, toegepast op de rij der parti¨ele sommen neemt daarom volgende vorm aan. Stelling 1.2.2 (convergentiecriterium van Cauchy) Een numerieke reeks ∑∞ n=1 un is convergent dan en slechts dan als ∀ε > 0, ∃N : ∀n > N, k ≥ 1 : |un+1 + un+2 + · · · + un+k | < ε ∞
Definitie 1.2.3 Een reeks
∑ un wordt een positieve reeks genoemd indien un > 0 voor elke n.
n=1 ∞
Indien de reeks
∑ un positief is, dan geldt duidelijk dat s1 < s2 < s3 < · · · en de rij van de parti¨ele
n=1
sommen is dus stijgend. We kunnen dus concluderen: Stelling 1.2.4 (basiskenmerk voor positieve reeksen) ∞
Een positieve reeks
∑ un is convergent als en alleen als de rij der parti¨ele sommen (sn) begrensd
n=1
is. Als de rij der parti¨ele sommen niet begrensd is, dan divergeert de reeks naar +∞.
8
Hoofdstuk 2 Reeksen van functies 2.1
Uniforme convergentie van een rij functies
Beschouw een rij functies un : (a, b) → R. Voor elke x ∈ (a, b) kunnen we de numerieke rij (un (x)) bekijken. We kunnen dan nagaan of deze rij een limiet heeft. Vandaar de volgende definitie. Definitie 2.1.1 De rij functies un convergeert puntsgewijs naar de functie ` : (a, b) → R indien voor elke x ∈ (a, b) geldt dat lim un (x) = `(x) n→∞
of, met andere woorden, indien ∀ε > 0, ∀x ∈ (a, b), ∃Nx : n > Nx =⇒ |un (x) − `(x)| < ε
(2.1)
Een analoge definitie kunnen we opschrijven voor rijen van functies gedefinieerd op een gesloten interval. Puntsgewijze limieten kunnen gemakkelijk berekend worden: beschouw de veranderlijke x als een parameter, en bereken de limiet alsof het de limiet van een gewone numerieke rij betreft. Puntsgewijze limieten hebben echter geen mooie eigenschappen. Zo is het mogelijk dat de puntsgewijze limiet van een rij continue functies niet langer continu is. Dit blijkt uit de volgende voorbeelden. Voorbeelden 2.1.2 1) Beschouw de rij functies un : R → R gedefinieerd als volgt: 1 −1 als x < − n ; un (x) = nx als − 1n < x < 1n ; 1 als x > n1 . De functies un zijn overal continu. De puntsgewijze limiet is echter niet continu in 0: lim un (x) = `(x)
n→∞
9
met
−1 `(x) = 0 1
als x = 0; als x > 0.
y=u_5(x)
2 1.5
1.5
1
1
0.5
0.5
0
0
-0.5
-0.5
-1
-1
-1.5
-1.5
-2 -2
-1.5
-1
-0.5
0
y=l(x)
2
y-as
y-as
als x < 0;
0.5
1
1.5
-2 -2
2
-1.5
-1
x-as
-0.5
0
0.5
1
1.5
2
x-as
Figuur 2.1: Een puntsgewijze limiet van continue functies is niet noodzakelijk continu 2) De puntsgewijze limiet van de rij functies un : R → R gedefinieerd door un (x) = 1 + x2 −
1 (1 + x2 )n
wordt gegeven door `(x) = lim un (x) = n→∞
1 + x2
als x 6= 0;
0
als x = 0.
Om dat te verhelpen is een nieuwe definitie van convergentie van een rij functies nodig. Merk op dat in (2.1) de index Nx in het algemeen afhangt van x. Indien Nx onafhankelijk van x kan gekozen worden, dan zeggen we dat de rij functies uniform convergeert. Definitie 2.1.3 De rij functies un convergeert uniform naar de functie ` : (a, b) → R indien ∀ε > 0, ∃N : n > N =⇒ ∀x ∈ (a, b) : |un (x) − `(x)| < ε
(2.2)
Een analoge definitie kunnen we opschrijven voor rijen van functies gedefinieerd op een gesloten interval. Meetkundig kunnen we dit als volgt zien: voor elke ε > 0 en voor n groot genoeg ligt de grafiek van de kromme y = un (x) tussen de grafieken van de functies y = `(x) − ε en y = `(x) + ε. Deze eigenschap geldt duidelijk niet voor de twee bovenstaande voorbeelden.
10
y=u_10(x)
9
9
8
8
7
7
6
6
5
5
4
4
3
3
2
2
1
1
0 -3
-2
-1
0
y=u_50(x)
10
y-as
y-as
10
1
2
0 -3
3
-2
-1
x-as y=u_500(x)
9
9
8
8
7
7
6
6
5
4
3
3
2
2
1
1 -2
-1
0
2
3
1
2
3
5
4
0 -3
1
y=l(x)
10
y-as
y-as
10
0 x-as
1
2
0 -3
3
x-as
-2
-1
0 x-as
Figuur 2.2: Een puntsgewijze limiet van continue functies is niet noodzakelijk continu Stelling 2.1.4 Als een rij functies un (x) uniform convergeert naar `(x) over een interval (a, b), en indien de functies un (x) continu zijn over (a, b) vanaf een zekere index M, dan is `(x) ook continu over (a, b). Een analoge eigenschap geldt voor gesloten intervallen. Bewijs. Kies ε > 0 willekeurig, en x ∈ (a, b). Dan geldt: |`(x + h) − `(x)| ≤ |`(x + h) − un (x + h)| + |un (x + h) − un (x)| + |un (x) − `(x)| Omdat (un (x)) uniform convergeert, bestaat er een index N zo dat voor elke y ∈ (a, b) en n > N geldt dat ε |un (y) − `(y)| < 3 in het bijzonder ε ε |`(x + h) − un (x + h)| < en |un (x) − `(x)| < (2.3) 3 3 Kies nu n vast en groter dan de indices M en N. Dan is un continu in x, en er bestaat een δ > 0 zo dat ε (2.4) |h| < δ =⇒ |un (x + h) − un (x)| < 3
11
l(x) + e l(x) l(x) – e
a
b
De kromme un(x) ligt overal dicht bij l(x)
x
Figuur 2.3: Uniforme convergentie Combineren van (2.3) en (2.4) geeft |h| < δ =⇒ |`(x + h) − `(x)| < ε en dus is ` continu in x.
Stelling 2.1.5 Onderstel dat un : [a, b] → R continu is, en dat un (x) uniform convergeert naar l(x) over [a, b]. Dan geldt Z b
lim
n→∞ a
Z b
un (x)dx =
l(x)dx a
Met andere woorden, we mogen dan limiet en integraal verwisselen. Bewijs. Uit stelling 2.1.4 volgt dat l continu is over [a, b], en bijgevolg kan ge¨ıntegreerd worden. Omdat un (x) uniform convergeert naar l(x) hebben we ∀ε > 0, ∃N : ∀n > N, ∀x ∈ [a, b] : |l(x) − un (x)| <
ε . b−a
Voor alle n > N hebben we dus Z b Z b Z b un (x)dx − l(x)dx = (un (x) − l(x))dx a a a ≤
Z b a
|un (x) − l(x)|dx <
ε(b − a) =ε b−a
Stelling 2.1.6 Onderstel dat un : (a, b) → R continu afleidbaar is, dat un (x) puntsgewijs convergeert naar l(x) over [a, b], en u0n (x) uniform naar f (x) over (a, b). Dan is l(x) afleidbaar, en l 0 (x) = f (x) = lim u0n (x) n→∞
Met andere woorden, we mogen dan limiet en afgeleide verwisselen.
12
Bewijs. Omdat u0n uniform convergeert naar f over elk gesloten deelinterval [x0 , x] van (a, b) hebben we, vanwege stelling 2.1.5: Z x
Z x
f (t)dt = lim x0
n→∞ x0
u0n (t)dt = lim un (x) − un (x0 ) = l(x) − l(x0 ) n→∞
Beide leden afleiden naar x geeft d f (x) = dx
Z x
f (t)dt = l 0 (x).
x0
Net zoals in het geval van een numerieke rij bestaat er een criterium voor uniforme convergentie van een rij functies waarin de limietfunctie `(x) niet voorkomt: Stelling 2.1.7 Criterium van Cauchy De rij functies un convergeert uniform over het interval (a, b) dan en slechts dan als ∀ε > 0, ∃N : m, n > N =⇒ ∀x ∈ (a, b) : |un (x) − um (x)| < ε
(2.5)
Een analoge eigenschap kunnen we opschrijven voor rijen van functies gedefinieerd op een gesloten interval. Bewijs. Onderstel dat un uniform convergeert over het interval (a, b) naar een functie `(x). Voor elke ε > 0 bestaat dan een index N zo dat voor n > N en x ∈ (a, b) geldt dat |un (x) − `(x)| <
ε 2
Voor n, m > N en x ∈ (a, b) geldt dus: |un (x) − um (x)| ≤ |un (x) − `(x)| + |`(x) − um (x)| <
ε ε + =ε 2 2
en dus is voldaan aan (2.5). Omgekeerd, onderstel dat de rij (un ) voldoet aan (2.5). Dan is voor elke x ∈ (a, b) de numerieke rij (un (x)) een Cauchyrij. De rij (un (x)) is dus convergent; noem de limiet `(x). un (x) convergeert dan puntsgewijs naar `(x), en de stelling is bewezen als we kunnen aantonen dat deze convergentie ook uniform is. Bij onderstelling is |un (x) − um (x)| < ε zodra n, m > N, en ongeacht x. Dus geldt voor elke x ∈ (a, b): lim |un (x) − um (x)| ≤ ε
m→∞
en dus |un (x) − `(x)| ≤ ε < 2ε voor elke n > N en x ∈ (a, b). De convergentie is dus uniform.
13
Opmerking 2.1.8 Voor een vector ~a ∈ Rn kan men de lengte (of norm) invoeren: s n
∑ a2i
k~ak =
i=1
Voor een functie die begrensd is over een interval [a, b] kan men ook een norm defini¨eren. Men doet dit als volgt: De norm van een begrensde functie f : [a, b] → R wordt gegeven door de formule k f k = sup{| f (x)| : a ≤ x ≤ b}
(2.6)
De definitie van uniforme convergentie van een rij functies kan nu als volgt herschreven worden: de rij functies un convergeert uniform over [a, b] naar de functie ` : [a, b] → R indien ∀ε > 0, ∃N : n > N =⇒ kun − `k < ε
(2.7)
Herschrijf zelf het convergentiecriterium van Cauchy met behulp van de norm.
2.2
Uniforme convergentie van een reeks functies
Zij A een deel van R. We zullen meestal aannemen dat A convex is, d.w.z. A is een (open, halfopen of gesloten) interval, een (open of gesloten) halve rechte, of gans R. Beschouw een rij functies un : A ⊂ R → R. De uitdrukking ∞
∑ un(x)
n=1
noemen we een reeks van functies. De n-de parti¨ele som n
sn (x) = ∑ ui (x) i=1
is dan ook een functie van x. ∞
We zeggen dat
∑ un(x) puntsgewijs convergent is over A met som s(x) als
n=1
lim sn (x) = s(x)
n→∞
voor elke x ∈ A, of, met andere woorden, indien ∀ε > 0, ∀x ∈ A, ∃Nε,x : n > Nε,x =⇒ |s(x) − sn (x)| < ε ∞
Om de puntsgewijze convergentie van de reeks
∑ un(x) te onderzoeken volstaat het dus om voor
n=1
∞
elke x0 ∈ A de convergentie van de numerieke reeks
∑ un(x0).
n=1
14
Voorbeeld 2.2.1 De reeks 1 + x + x2 + x3 + · · · convergeert puntsgewijs over (−1, 1). De reekssom is s(x) = 1/(1 − x). Een reeks convergeert puntsgewijs indien de rij van de parti¨ele sommen puntsgewijs convergeert. In § 2.1 hebben we een sterkere vorm van convergentie van rijen van functies besproken: uniforme convergentie. Uniforme convergentie heeft een aantal interessante eigenschappen: de uniforme limiet van een rij continue functies is opnieuw continu, en deze eigenschap geldt niet voor puntsgewijze convergentie. ∞
∑ un(x) uniform convergent over A met som s(x) indien de rij der parti¨ele
We noemen de reeks
n=1
sommen (sn (x)) uniform convergeert naar s(x) over A, of, met andere woorden, indien ∀ε > 0, ∃Nε : ∀x ∈ A : n > Nε =⇒ |s(x) − sn (x)| < ε We hebben onmiddellijk de volgende eigenschappen. ∞
Stelling 2.2.2 Als
∑ un(x) uniform convergeert naar s(x) over het interval (a, b), en de functies
n=1
un zijn continu over (a, b), dan is s(x) continu over (a, b). Bewijs. s(x) is de uniforme limiet van de rij continue functies sn (x), en is dus continu vanwege stelling 2.1.4. ∞
Stelling 2.2.3 Onderstel dat de functies un : [a, b] → R continu zijn, en dat
∑ un(x) = s(x) uniform
n=1
convergent is over [a, b]. Dan geldt dat Z b
∞ Z b
s(x)dx =
∑
a
n=1 a
un (x)dx
met andere woorden, een uniform convergente reeks mag term per term ge¨ıntegreerd worden.
Bewijs. Dit volgt onmiddellijk uit stelling 2.1.5.
Stelling 2.2.4 Onderstel dat de functies un : (a, b) → R een continue afgeleide bezitten, en dat ∞
de reeks
∞
∑ un(x) = s(x) puntsgewijs convergeert over (a, b). Indien de reeks
n=1
∑ u0n(x) uniform
n=1
∞
convergeert over (a, b), dan is
∑ u0n(x) = s0(x), m.a.w. men mag de reeks term per term afleiden.
n=1
Bewijs. Dit volgt onmiddellijk uit stelling 2.1.6. ∞
Stelling 2.2.5 (criterium van Cauchy) De reeks
∑ un(x) convergeert uniform over A als en al-
n=1
leen als ∀ε > 0, ∃Nε : ∀x ∈ A, ∀n > Nε , ∀p > 0 : |un+1 (x) + un+2 (x) + · · · + un+p (x)| < ε
15
Bewijs. Dit volgt onmiddellijk uit stelling 2.1.7.
In de praktijk is het criterium van Cauchy niet erg handig. In de volgende twee stellingen geven we twee meer eenvoudige criteria. Stelling 2.2.6 (criterium van Weierstrass) Beschouw een rij functies un : A ⊂ R → R. Indien |un (x)| ≤ mn voor elke n > N en x ∈ A, en indien ∑ mn een convergente numerieke reeks is, dan is ∑ un (x) uniform convergent over A. Bewijs. Omdat ∑ mn convergeert, geldt, vanwege stelling 1.2.2 ∀ε > 0, ∃Nε : ∀n > Nε , p > 0 : |mn+1 + mn+2 + · · · + mn+p | < ε Voor elke n > max{N, Nε } en x ∈ A geldt nu dat |un+1 (x) + un+2 (x) + · · · + un+p (x)| ≤ mn+1 + mn+2 + · · · + mn+p < ε en uit stelling 2.2.5 volgt dat ∑ un (x) uniform convergent over A.
∞
Stelling 2.2.7 Als
∑ un(x) = s(x) uniform convergeert over A, en f :
A ⊂ R → R een begrensde
n=1 ∞
functie is, dan is
∑ f (x)un(x) uniform convergent met som f (x)s(x).
n=1
Bewijs. Onderstel dat | f (x)| < m, voor elke x ∈ A. Voor elke ε > 0 bestaat een index N zodat voor elke n > N en x ∈ A geldt dat ε |s(x) − sn (x)| < m en hieruit volgt dat ε | f (x)s(x) − f (x)sn (x)| < m = ε m ∞
en dit betekent dat
∑ f (x)un(x) uniform convergent is met som f (x)s(x).
n=1
2.3
Machtreeksen
Definitie 2.3.1 Een machtreeks is een reeks van functies van de vorm ∞
∑ an(x − a)n = a0 + a1(x − a) + a2(x − a)2 + · · ·
n=0
waarbij a0 , a1 , a2 , · · · ∈ R.
16
Voor wat de convergentie van een machtreeks betreft volstaat het het geval a = 0 te bestuderen. Immers, als we de substitutie u = x − a uitvoeren, dan vinden we ∞
∞
n=0
n=0
∑ an(x − a)n = ∑ anun
Lemma 2.3.2 Onderstel dat r > 0 en lim an rn = 0 en |x| < r
n→∞ ∞
Dan is
∑ anxn absoluut convergent.
n=0
Bewijs. Er bestaat een index N zodat voor elke n > N geldt dat |an rn | < 1. Kies α ∈ [0, 1) zodat |x| < αr. Dan geldt voor elke n > N dat |an xn | < |an |αn rn < αn ∞
De meetkundige reeks
∑ αn is convergent, en uit het vergelijkend criterium volgt dat
n=0
Onderstel dat lim an rn = 0
n→∞
Indien 0 ≤ s ≤ r, dan geldt voor elke n dat 0 ≤ |an |sn ≤ |an |rn en dus is ook lim an sn = 0
n→∞ ∞
Stelling 2.3.3 Beschouw een machtreeks
∑ anxn, en stel
n=0
R = sup{r ≥ 0| lim an rn = 0} n→∞
Dan geldt ∞
∑ anxn absoluut convergent;
n=0 ∞
2. |x| > R =⇒
∑ |anxn|
n=0
convergent is.
1. |x| < R =⇒
∞
∑ anxn divergent.
n=0
17
Bewijs. 1) Kies r zodat |x| < r < R. Uit de opmerking hierboven volgt dat lim an rn = 0
n→∞ ∞
en uit lemma 2.3.2 volgt dat
∑ anxn absoluut convergeert.
n=0
2) Als |x| > R, dan is het onmogelijk dat lim an |xn | = 0, want dat zou impliceren dat n→∞
|x| ≤ sup{r ≥ 0| lim an rn = 0} = R n→∞
Dus is lim an xn 6= 0, en de reeks is divergent.
n→∞
De machtreeks convergeert dus absoluut op het interval (−R, R), en divergeert op (−∞, −R) ∪ (R, +∞). Enkel in de punten R en −R is er nog onzekerheid. Alles is er mogelijk: convergentie in beide punten, divergentie in beide punten, en convergentie in het een en divergentie in het ander. R noemen we de convergentiestraal van de machtreeks, en (−R, R) het convergentieinterval. In de volgende twee stellingen zullen we formules zien die toelaten om in veel gevallen R te berekenen. ∞
Stelling 2.3.4 Beschouw een machtreeks
∑ anxn. Als
n=0
a n L = lim n→∞ an+1 bestaat, dan is L = R, de convergentiestraal van de machtreeks. Bewijs. De stelling volgt uit het convergentiekenmerk van d’Alembert (ook genaamd de verhoudingstest). Als |x| < L, dan is a x |x| a xn+1 n+1 n+1 = lim <1 lim = n n→∞ n→∞ an x an L ∞
en dit impliceert dat
∑ anxn absoluut convergeert.
n=0
∞
Op dezelfde manier volgt dat
∑ anxn divergeert als |x| > L. Het kan dus niet anders zijn dan dat
n=0
R = L.
∞
Stelling 2.3.5 Beschouw een machtreeks
∑ anxn. Als
n=0
Λ = lim
n→∞
p n |an |
bestaat, dan is 1/Λ = R, de convergentiestraal van de machtreeks.
18
Bewijs. De stelling volgt uit het convergentiekenmerk van Cauchy (ook genaamd de worteltest). Bewijs dit zelf als oefening. Voorbeelden 2.3.6 1) De machtreeks ∞
x2
xn
x3
∑ n! = 1 + x + 2 + 6 + · · ·
n=0
convergeert voor elke x ∈ R, want a (n + 1)! n = +∞ R = lim = lim n→∞ n→∞ an+1 n! 2) Beschouw de machtreeks ∞
∑ (n + 1)xn = 1 + 2x + 3x2 + 4x3 + · · ·
n=0
De convergentiestraal is a n+1 n R = lim =1 = lim n→∞ an+1 n→∞ n + 2 De reeks convergeert dus absoluut voor x ∈ (−1, 1). In de eindpunten x = ±1 is de reeks divergent. 3) Beschouw de machtreeks ∞ n x x2 x3 = x + + +··· ∑ 2 3 n=1 n De convergentiestraal is a 1/n n =1 R = lim = lim n→∞ an+1 n→∞ 1/n + 1 De reeks convergeert dus absoluut voor x ∈ (−1, 1). Voor x = 1 is de reeks divergent (harmonische reeks), en voor x = −1 is ze (relatief) convergent (alternerende harmonische reeks). Het convergentieinterval is dus [−1, 1). 4) Beschouw de machtreeks ∞
∑ nnxn = 1 + x + 4x2 + 27x3 + · · ·
n=0
Het invers van de convergentiestraal is p 1 = lim n |an | = lim n = +∞ n→∞ R n→∞ en de convergentiestraal R = 0. De machtreeks convergeert dus alleen voor x = 0.
19
Uniforme convergentie van machtreeksen ∞
Stelling 2.3.7 Onderstel dat R de convergentiestraal is van de machtreeks
∑ anxn. De machtreeks
n=0
convergeert uniform over elk gesloten interval [−r, r], met r < R. ∞
Bewijs. Kies r < R. De reeks
∑ |an|rn is convergent, en voor elke x ∈ [−r, r] hebben we dat
n=0
|an xn | ≤ |an |rn
Het gestelde volgt nu uit het kenmerk van Weierstrass (stelling 2.2.6). ∞
Gevolg 2.3.8 De som van een machtreeks s(x) =
∑ anxn is continu over (−R, R).
n=0
Bewijs. Onderstel x ∈ (−R, R), en neem r zodat |x| < r < R. De machtreeks is uniform convergent over [−r, r] (stelling 2.3.7) en dus is s continu over (−r, r) (stelling 2.2.2). In het bijzonder is s continu in x. Zonder bewijs vermelden we ook de volgende stelling. ∞
Stelling 2.3.9 (stelling van Abel) Indien de machtreeks s(x) =
∑ anxn convergeert in x = R (of
n=0
x = −R), dan is s(x) linkscontinu in R (rechtscontinu in −R). Afleiden en integreren van machtreeksen ∞
Lemma 2.3.10 De machtreeks
∑ nanxn−1 heeft dezelfde convergentiestraal als
n=1
∞
∑ an x n .
n=0
∞
Bewijs. Zij R de convergentiestraal van de machtreeks
∑ anxn. We weten dat
n=0
R = sup{r ≥ 0| lim an rn = 0} n→∞
∞
(stelling 2.3.3). Het volstaat nu aan te tonen dat de reeks
∑ nanxn−1 absoluut convergent is voor
n=1
|x| < R en divergent voor |x| > R. Onderstel eerst dat |x| < R, en kies r > 0 en α ∈ [0, 1) zo dat |x| < αr < r < R
20
Omdat lim an rn = 0
n→∞
bestaat een index N zodat voor elke n > N geldt dat |an |rn < 1 en dus n|an ||xn−1 | < n|an |αn−1 rn−1 <
nαn−1 r ∞
Omdat α < 1 is ∑ nαn−1 convergent, en uit het vergelijkend kenmerk volgt dat
∑ nanxn−1 absoluut
n=1
convergent is. Indien |x| > R, dan is lim an |x|n 6= 0
n→∞
en dus is lim nan |x|n−1 6= 0
n→∞ ∞
en bijgevolg is
∑ nanxn−1 divergent.
n=1
We kunnen nu aantonen dat machtreeksen binnen het convergentieinterval term per term mogen afgeleid en ge¨ıntegreerd worden. ∞
Stelling 2.3.11 Beschouw een machtreeks s(x) =
∑ anxn met convergentiestraal R.
Voor elke
n=0
x, x0 ∈ (−R, R) geldt dat 0
s (x) =
∞
∑ nanxn−1
n=1 ∞
Z x
s(t)dt = x0
(2.8)
an
∑ n + 1 (xn+1 − x0n+1)
(2.9)
n=0
Bewijs. Kies r < R zo dat x, x0 ∈ (−r, r). Uit stelling 2.3.7 volgt dat de machtreeks uniform convergent is over [−r, r], en a fortiori over (−r, r). (2.9) volgt nu onmiddellijk uit stelling 2.2.3. ∞
In stelling 2.3.10 hebben we gezien dat R ook de convergentiestraal van de reeks
∑ nanxn is. Uit
n=1
stelling 2.3.7 volgt nu dat ook deze reeks uniform convergent is over (−r, r), en (2.8) volgt nu onmiddellijk uit stelling 2.2.4. ∞
Gevolg 2.3.12 s(x) =
∑ anxn is onbeperkt afleidbaar over het interval (−R, R). De i-de afgeleide
n=0
is de som van de machtreeks s(i) (x) =
∞
(n + i)! an+i xn n! n=0
∑
21
Voorbeelden 2.3.13 1) Voor elke x ∈ (−1, 1) geldt dat 1 = 1 + x + x2 + x3 + · · · 1−x en (vervang x door −x): 1 = 1 − x + x2 − x3 + · · · 1+x We mogen term per term integreren, en dus ln (1 + x) = x −
(2.10)
x2 x3 x4 + − +··· 2 3 4
De convergentiestraal van deze laatste machtreeks is dus eveneens 1. Voor x = −1 is de reeks divergent (harmonische reeks), en voor x = 1 is ze convergent (alternerende harmonische reeks). Het convergentieinterval van de machtreeks is dus (−1, 1] en uit de stelling van Abel volgt dat ∞ 1 1 1 (−1)n+1 ln 2 = 1 − + − + · · · = ∑ 2 3 4 n n=1
2) Uit (2.10) volgt, na vervanging van x door x2 , dat voor elke x ∈ (−1, 1): 1 = 1 − x2 + x4 − x6 + · · · 2 1+x Integreren geeft dat x3 x5 x7 + − +··· 3 5 7 Het convergentieinterval van deze machtreeks is [−1, 1]. Uit de stelling van Abel volgt dus dat bgtg x = x −
∞ π 1 1 1 (−1)n = 1− + − +··· = ∑ 4 3 5 7 n=0 2n + 1
3) Term per term afleiden van de meetkundige reeks geeft ∞ 1 2 3 4 = 1 + 2x + 3x + 4x + 5x + · · · = ∑ (n + 1)xn 2 (1 − x) n=0
voor elke x ∈ (−1, 1). Men kan deze formule ook rechtstreeks bewijzen. Doe dit zelf als oefening. Vermenigvuldiging van machtreeksen ∞
Eerst geven we nog een stelling over absoluut convergente reeksen. Onderstel dat
∞
∑ an en ∑ bn
n=0
n=0
twee numerieke reeksen zijn. Het product van deze twee reeksen is per definitie de reeks met algemene term n
cn = a0 bn + a1 bn−1 + a2 bn−2 + · · · + an b0 = ∑ ai bn−i i=0
22
met andere woorden, ∞
∞
∑ cn = a0b0 + a0b1 + a1b0 + a0b2 + a1b1 + a2b0 + · · · =
n=0
n
∑ ∑ aibn−i
n=0 i=0
Zowel de haakjes als de volgorde van de termen is van belang. ∞
Stelling 2.3.14 (stelling van Cauchy) Als de reeksen
∞
∑ an en ∑ bn absoluut convergent zijn, dan
n=0
n=0
∞
is ook de productreeks
∑ cn absoluut convergent, en de som van de productreeks is het product
n=0
van de sommen van beide reeksen. ∞
∞
∞
∑ cn = ∑ an ∑ bn
n=0
n=0
n=0
Bewijs. Dit bewijs is opgenomen voor de liefhebbers. De anderen kunnen opgelucht adem halen. ∞
We bewijzen eerst dat
∑ cn absoluut convergent is. Schrijf
n=0
n
n
i=0
i=0
s0n = ∑ |ai | en tn0 = ∑ |bi | ∞
∑ |cn| begrensd zijn.
We zullen aantonen dat de parti¨ele sommen van de reeks
n=0 n
j a b ≤ i j−i ∑∑ n
∑ |c j | =
j=0
j=0 i=0 n
≤
∑
n
∑ |ai||bk | =
j
n
∑ ∑ |aib j−i| = ∑
j=0 i=0
i=0 k=0
|ai bk |
0≤i,k≤n i+k≤n
n
n |a | ( |b | ∑ i ∑ k = s0ntn0
i=0
k=0
Omdat de parti¨ele sommen s0n en tn0 zijn begrensd, zijn de parti¨ele sommen van de productreeks ∞
∑ |cn| begrensd,en de productreeks is absoluut convergent.
n=0
Schrijf nu n
n
sn = ∑ ai en tn = i=0
zodat
∑ bk k=0
n
sn tn = ∑
n
∑ ai bk
i=0 k=0
We hebben ook dat 2n
2n
j
∑ c j = ∑ ∑ aib j−i = ∑
j=0
j=0 i=0
0≤i,k≤2n i+k≤2n
23
ai bk
en dus is 2n ∑ c j − sntn =
ai bk +
∑
j=0
n
≤
|ai bk | +
∑
n
≤
∑ n
∑ |aibk | +
i=n+1 k=0 2n
= (
|ai |)s0n + (
∑
|ai bk |
∑
n
∑
ai bk
i=n+1
n
∑ ∑ |aibk |
k=n+1 i=0 2n
∑
|bi |)tn0
i=n+1
Voor n groot genoeg worden 2n
∑
2n
|ai | en
i=n+1
willekeurig klein (criterium van Cauchy);
∑
|bi |
i=n+1
s0n
en
tn0
zijn begrensd, zodat 2n ∑ c j − sntn j=0
willekeurig klein kan gemaakt worden door n groot genoeg te nemen. Hieruit volgt dat ∞ ∞ ∞ c = lim s t = lim s lim t = nn n n ∑ n ∑ an ∑ bn n=0
n→∞
n→∞
n→∞
n=0
n=0
Onderstel nu dat twee machtreeksen ∞
∑ anxn en
n=0
∞
∑ bnxn
n=0
gegeven zijn. Het product van deze twee reeksen is de machtreeks ∞
n
∑ ∑ aibn−i
n x = a0 b0 + (a0 b1 + a1 b0 )x + (a0 b2 + a1 b1 + a2 b0 )x2 + · · ·
n=0 i=0
Uit stelling 2.3.14 volgt onmiddellijk Stelling 2.3.15 De convergentiestraal van het product van twee machtreeksen is tenminste gelijk aan het minimum van de convergentiestralen van deze twee machtreeksen. Bewijs. Onderstel van niet. Neem x zodanig dat |x| groter is dan de convergentiestraal van de productreeks, maar kleiner dan de convergentiestralen van de twee gegeven reeksen. Uit stelling 2.3.14 volgt dat de productreeks absoluut convergeert in x, en dit is strijdig met de onderstelling dat |x| groter is dan de convergentiestraal van de productreeks.
24
2.4
Taylorreeksen
In de vorige paragraaf hebben we gezien dat de som van een machtreeks een onbeperkt afleidbare functie op het convergentieinterval van de machtreeks is. In enkele gevallen konden we deze functie bepalen. We bekijken nu het omgekeerde probleem: gegeven is een functie f (x) die onbeperkt afleidbaar is in een omgeving van het punt a. Kan deze som geschreven worden als een machtreeks? Om dit probleem op te lossen herhalen we de formule van Taylor, zie volume 1. Stelling 2.4.1 (formule van Taylor) Onderstel dat f : (a − r, a + r) → R een onbeperkt afleidbare functie is en dat x ∈ (a − r, a + r). Voor elke n ∈ N bestaat een ξn (x) ∈ (a, x) zodanig dat f (x) = sn (x) + rn (x) met n f (i) (a) (x − a)i sn (x) = ∑ i! i=0 met rn (x) =
f (n+1) (ξn (x)) (x − a)n+1 (n + 1)!
sn (x) is de n-de parti¨ele som van de machtreeks ∞
∑
n=0
f (n) (a) (x − a)n n!
Deze machtreeks noemt men de Taylorreeks van f in het punt a. Als lim rn (x) = 0
n→∞
dan is lim sn (x) = f (x)
n→∞
of ∞
f (x) =
∑
n=0
f (n) (a) (x − a)n n!
en in dit geval kan f (x) geschreven worden als een machtreeks. Indien de Taylorreeks van de functie f convergeert naar f (x) voor x in een omgeving van het punt a, dan zeggen we dat f analytisch is in het punt a. In hetgeen volgt zullen we de Taylorreeks van enkele belangrijke functies bestuderen. De exponenti¨ele functie We schrijven de formule van Taylor op voor de functie f (x) = ex in het punt 0. We vinden gemakkelijk dat f (n) (x) = ex en f (n) (0) = 1 zodat ex = sn (x) + rn (x) met sn (x) = 1 + x +
x2 x3 xn + +···+ 2! 3! n!
25
en rn (x) =
xn+1 ξn e n + 1!
met ξn ∈ (0, x). Als x ≥ 0, dan is eξn < ex , en als x ≤ 0, dan is eξn < 1. In beide gevallen geldt |rn (x)| ≤ zodat
|x|n+1 |x| n→∞ e −→0 n + 1!
xn x2 x3 = 1+x+ + +··· e =∑ 2! 3! n=0 n! ∞
x
(2.11)
voor elke x ∈ R. De hyperbolische functies Als we x vervangen door −x in (2.11), dan vinden we −x
e
x2 x3 (−1)n xn = 1−x+ − +··· n! 2! 3! n=0 ∞
=
∑
(2.12)
voor elke x ∈ R. Door (2.11) en (2.12) op te tellen en af te trekken en te delen door 2 vinden we sh x = x +
∞ x3 x5 x2n+1 + +··· = ∑ 3! 5! n=0 (2n + 1)!
(2.13)
ch x = 1 +
∞ x2 x4 x2n + +··· = ∑ 2 4! n=0 (2n)!
(2.14)
De goniometrische functies Als we de formule van Taylor opschrijven voor de functie f (x) = sin x vinden we sin x = sn (x) + rn (x) met (−1)n−1 x2n−1 x3 x5 s2n (x) = x − + + · · · + 3! 5! (2n − 1)! en x2n+1 r2n (x) = (−1)n cos(ξ2n ) (2n + 1)! met ξ2n ∈ (0, x). We zien gemakkelijk dat |x|2n+1 n→∞ |r2n (x)| ≤ −→0 (2n + 1)! voor elke x ∈ R, zodat sin x = x −
∞ (−1)n x2n+1 x3 x5 + −··· = ∑ 3! 5! n=0 (2n + 1)!
26
(2.15)
voor elke x ∈ R. Op volledig analoge wijze vinden we cos x = 1 −
∞ x2 x4 (−1)n x2n + −··· = ∑ 2 4! n=0 (2n)!
(2.16)
De binomiaalreeks Neem m ∈ R en beschouw de functie f (x) = (1 + x)m Toepassing van de formule van Taylor levert f (x) = sn (x) + rn (x), met sn (x) = 1 + mx +
m(m − 1) 2 m(m − 1)(m − 2) · · · (m − n + 1) n x +···+ x 2 n!
en
m(m − 1)(m − 2) · · · (m − n) n+1 x (1 + θx)m−n−1 (n + 1)! Als m ∈ N, dan is rn = 0 voor n > m. We krijgen dan m m i m f (x) = (1 + x) = ∑ x i=0 i rn (x) =
en dit is niets anders dan het binomium van Newton. Immers, de binomiaalco¨effici¨enten worden gegeven door de formule m(m − 1)(m − 2) · · · (m − i + 1) m m! = = i!(m − i)! i! i Voor m ∈ R defini¨eren we nu de binomiaalco¨effici¨enten als volgt m m(m − 1)(m − 2) · · · (m − i + 1) = i i! De Taylorreeks van f (x) is dan m ∑ n xn n=0 ∞
De convergentiestraal van deze machtreeks is n+1 a n lim = lim =1 n→∞ m − n n→∞ an+1 Men kan aantonen dat lim rn (x) = 0 als −1 < x < 1 en we kunnen dus besluiten dat n→∞ ∞ m n m (1 + x) = ∑ x n=0 n voor −1 < x < 1. Voor m = 1/2 vinden we bijvoorbeeld √ x x2 3x3 3.5x4 1+x = 1+ − + − +··· 2 2.4 2.4.6 2.4.6.8
27
(2.17)
Samenvatting ∞ 1 = 1 + x + x2 + x3 + · · · = ∑ xn 1−x n=0
ln (1 + x) = x −
(−1 ≤ x < 1)
∞ x2 x3 x4 (−1)n+1 xn + − +··· = ∑ 2 3 4 n n=1
ex = 1 + x +
∞ n x x2 x3 + +··· = ∑ 2! 3! n=0 n!
(x ∈ R)
sh x = x +
∞ x3 x5 x2n+1 + +··· = ∑ 3! 5! n=0 (2n + 1)!
ch x = 1 +
∞ x2 x4 x2n + +··· = ∑ 2 4! n=0 (2n)!
sin x = x −
∞ x3 x5 (−1)n x2n+1 + −··· = ∑ 3! 5! n=0 (2n + 1)!
(x ∈ R)
(x ∈ R)
∞ x2 x4 (−1)n x2n + −··· = ∑ 2 4! n=0 (2n)! ∞ m n = ∑ x (−1 < x < 1) n=0 n
cos x = 1 − (1 + x)m
(−1 < x ≤ 1)
(x ∈ R) (x ∈ R)
Opmerkingen 2.4.2 1) Men zou kunnen denken dat, indien de Taylorreeks van f convergeert, dat ze dan automatisch naar de functie convergeert. Dit is niet noodzakelijk het geval. Beschouw bijvoorbeeld de functie ( −
f (x) =
e
1 x2
als x 6= 0;
als x = 0.
0 Voor x 6= 0 is f (n) (x) =
P(x) − 12 e x Q(x)
waarbij P en Q veeltermen zijn (verifieer dit zelf). Per inductie kunnen we dan aantonen dat f (n) (0) = 0 voor elke n. Immers, als f (n) (0) = 0, dan is P(h) − 12 e h h→0 hQ(h)
f (n+1) (0) = lim
Als we nu de formule van Taylor toepassen, dan vinden we sn (x) = 0 voor elke n. De Taylorreeks is dus 0, en de Taylorreeks convergeert, maar niet naar de functie f ! 2) In het algemeen zal de Taylorreeks sneller convergeren naarmate x dichter bij het punt a komt te liggen. Als we x ver van a nemen, dan kan de convergentie zeer traag zijn, en in vele gevallen
28
totaal ongeschikt voor praktische berekeningen. Neem bijvoorbeeld de reeksontwikkeling voor de sinusfunctie, en schrijf deze neer voor x = 1000: 106 109 + −··· 6 120 en dit is natuurlijk geen erg goede benadering voor een getal dat tussen −1 en 1 moet liggen. De 1500 algemene term in de reeks begint pas vanaf de 500ste term te dalen! De 500ste term is 10500! , en dit is een astronomisch groot getal. sin 1000 = 103 −
De voorbeelden die we besproken hebben zijn natuurlijk de eenvoudigste die we konden kiezen. Als een willekeurige functie f gegeven is, dan is het gewoonlijk niet mogelijk om een algemene formule voor de n-de afgeleide op te stellen, zodat we de algemene term in de Taylorreeks niet kunnen opschrijven. In veel gevallen kunnen we stelling 2.3.15 gebruiken om de Taylorreeks te bepalen. Voorbeelden 2.4.3 1) Voor elke x ∈ R hebben we x3 x5 x2 x3 ex sin x = 1 + x + + + · · · x − + − · · · 2! 3! 3! 5! 1 1 1 1 1 1 1 5 = x + x2 + − x3 + − x4 + − + x +··· 2 6 6 6 24 12 120 x3 x5 = x + x2 + − + · · · 3 30 2) Voor elke x ∈ R geldt x2 x4 cos x = 1 − + − · · · = 1 − u 2 4! met x2 x4 u = − +··· 2 4! Als x ∈ − π2 , π2 , dan is 0 < cos x ≤ 1, zodat u ∈ (0, 1) en 1 1 = = 1 + u + u2 + u3 + · · · cos x 1 − u en sin x cos x x3 x5 = x− + −··· 3! 5! x2 x4 x6 2 x2 x4 3 x2 x4 x6 1+ − + +··· + − + +··· + − +··· +··· 2 4! 6 2 4! 6 2 4! 1 1 1 1 1 1 = x+ − x3 + − − + + x5 2 6 24 12 120 4 3 5 x 2x = x+ + +··· 3 15
tg x =
29
π π voor elke x ∈ − 2 , 2 . 3) Beschouw de functie f (x) =
x ex − 1
Onderstel dat de taylorreeks van f co¨effici¨enten a0 , a1 , · · · heeft. Dan is f (x) =
x ex − 1
= a0 + a1 x + a2 x2 + · · ·
en
x2 x3 (a0 + a1 x + a2 x2 + · · ·) x + + + · · · = x 2 6 Identificeren van co¨effici¨enten in gelijke machten van x geeft a0 = 1 a0 a1 + = 0 2 1 a0 a1 + + a2 = 0 6 2 ···
Hieruit kunnen a0 , a1 , a2 , · · · recursief berekend worden. Men vindt x x2 x4 x = 1 − + − +··· ex − 1 2 12 720 Het nadeel van deze methode is dat we niet weten wanneer de Taylorreeks convergeert.
2.5
Goniometrische reeksen
Een goniometrische reeks is een reeks van de vorm ∞ a0 + ∑ (an cos nx + bn sin nx) 2 n=1
(2.18)
Opmerkingen 2.5.1 1) De factor 1/2 werd ingevoerd om technische redenen; sommige formules die we verderop zullen zien vereenvoudigen hierdoor een beetje. 2) In hoofdstuk 2 hebben we reeds machtreeksen besproken; hier krijgen we dus een tweede belangrijke klasse van reeksen van functies. We zullen zien dat sommige functies kunnen geschreven worden als de som van een goniometrische reeks. Merk op dat een parti¨ele som van de goniometrische reeks (2.18) kan geschreven worden als een veelterm in cos x en sin x (men spreekt van een goniometrische veelterm). 3) Alle termen in de reeks (2.18) zijn periodiek met periode 2π. Als de reeks (2.18) (puntsgewijs) convergeert, dan is ook de reekssom f (x) periodiek met periode 2π, dit wil zeggen f (x + 2kπ) = f (x)
30
voor elke x ∈ R en k ∈ Z. Het zal blijken dat goniometrische reeksen een geschikt instrument vormen om periodieke functies te bestuderen. 4) Men kan ook goniometrische reeksen met periode T (in plaats van periode 2π) bestuderen. Zulk een reeks is van de vorm ∞ a0 2nπt 2nπt + ∑ (an cos + bn sin ) (2.19) 2 n=1 T T Men noemt T1 de frequentie en ω = 2π T de pulsatie van de goniometrische reeks. (2.19) kan dus nog geschreven worden als ∞ a0 + ∑ (an cos nωt + bn sin nωt) (2.20) 2 n=1 (2.19) en (2.20) kunnen herleid worden tot (2.18) met de substitutie x = ωt =
2π t T
(2.20) kan ook nog als volgt geschreven worden. We stellen ρ0 = a0 /2 en an = ρn cos φn bn = ρn sin φn en we verkrijgen de reeks ∞
ρ0 + ∑ ρn cos(nωt − φn )
(2.21)
n=1
De termen ρn cos(nωt − φn ) noemt men de harmonischen. ρn noemt men de amplitude en φn de faze. 5) Een goniometrische reeks kan ook herschreven worden in complexe vorm: ∞ a0 + ∑ (an cos nx + bn sin nx) 2 n=1 ∞ einx − e−inx a0 einx + e−inx + ∑ (an + bn ) 2 n=1 2 2i ∞ a0 an − ibn inx an + ibn −inx = +∑ e + e 2 n=1 2 2
=
=
∞ ∞ a0 + ∑ (αn einx + αn e−inx ) = ∑ αn einx , 2 n=1 n=−∞
waarbij α−n = αn . Onderstel nu dat de reeks (2.18) uniform convergeert over een interval [a, a + 2π], en schrijf f voor de som van de reeks. We verkrijgen zo een functie f : R → R waarvan we weten dat ze periodiek is met periode 2π (zie opmerking 3) hierboven). Bovendien is f continu (zie stelling 2.2.2). We hebben ∞ a0 f (x) = + ∑ (an cos nx + bn sin nx) (2.22) 2 n=1
31
We zullen nu aantonen dat de co¨effici¨enten an en bn kunnen bepaald worden uit de functie f . Vermenigvuldig (2.22) met cos mx en integreer van a tot a + 2π. Gebruik makend van stelling 2.2.3 vinden we Z a+2π Z a+2π ∞ a0 f (x) cos mxdx = + ∑ (an cos nx + bn sin nx) cos mxdx 2 n=1 a a =
a0 2
Z a+2π
∞ Z a+2π
cos mxdx + ∑
a
n=1 a
(an cos nx + bn sin nx) cos mxdx
Aangezien voor m, n ∈ Z geldt (reken zelf na): Z a+2π
Z a+2π
sin nxdx =
cos nxdx = 0
Za a+2π
a
sin nx cos mxdx = 0 Za a+2π
Z a+2π
sin nx sin mxdx = a
vinden we
cos nx cos mxdx = πδnm
a
Z a+2π
f (x) cos mxdx = πam
a
Op analoge wijze vinden we (vermenigvuldig met sin mx) Z a+2π
f (x) sin mxdx = πbm
a
en (integreer f van a tot a + 2π): Z a+2π
f (x)dx = πa0
a
We kunnen dus besluiten: Stelling 2.5.2 Als de reeks (2.22) uniform convergeert, dan worden de co¨effici¨enten an en bn gegeven door de formules an =
1 π
bn =
1 π
Z a+2π
f (x) cos nxdx Za a+2π
(2.23) f (x) sin nxdx
a
Dankzij de factor 1/2 in (2.18) geldt de eerste formule van (2.23) zowel voor n = 0 als voor n > 0. Opmerking 2.5.3 De co¨effici¨enten αn van de goniometrische reeks in complexe vorm kunnen ook gemakkelijk bepaald worden: voor n ≥ 0 hebben we 1 2π i an − ibn = = f (x) cos nxdx − 2 2π 0 2π Z 2π 1 f (x)e−inx dx = 2π 0 Z 1 2π = f (x)einx dx 2π 0 Z
αn
α−n
32
Z 2π
f (x) sin nxdx 0
We hebben dus voor alle n ∈ Z: 1 αn = 2π
Z 2π
f (x)e−inx dx.
(2.24)
0
Onderstel nu dat f : R → R een functie is met periode 2π, en dat f een stuksgewijs continue functie is over het interval [a, a + 2π]. f hoeft dus niet noodzakelijk continu te zijn. We berekenen nu de co¨effici¨enten an en bn met behulp van de formules (2.23), en schrijven de reeks (2.18) op ∞ a0 + ∑ (an cos nx + bn sin nx) 2 n=1
Men noemt deze reeks de Fourierreeks van de functie f . A priori weten we niet of deze reeks al dan niet convergeert. En als de reeks convergeert, dan weten we ook niet of de reekssom f (x) is. Dat dit laatste niet steeds het geval is blijkt uit het volgende voorbeeld: neem een uniform convergente goniometrische reeks, en de reekssom f . Wijzig f in e´ e´ n punt a, en zet deze wijziging periodiek voort. De Fourierreeks van de aldus bekomen functie is nog steeds de oorspronkelijke reeks, maar in het punt a convergeert deze niet naar de functiewaarde. We gaan nu volgende vragen bestuderen. 1) Wanneer convergeert de reeks (2.18) puntsgewijs naar de functie f ? 2) Wanneer convergeert de reeks (2.18) uniform naar de functie f ? Het antwoord op de tweede vraag kan uiteraard alleen maar positief zijn als f een continue functie is. Opmerkingen 2.5.4 1) We zijn vertrokken van een functie f die periodiek is met periode 2π. Men kan uiteraard ook een stuksgewijs continue functie gedefinieerd op het interval [a, a + 2π] beschouwen en deze periodiek voortzetten. 2) Als de functie f een oneven functie is (m.a.w. f (−x) = − f (x)), dan bevat de Fourierreeks van f slechts termen in sinus. Immers, 1 an = π
Z π
f (x) cos nxdx = 0 −π
omdat de integraal de integraal is van een oneven functie over een symmetrisch interval. Voor de co¨effici¨enten bn geldt in dit geval 2 bn = π
Z π
f (x) sin nxdx 0
3) Als de functie f een even functie is (m.a.w. f (−x) = f (x)), dan bevat de Fourierreeks van f slechts termen in cosinus. Immers, 1 bn = π
Z π
f (x) sin nxdx = 0 −π
omdat de integraal de integraal is van een oneven functie over een symmetrisch interval. Voor de co¨effici¨enten an geldt in dit geval 2 an = π
Z π
f (x) cos nxdx 0
33
We zeggen dat een functie f voldoet aan een rechtervoorwaarde van Lipschitz in het punt a als de rechterlimiet f (a+) bestaat en als er een constante c > 0 en δ > 0 bestaan zodat 0 < h < δ =⇒ | f (a + h) − f (a+)| ≤ ch Op dezelfde manieren defini¨eren we de linkervoorwaarde van Lipschitz. f voldoet aan zulk een voorwaarde in a als f (a−) bestaat en als er een constante c > 0 en δ > 0 bestaan zodat 0 < h ≤ δ =⇒ | f (a − h) − f (a−)| ≤ ch Als f een eindige rechterafgeleide bezit in het punt a (eventueel nadat men in a de functiewaarde vervangen heeft door de rechterlimiet in a), dan voldoet f in a aan een rechtervoorwaarde van Lipschitz. Immers, de limiet f+0 (a) = lim
h→0+
f (a + h) − f (a+) h
bestaat en is eindig. Er bestaat dus een δ > 0 zodanig dat voor 0 < h < δ geldt: |
f (a + h) − f (a+) − f+0 (a)| < 1 h
of ( f+0 (a) − 1)h < f (a + h) − f (a+) < ( f+0 (a) + 1)h en hieruit volgt dat | f (a + h) − f (a+)| < (| f+0 (a)| + 1)h Zonder bewijs formuleren we volgende belangrijke stellingen. Stelling 2.5.5 (Stelling van Dirichlet) Onderstel dat de functie f : R → R stuksgewijs continu is, en periodiek is met periode 2π, en dat ze in elk punt voldoet aan een linker- en een rechtervoorwaarde van Lipschitz. In elk punt x ∈ R convergeert de Fourierreeks (2.18) van f naar f (x+) + f (x−) 2 In het bijzonder convergeert de Fourierreeks naar f (x) in de punten waar f continu is. De co¨effici¨enten an en bn van de Fourierreeks voldoen aan de gelijkheid van Parseval: ∞ a20 1 + ∑ (a2n + b2n ) = 2 n=1 π
Z 2π
f (x)2 dx
(2.25)
0
Stelling 2.5.6 Onderstel dat f : R → R overal continu is, en periodiek is met periode 2π. Als de afgeleide functie f 0 stuksgewijs continu is, dan convergeert de Fourierreeks van f uniform naar f over R.
34
Voorbeeld 2.5.7 Gevraagd wordt om de periodieke functie f : R → R gedefinieerd door als 0 < x < π 1 f (x) = 0 als x = 0, ±π −1 als −π < x < 0 te ontwikkelen in een Fourierreeks. f wordt geschetst in Figuur 2.4. blokfunctie
2 1.5 1
y-as
0.5 0 -0.5 -1 -1.5 -2 -10
-8
-6
-4
-2
0
2
4
6
8
10
x-as
Figuur 2.4: De blokfunctie Uit de stelling van Dirichlet kunnen we concluderen dat de Fourierreeks puntsgewijs naar de functie f zal convergeren. Aangezien f een oneven functie is, kunnen we onmiddellijk concluderen dat an = 0. De co¨effici¨enten bn berekenen we als volgt. bn =
2 π
Z π
f (x) sin nxdx " #π 2 cos nx = − π n 0
0
2 = − (cos nπ − cos 0) ( nπ 0 als n even 4 = als n oneven nπ Voor elke x ∈ R hebben we dus f (x) =
4 ∞ sin(2n + 1)x ∑ 2n + 1 π n=0
(2.26)
Enkele parti¨ele sommen van de reeks (2.26) worden geschetst in Figuur 2.5 In de omgeving van het discontinu¨ıteitspunt is de convergentie - zoals verwacht - trager dan elders. Men treft in de onmiddellijke omgeving van het discontinu¨ıteitspunt steeds een “piekafwijking”waar. Men noemt deze afwijking het Gibbsverschijnsel.
35
blokfunctie: eerste en tweede partiele som
1
1
0.5
0.5
0
0
-0.5
-0.5
-1
-1
-1.5 -4
-3
-2
-1
0
1
2
blokfunctie: derde en vierde partiele som
1.5
y-as
y-as
1.5
3
-1.5 -4
4
-3
-2
-1
x-as blokfunctie: vijfde en zesde partiele som
1
1
0.5
0.5
0
-0.5
-1
-1
-3
-2
-1
0
2
3
4
2
3
4
0
-0.5
-1.5 -4
1
1
blokfunctie: twaalfde partiele som
1.5
y-as
y-as
1.5
0 x-as
2
3
-1.5 -4
4
x-as
-3
-2
-1
0
1
x-as
Figuur 2.5: Parti¨ele sommen van de Fourierreeks van een blokfunctie Voorbeeld 2.5.8 Beschouw de functie f : [0, π] → R : x → x. Gevraagd wordt om deze te schrijven als een Fourierreeks met periode 2π. We kunnen dit op oneindig veel manieren doen: we kunnen de functie f willekeurig kiezen op [−π, 0], en dan de functie periodiek voortzetten. Er zijn twee mogelijke natuurlijke keuzes: we kunnen f op [−π, 0] zo kiezen dat f een oneven functie is, en we kunnen f op [0, π] zo bepalen dat f een even functie is. Grafisch ziet f er in de twee gevallen eruit als geschetst in Figuur 2.6. In het eerste geval is de Fourierreeks van f van de vorm ∞
∑ bn sin nx
n=1
( f is een oneven functie met periode 2π). De Fourierco¨effici¨enten van f kunnen gemakkelijk berekend worden (doe dit zelf). We vinden sin 2x sin 3x + −··· (2.27) f (x) = 2 sin x − 2 3 Enkele parti¨ele sommen worden geschetst in Figuur 2.7. Merk de afwijkende piek op in de buurt
36
zaagtand 1
3
3
2
2
1
1
0
0
-1
-1
-2
-2
-3
-3
-4 -10
-8
-6
-4
-2
0
zaagtand 2
4
y-as
y-as
4
2
4
6
8
-4 -10
10
-8
-6
-4
-2
x-as
0
2
4
6
8
10
x-as
Figuur 2.6: De zaagtandfunctie van de discontinu¨ıteiten. Dit is weer het Gibbsverschijnsel. Op analoge manier bepalen we de co¨effici¨enten van de tot een even functie voortgezette functie f . Dit keer vinden we een Fourierreeks van de vorm ∞ a0 + ∑ an cos nx 2 n=1
Bereken zelf de co¨effici¨enten ai . We vinden nu π 4 cos 3x cos 5x cos 7x f (x) = − cos x + 2 + 2 + 2 + · · · 2 π 3 5 7
(2.28)
Enkele parti¨ele sommen worden geschetst in Figuur 2.8 Merk op dat deze Fourierreeks “sneller”convergeert dan de vorige. Hier zijn twee verklaringen voor: in het tweede geval is de functie continu, en een functie die overal continu is, kan natuurlijk gemakkelijker benaderd worden door een goniometrische veelterm dan een discontinue. We zien ook dat de co¨effici¨enten in (2.28) sneller naar nul gaan dan die in (2.27), zodat we een snellere convergentie kunnen verwachten. Door in (2.27) of (2.28) specifieke waarden voor x te substitueren vinden we de som van een numerieke reeks. Als we in (2.27) x = π/2 stellen, dan vinden we π 1 1 1 1− + − +··· = 3 5 7 4 Als we in (2.28) x = π/4 stellen, dan vinden we √ 1 1 1 1 1 π 2 1 1− 2 − 2 + 2 + 2 − 2 − 2 ··· = 3 5 7 9 11 13 16 Ook de gelijkheid van Parseval levert interessante formules. Als we de gelijkheid van Parseval opschrijven voor de reeks (2.27) vinden we 1 π2 = ∑ 2 6 n=1 n ∞
37
zaagtand 1: eerste en tweede partiele som
3
3
2
2
1
1
0
0
-1
-1
-2
-2
-3
-3
-4 -4
-2
0
2
4
zaagtand 1: derde en vierde partiele som
4
y-as
y-as
4
6
8
-4 -4
10
-2
0
2
x-as zaagtand 1: vijfde en zesde partiele som
3
3
2
2
1
1
0
-1
-2
-2
-3
-3
-2
0
2
8
10
8
10
0
-1
-4 -4
6
4
zaagtamd 1: elfde en twaalfde partiele som
4
y-as
y-as
4
4 x-as
6
8
-4 -4
10
-2
0
x-as
2
4
6
x-as
Figuur 2.7: Parti¨ele sommen van de Fourierreeks van de eerste zaagtandfunctie
2.6
De Fourierintegraal
In § 2.5 hebben we de Fourierreeks bekeken van een periodische functie, of, equivalent, een functie gedefinieerd op een gesloten interval. We laten nu de periode naar oneindig gaan, en bekijken dus een functie gedefinieerd op gans de re¨ele rechte: f : R → R. We noemen F(p) = F { f }(p) =
Z +∞
f (t)e−ipt dt
(2.29)
−∞
de Fouriergetransformeerde van f . f kan terug berekend worden uit F via de formule 1 f (t) = 2π met andere woorden, 1 f (x) = 2π
Z +∞
F(p)eipt dp,
Z +∞ Z +∞
−ipt
f (t)e −∞
(2.30)
−∞
−∞
38
dt eipx dp.
(2.31)
zaagtand 2: eerste en tweede partiele som
3.5
3.5
3
3
2.5
2.5
2
2
1.5
1.5
1
1
0.5
0.5
0 -4
-2
0
2
4
zaagtand 2: derde en vierde partiele som
4
y-as
y-as
4
6
8
0 -4
10
-2
0
2
x-as
4
6
8
10
x-as
Figuur 2.8: Parti¨ele sommen van de Fourierreeks van de tweede zaagtandfunctie We zullen dit nu iets preciezer formuleren. Voor een functie f = u + iv : R → C defini¨eren we Z b
Z b
f (t)dt = a
Z b
u(t)dt + i
v(t)dt.
a
a
De oneigenlijke integraal wordt gedefinieerd door Z +∞
Z S
f (t)dt = lim
R→+∞ S→+∞
−∞
f (t)dt. −R
Indien deze limiet bestaat in C, dan noemen we de integraal convergent. De volgende eigenschap is een continue versie van de eigenschap dat een absoluut convergente reeks convergent is. Stelling 2.6.1 Neem f : R → C stuksgewijs continu. Als R ook 0+∞ f (t)dt.
R +∞ 0
| f (t)|dt convergeert, dan convergeert
Bewijs. We herinneren eraan dat q | f (t)| = u(t)2 + v(t)2 ; merk ook op |u(t)| ≤ | f (t)| en |v(t)| ≤ | f (t)|. Stel Z +∞
| f (t)|dt = I.
0
Dan hebben we, uit de definitie van limiet: ∀ε > 0, ∃T0 : ∀T > T0 : |I −
Z T 0
ε | f (t)|dt| < . 2
Dan volgt ook, voor alle S > T > T0 : Z S T
| f (t)|dt =
Z S
| f (t)|dt −
0
Z T 0
39
| f (t)|dt < ε.
We beschouwen nu de rij
Z n
un =
u(t)dt. 0
Voor m ≥ n > T0 hebben we |um − un | = |
Z m
u(t)dt| ≤
Z m
n
|u(t)|dt ≤
Z m
n
| f (t)|dt < ε,
n
en bijgevolg is (un ) een Cauchy rij, en een convergente rij. Stel lim un = I1 .
n→∞
Dan hebben we
ε ∀ε > 0, ∃N : ∀n > N : |un − I1 | < . 2 Als T > max{N + 1, T0 + 1}, dan is [T ] = n > N, en n, T > T0 , zodat |
Z T 0
u(t)dt − I1 | ≤ |
Z n 0
u(t)dt − I1 | + |
≤ |un − I1 | + ≤ |un − I1 | + We concluderen dat
Z T
Z T
u(t)dt| n
|u(t)|dt
n
Z T
| f (t)|dt <
n
ε ε + = ε. 2 2
Z +∞ 0
u(t)dt = I1
convergeert. Op dezelfde wijze kunnen we aantonen dat de integraal van −∞ tot 0 verloopt geheel analoog.
R +∞ 0
v(t)dt convergeert. Het bewijs voor
Aangezien |eipt | = 1 concluderen we: +∞ Gevolg 2.6.2 Neem f : R → C stuksgewijs continu. Als −∞ | f (t)|dt convergeert, dan convergeert R +∞ −ipt ook −∞ f (t)e dt. De Fouriergetransformeerde van f bestaat dus in elk punt p ∈ R.
R
Zonder bewijs formuleren we volgende stelling: Stelling 2.6.3 Onderstel dat f : R → C stuksgewijs continu is,R en in elk punt voldoet aan de +∞ linker- en rechtervoorwaarde van Lipschitz. Onderstel ook dat −∞ | f (t)|dt convergent is. Dan geldt formule (2.31) in elk punt x waarin f continu is. Als f discontinu is in x, dan is het rechterlid van (2.31) gelijk aan f (x+) + f (x−) . 2
40
Opmerkingen 2.6.4 1) Als f : R → R een re¨ele functie is, dan is F(p) = F(−p). 2) Als f : R → R een even re¨ele functie is (m.a.w. ( f (−x) = f (x)), dan is F(p) ∈ R, en Z 0
F(p) =
f (t)e−ipt dt +
Z +∞
−∞
0
Z +∞
=
f (t)e−ipt dt
f (u)eipu du +
Z +∞
0
f (t)e−ipt dt
0
Z +∞
= 2
f (t) cos ptdt. 0
3) Stel zelf een analoge formule op in het geval dat f : R → R een oneven re¨ele functie is. Voorbeeld 2.6.5 We berekenen de Fouriergetransformeerde van een signaal dat constant is gedurende het tijdsinterval [−R, R], en nul erbuiten: 1 als −R ≤ t ≤ R f1 (t) = 0 als |t| > R Z +∞
F1 (p) = 2
Z R
f (t) cos ptdt = 2 sin pt R sin pR = 2 =2 . p 0 p 0
cos ptdt 0
De sinc-functie wordt gegeven door de formule sinc(x) =
sin x . x
We krijgen dan F1 (p) = 2Rsinc(pR). Merk op: hoe groter R (dus hoe breder het signaal), hoe “smaller” de grafiek van het getransformeerde signaal: de afstand van de oorsprong tot het eerste nulpunt van F1 (p) is π/R. Als we de formule (2.30) toepassen in het geval t = 0 en R = 1, dan vinden we 1 1= 2π of
Z +∞ sin p
2
−∞
Z +∞ sin p −∞
p
p
dp,
dp = π.
Merk op dat de klassieke methodes (substitutie, parti¨ele integratie,...) niet toelaten om een primitieve functie van de sinc functie te bepalen.
41
1
0.8
0.6
0.4
0.2
0
!0.2
!0.4 !25
!20
!15
!10
!5
0
5
10
15
20
25
Figuur 2.9: De Fouriergetransformeerde van een rechthoekige puls (R = 1) Voorbeeld 2.6.6 In plaats van een rechthoekig signaal nemen we nu een driehoekig signaal, beschreven door de formule t 1 − R als 0 ≤ t ≤ R f2 (t) = 1 + Rt als −R ≤ t ≤ 0 0 als |t| > R We vinden nu, na parti¨ele integratie +∞ R t F2 (p) = 2 f (t) cos ptdt = 2 (1 − ) cos ptdt R 0 0 R Z R 2 t 2 = (1 − ) sin pt + sin(pt)dt p R pR 0 0 2(1 − cos pR) 2 = − 2 [cos(pt)]R0 = p R p2 R 4 sin2 (pR/2) pR = = Rsinc2 . 2 p R 2
Z
Z
Merk op dat F2 (0) = R. Het kleinste positieve nulpunt van F2 (p) is p = 2π R . Ook hier zien we het verschijnsel dat de breedte van het signaal omgekeerd evenredig is met de breedte van het getransformeerde signaal. We passen weer (2.30) toe, in het geval t = 0 en R = 2. We hebben dan F2 (p) = 2sinc2 p, en 1 1 = f (0) = 2π
Z +∞ −∞
42
2sinc2 pdp,
of
Z +∞ sin2 p
p2
−∞
dp = π.
(2.32)
1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 !25
!20
!15
!10
!5
0
5
10
15
20
25
Figuur 2.10: De Fouriergetransformeerde van een driehoekige puls (R = 1) We geven nu een aantal elementaire eigenschappen van de Fourierintegraal. In de volgende stellingen is f : R → C een functie die aan de voorwaarden van stelling 2.6.3 voldoet. We noteren F { f } = F. Stelling 2.6.7 Neem a ∈ R, en definieer fa : R → C door de formule fa (t) = f (t − a). De grafiek van fa is dan dezelfde als die van f , maar over een afstand a verschoven. We hebben dan F { fa }(p) = e−ipa F(p). (2.33) Bewijs.
F { fa }(p) =
Z +∞ −∞ Z +∞
=
f (t − a)e−ipt dt f (u)e−ipu e−ipa du = e−ipa F(p).
−∞
We voerden hierbij de substitutie u = t − a uit.
Stelling 2.6.8 Neem a > 0. Dan is 1 a
p a
F { f (at)}(p) = F( ).
43
(2.34)
Bewijs. Z +∞
F { f (at)}(p) =
f (at)e−ipt dt
−∞
ipu 1 +∞ = f (u)e− a du a −∞ 1 p F( ). = a a
Z
Ditmaal gebruikten we de substitutie u = at.
Als a > 1 in stelling 2.6.8, dan is het signaal f (at) “breder”dan het signaal f (t). Het getransformeerde signaal wordt dan smaller. Dit is het fenomeen dat we reeds opmerkten in de voorbeelden 2.6.5 en 2.6.6. Stelling 2.6.9
F {t n f (t)}(p) = in
dn F (p). dpn
(2.35)
Bewijs.
F {t n f (t)}(p) =
Z +∞
t n f (t)e−ipt dt
−∞ n
Z +∞
= i
−∞
∂n e−ipt f (t) dt ∂pn
dn F = in n (p). dp Stelling 2.6.10
F { f (n) (t)}(p) = (ip)n F(p).
(2.36)
Bewijs. 0
F { f (t)}(p) = =
Z +∞
f 0 (t)e−ipt dt
−∞
+∞ f (t)e−ipt −∞ + ip
Z +∞
f (t)e−ipt dt = ipF(p),
−∞
waarbij we parti¨ele integratie gebruikten. We passen deze formule n keer toe.
Beschouw twee functies f , g : R → C. Het convolutieproduct f ∗ g van f en g wordt gedefinieerd door de formule Z +∞ ( f ∗ g)(t) = f (t − x)g(x)dx. −∞
Merk op dat f ∗ g = g ∗ f : dit volgt als we de substitutie u = t − x uitvoeren. De Fouriergetransformeerde van het convolutieproduct is het product van de Fouriergetransformeerden.
44
Stelling 2.6.11
F { f ∗ g}(p) = F(p)G(p).
(2.37)
Bewijs.
F { f ∗ g}(p) =
Z +∞
Z +∞
dt −∞
−∞ +∞
Z +∞
=
dx −∞
dx −∞ Z +∞
=
Z
−∞ +∞
Z +∞
=
f (t − x)g(x)e−ipt dx
Z
f (t − x)e−ip(t−x) dt g(x)e−ipx f (s)e−ips ds g(x)e−ipx
−∞
dxF(p)g(x)e−ipx = F(p)G(p).
−∞
We verwisselden eerst de integratievolgorde, en voerden dan de substitutie s = t − x uit.
Beschouw weer twee functies f , g : R → C. Het inwendig product h f , gi wordt gedefinieerd door de formule Z +∞ h f , gi = (2.38) f (t)g(t)dt. −∞
Men kan aantonen dat de integraal (2.38) convergeert als f en g kwadratisch integreerbaar zijn, dit betekent dat Z +∞ Z +∞ | f (t)|2 dt, |g(t)|2 dt < +∞. −∞
−∞
De norm k f k van f wordt dan gedefinieerd als volgt: k f k2 = h f , f i =
Z +∞
| f (t)|2 dt.
−∞
Stelling 2.6.12 Onderstel f , G : R → C kwadratisch integreerbaar. hF { f }, Gi = 2πh f , F −1 {G}i.
(2.39)
Bewijs. 2πh f , F
−1
{G}i =
Z +∞
Z +∞
f (t)dt −∞ Z +∞
=
−∞ Z +∞
f (t)dt −∞
G(p)eipt dp G(p)e−ipt dp
−∞
Z +∞
=
F(p)G(p)dp. −∞
Een onmiddellijk gevolg van stelling 2.6.12 is het feit dat de Fouriertransformatie het inwendig product en de norm bewaren (op een factor 2π na).
45
Gevolg 2.6.13 (Formule van Plancherel) Onderstel f , g : R → C kwadratisch integreerbaar. hF, Gi = 2πh f , gi, kFk2 = 2πk f k2 . Bewijs.
(2.38)
hF { f }, F {g}i = 2πh f , F −1 {F {g}}i = 2πh f , gi. Voorbeeld 2.6.14 We passen de formule van Plancherel toe op voorbeelden 2.6.5 en 2.6.6. Neem eerst R = 1 in voorbeeld 2.6.5. We vinden formule (2.32) terug: 4
Z +∞ sin2 p −∞
p2
Z +∞
Z 1
2
dp = 2π
f1 (t) dt = 2π
−∞
dp = 4π. −1
Nu nemen we R = 2 in voorbeeld 2.6.6. We krijgen 4
Z +∞ sin4 p −∞
p4
Z +∞
dp = 2π −∞
Z 2
2
f2 (t) dt = 4π
zodat
Z +∞ sin4 p
p4
−∞
0
2 8π t 8π t 2 3 ( − 1) = . (1 − ) dt = 2 3 2 3 0
dp =
2π . 3
Uit de formule van Plancherel volgt ook dat hF1 , F2 i = 4
Z +∞ sin3 p
p3
−∞
dp
1 t = 2πh f1 , f2 i = 4π (1 − )dt 2 0 1 t2 = 4π t − = 3π, 4 0
Z
en hieruit volgt dat Z +∞ sin3 p
p3
−∞
dp =
3π . 4
We beschouwen tenslotte de constante functie g(t) =
1 2π
De Fouriergetransformeerde bestaat niet, omdat de integraal 1 π
Z ∞
costdt 0
divergeert. We onderstellen - tegen beter weten in - dat de Fouriergetransformeerde wel bestaat, en schrijven 1 F { }(p) = δ(p). 2π
46
We kunnen dan (2.38) toepassen. Neem f : R → C, and stel F = F { f }. We vinden dan Z +∞ −∞
1 F(p)δ(p)dp = hF { f }, δi = 2πh f , i = 2π
Z +∞
f (t)dt = F(0). −∞
δ heeft dus de volgende eigenschap, die men de zifteigenschap noemt: Z +∞
F(p)δ(p)dp = F(0).
(2.40)
−∞
Er is geen enkele functie δ die de zifteigenschap heeft; dit wekt natuurlijk geen verwondering, want de Fouriergetransformeerde die we bekijken bestaat niet. Wel bestaan er functies die in benadering de zifteigenschap hebben: voor elke R > 0 bekijken we de functie δR : R → C gegeven door ( 1 1 R als − 2R ≤ p ≤ 2R δR (p) = 1 0 als |p| > 2R Met behulp van de stelling van het gemiddelde berekenen we, voor elke (continue) functie F : R → C dat Z +∞ Z 1/2R F(p)δR (p)dp = R F(p)dp = F(ξ), −∞
−1/2R
waarbij ξ ∈ (−R/2, R/2). Louter formeel kunnen we dus stellen: δ(p) = lim δR (p). R→∞
Een meer correcte interpretatie is de volgende. Neem een deelruimte V van de vectorruimte bestaande uit alle absoluut integreerbare functies g : R → C. Een veralgemeende functie of distributie is per definitie een lineaire afbeelding ϕ : V → C. Aan een functie f : R → C kunnen we een distributie associ¨eren, namelijk ϕ f : V → C, ϕ f (g) =
Z +∞
f (x)g(x)dx. −∞
Maar niet elke distributie is afkomstig van een functie; bekijk bijvoorbeeld δ : V → C, δ(g) = g(0). Deze distributie noemt men de Dirac distributie. De formule voor de distributie geassocieerd aan een functie wordt uitgebreid naar willekeurige distributies: men schrijft Z +∞
ϕ(g) =
−∞
ϕ(x)g(x)dx.
47
Deel II Gewone en parti¨ele differentiaalvergelijkingen
48
Hoofdstuk 3 Differentiaalvergelijkingen 3.1
Definitie en voorbeelden
Een differentiaalvergelijking is een vergelijking waarin de onbekende een functie y van e´ e´ n veranderlijke x is, en waarin naast x en y ook de afgeleiden van de functie y naar x voorkomen. Het is dus een vergelijking van de vorm f (x, y, y0 , y00 , · · · , y(n) ) = 0
(3.1)
Het getal n, de hoogste afgeleide die in de vergelijking voorkomt, noemt men de orde van de differentiaalvergelijking. Een integraal of oplossing van de differentiaalvergelijking is een functie y = y(x) die aan (3.1) voldoet voor elke x gelegen in een open interval. De grafiek van een oplossing y = y(x) noemt men ook wel een integraalkromme van de differentiaalvergelijking. Een differentiaalvergelijking integreren of oplossen is er al de oplossingen van bepalen. Een normale differentiaalvergelijking is een differentiaalvergelijking van de vorm y(n) (x) = f (x, y, y0 , y00 , · · · , y(n−1) ) Voorbeelden 3.1.1 1) y0 = g(x), waarbij g een gegeven continue functie is. Als G een primitieve van g is, dan worden de oplossingen van y0 = g(x) gegeven door y = G(x) + c waarbij c een willekeurige constante is. 2) y0 = y. De oplossingen van deze vergelijking zijn y = cex waarbij alweer c een willekeurige constante is. 3) y00 = 0. Dan is y0 = c, en y = cx + d
49
waarbij c en d twee willekeurige constanten zijn. 4) y00 + y = 0. De oplossingen zijn nu y = c sin x + d cos x waarbij c en d twee willekeurige constanten zijn. Uit deze voorbeelden blijkt dat de oplossing van een differentiaalvergelijking niet uniek is. Voor onze voorbeelden van orde 1 vinden we oplossingen afhangende van e´ e´ n constante. Voor onze voorbeelden van orde 2 vinden we oplossingen afhangende van twee constanten. We voeren daarom volgende terminologie in. Beschouw een differentiaalvergelijking van n-de orde. Een algemene integraal is een oplossing van de differentiaalvergelijking die afhangt van n willekeurige constanten. Een particuliere integraal is een oplossing die men verkrijgt door aan de constanten in de algemene integraal bijzondere waarden te geven. Een singuliere integraal is een oplossing die men niet kan verkrijgen door aan de constanten in de algemene integraal bijzondere waarden te geven. Voorbeeld 3.1.2 y02 = x2 y. y=
x2
2 +c
4 is een algemene integraal van de vergelijking. Immers, y0 = 2
x2
2x x2 +c =x +c 4 4 4
en men heeft dus y0 (x)2 = x2 y(x), voor elke x ∈ R. De oplossing hangt af van de willekeurige constante c. y = x4 /16 is een particuliere integraal (stel c = 0). y = 0 is een singuliere integraal. Immers, het is een oplossing die niet kan verkregen door aan c een particuliere waarde te geven. Veel problemen uit de natuurkunde (en uit andere natuurwetenschappen en economische en sociale wetenschappen) geven aanleiding tot differentiaalvergelijkingen. We geven enkele zeer eenvoudige voorbeelden. Voorbeelden 3.1.3 1) Een deeltje met massa m is bevestigd aan een veer. Als men de veer uitrekt over een afstand x, dan zal de veer een terugroepende kracht uitoefenen op het deeltje die in eerste benadering rechtevenredig is met x: F = −kx waarbij k > 0 een constante is. Kracht is massa maal versnelling, en dus geldt m
d2 x = −kx dt 2
50
m
u
x
0 F
Figuur 3.1: een veer De algemene integraal van deze differentiaalvergelijking is r r k k x = c sin t + d cos t m m 2) In de buurt van het aardoppervlak is de zwaartekracht die werkt op een deeltje in goede benadering rechtevenredig met de massa van dit deeltje: F = −gm De constante g (de valversnelling) kan experimenteel bepaald worden en is in benadering gelijk aan 9.81 m/sec2 . Als x de hoogte voorstelt, dan wordt de valbeweging van het deeltje dus beschreven door de differentiaalvergelijking d2 x m 2 = −mg dt of d2 x = −g dt 2 De algemene integraal van deze differentiaalvergelijking is x = −g
t2 + ct + d 2
3) N(t) stelt een populatie op tijdstip t voor (bijvoorbeeld het aantal bijen in een bijenkorve, of het aantal inwoners van de stad Antwerpen). Als we onderstellen dat de aangroei van de populatie (de afgeleide van N naar de tijd) rechtevenredig is met de populatie zelf, dan vinden we dat N een oplossing is van de differentiaalvergelijking N 0 (t) = λN(t) waarbij λ een constante is. De oplossing van deze vergelijking is N(t) = ceλt waarbij c een constante is. Als we t = 0 invullen vinden we c = N(0), en dus N(t) = N(0)eλt We spreken van een exponenti¨ele groei.
51
Differentiaalvergelijkingen hebben ook meetkundige toepassingen. Beschouw een familie vlakke krommen die afhangt van twee constanten: ϕ(x, y, c, d) = 0
(3.2)
Tweemaal afleiden naar x geeft dϕ ∂ϕ ∂ϕ 0 = + y =0 dx ∂x ∂y d2 ϕ ∂2 ϕ ∂2 ϕ 0 ∂2 ϕ 0 2 ∂ϕ 00 = + 2 y + 2 (y ) + y = 0 dx2 ∂x2 ∂x∂y ∂y ∂y
(3.3) (3.4)
Eliminatie van c en d uit (3.2,3.3,3.4) geeft een betrekking tussen x, y, y0 en y00 , een differentiaalvergelijking van orde 2. f (x, y, y0 , y00 ) = 0 (3.5) De krommen (3.2) zijn allen integraalkrommen van (3.5). Men noemt (3.5) de differentiaalvergelijking van de familie krommen (3.2). Bepaalde meetkundige vraagstukken kunnen aanleiding geven tot differentiaalvergelijkingen. Voorbeeld 3.1.4 Neem een vlakke kromme met vergelijking y = y(x). Zij M een punt op de ~ k= kromme en T het snijpunt van de raaklijn door M en de x-as. Voor welke krommen geldt kMT ~ kOT k voor elke M op de kromme? y
M(x, y) x (x1, y1)
T
0
x
Figuur 3.2: Een meetkundig vraagstuk De vergelijking van de raaklijn in M(x, y) is Y − y = y0 (x)(X − x) Als Y = 0, dan is X = x −
y y ~ k2 wordt dus ~ k2 = kOT , en dus is T (x − 0 , 0). De voorwaarde kMT 0 y y y 2 y2 x − 0 = 02 + y2 y y
of y0 =
2xy x2 − y2
52
3.2
Beginvoorwaarden en randvoorwaarden
Zoals we gezien hebben hangt de algemene integraal van een differentiaalvergelijking van n-de orde af van n constanten. Door n bijkomende voorwaarden op te leggen kan men deze n constanten vastleggen. De oplossing wordt dan uniek. Het vraagstuk met beginvoorwaarden en de bestaansstelling Men zoekt een oplossing van de differentiaalvergelijking (3.1) die zodanig is dat voor een gegeven waarde x0 : y(x0 ) = y0 y0 (x0 ) = y00 .. . (n−1) (n−1) y (x0 ) = y0 (n−1)
waarbij y0 , y00 , · · · , y0
n gegeven getallen zijn.
Voorbeeld 3.2.1 Op tijdstip t = 0 laat men een steen vallen vanaf hoogte h. Op welk tijdstip t zal de steen de grond bereiken? We hebben hierboven gezien dat de valbeweging van beschreven wordt door de differentiaalvergelijking d2 x = −g dt De algemene integraal van deze differentiaalvergelijking is x = −g
t2 + ct + d 2
(3.6)
Op tijdstip t = 0 is x(0) = h en
dx (0) = 0 dt
(3.7)
Als we (3.7) invullen in (3.6) dan vinden we dat d = h en c = 0 zodat x(t) = h − g We vinden dat x = 0 als t =
t2 2
p 2h/g.
Voor een normale differentiaalvergelijking heeft het beginvoorwaardeprobleem onder bepaalde voorwaarden altijd een unieke oplossing. De onderliggende stelling is de volgende:
53
(n−1)
Stelling 3.2.2 (Bestaansstelling) Onderstel dat G ⊂ Rn+1 een gebied is, en dat (x0 , y0 , y00 , y000 , · · · , y0 een inwendig punt van G is. Als f : G ⊂ Rn+1 → R continu is over G en differentieerbaar over het inwendige van G, dan bestaat er een omgeving V = (x0 − c, x0 + c) van x0 en een unieke functie g : V → R zodat • g is een oplossing van de differentiaalvergelijking y(n) (x) = f (x, y, y0 , y00 , · · · , y(n−1) ) voor elke x ∈ V ; • g voldoet aan de beginvoorwaarden: g(x0 ) = y0 0 0 g00(x0 ) = y000 g (x0 ) = y0 . .. (n−1) (n−1) g (x0 ) = y0 Merk op dat de stelling ons vertelt dat er een oplossing bestaat, maar niet hoe die oplossing expliciet kan geconstrueerd worden. Zulk een situatie zijn we al eens tegengekomen: de stelling van de inverse functie, en de stelling van de impliciete functie vertelden ons onder welke voorwaarde een algebraische vergelijking een unieke oplossing heeft, maar vertelde er niet bij hoe die oplossing kon geconstrueerd worden. Hetzelfde zien we hier, maar dan voor een nieuw type vergelijkingen, namelijk differentiaalvergelijkingen. Het is dan ook niet toevallig dat de bewijzen van beide stellingen gebruik maken van dezelfde wiskundige technieken. We zien hier een fenomeen dat in de moderne wiskunde meermaals opduikt: om de stelling van de impliciete functie en de bestaansstelling voor differentiaalvergelijkingen te kunnen bewijzen moeten we eerst een abstracte theorie ontwikkelen: men voert metrische ruimten, dit zijn verzamelingen waarop een afstand gedefinieerd is. Vervolgens bewijst men de dekpuntstelling, deze geeft nodige en voldoende voorwaarden opdat een functie van een metrische ruimte naar zichzelf een uniek dekpunt heeft. Tot slot merken we op dat de dekpuntstelling nog andere toepassingen heeft: ze kan gebruikt worden om de convergentie van bepaalde algoritmen om vergelijkingen numeriek op te lossen te bespreken. Verdere details laten we hier achterwege. Voorbeeld 3.2.3 In de meeste gevallen die men in de praktijk tegenkomt is aan de voorwaarden van de bestaansstelling voldaan. Toch zijn er voorbeelden van differentiaalvergelijkingen met een beginvoorwaarde waarvoor de oplossing niet uniek is. Neem bijvoorbeeld de vergelijking √ y0 = 2 y De functie
y=
(x − c)2
als x ≥ c;
0 als x < c. is een oplossing voor iedere c. De differentiaalvergelijking met de beginvoorwaarde y(x0 ) = 0 heeft een oneindig aantal oplossingen. Waarom is de bestaansstelling hier niet van toepassing?
54
)
Hoofdstuk 4 Differentiaalvergelijkingen van eerste orde In dit hoofdstuk bestuderen we normale differentiaalvergelijkingen van orde 1: y0 = g(x, y)
(4.1)
In een aantal gevallen kunnen we de algemene integraal van (4.1) expliciet uitrekenen. Merk op dat we (4.1) altijd kunnen herschrijven onder de vorm p(x, y)dx + q(x, y)dy = 0
4.1
(4.2)
Juiste differentiaalvergelijkingen
Totale differentialen Beschouw een differentieerbare functie van twee veranderlijken z = g(x, y) De differentiaal van g wordt gegeven door dz =
∂g ∂g dx + dy ∂x ∂y
Onderstel nu, omgekeerd dat een uitdrukking van de vorm p(x, y)dx + q(x, y)dy gegeven is. Onder welke voorwaarden is deze uitdrukking de totale differentiaal van een functie? Of, equivalent, wanneer is p(x, y)~u1 + q(x, y)~u2 de gradi¨ent van een functie z = g(x, y)?
55
Stelling 4.1.1 Onderstel dat p(x, y) en q(x, y) continue parti¨ele afgeleiden van eerste orde bezitten in een enkelvoudig samenhangend domein D. Dan is p(x, y)dx + q(x, y)dy een totale differentiaal over D als en alleen als ∂p ∂q = ∂y ∂x over D. Bewijs. Onderstel dat er een functie z : D → R bestaat zodat dz = p(x, y)dx + q(x, y)dy voor elke (x, y) ∈ D. Dan is ∂z ∂z p= en q = ∂x ∂y zodat ∂2 z ∂2 z ∂q ∂p = = = ∂y ∂y∂x ∂x∂y ∂x Omgekeerd, onderstel dat ∂p ∂q = ∂y ∂x over D. We zullen de stelling bewijzen in het geval waarin D het inwendige van een rechthoek is. Kies (x0 , y0 ) ∈ D willekeurig. We zoeken een functie z = g(x, y) zodat p=
∂z ∂z en q = ∂x ∂y
Uit de eerste vergelijking volgt dat Z x
z=
p(t, y)dt + ϕ(y) x0
waarbij ϕ een willekeurige functie in 1 veranderlijke is. Als we dit in de tweede vergelijking stoppen, dan vinden we, gebruik makend van de regel van Leibniz q(x, y) = = = = =
∂z ∂y Z ∂ x p(t, y)dt + ϕ(y) ∂y x0 Z x ∂p (t, y)dt + ϕ0 (y) x0 ∂y Z x ∂q (t, y)dt + ϕ0 (y) x0 ∂x q(x, y) − q(x0 , y) + ϕ0 (y)
zodat ϕ0 (y) = q(x0 , y) en
Z y
ϕ(y) = y0
q(x0 , u)du + c
56
en de gezochte functie g is dus Z x
g(x, y) =
Z y
p(t, y)dt + x0
y0
q(x0 , u)du + c
(4.3)
Door de rollen van x en y te verwisselen kunnen we z = g(x, y) ook schrijven onder de vorm Z y
g(x, y) =
Z x
q(x, u)du + y0
x0
p(t, y0 )dt + c
(4.4)
De constante c is in beide gevallen niets anders dan g(x0 , y0 ). Beide vormen zijn gelijkwaardig want twee functies die dezelfde totale differentiaal hebben in D zijn gelijk op een constante na. Waar in het bewijs gebruikten we dat D het inwendige van een rechthoek is? Voorbeeld 4.1.2 Bepaal z = g(x, y) zo dat grad g = (3x2 − 1)y~u1 + (x3 − x + 2y)~u2 of dg = (3x2 − 1)ydx + (x3 − x + 2y)dy Stel p(x, y) = (3x2 − 1)y en q(x, y) = x3 − x + 2y. We zien gemakkelijk dat ∂q ∂p = 3x2 − 1 = ∂y ∂x zodat aan de voorwaarden van stelling 4.1.1 voldaan is. Om g te bepalen kan men (4.3) of (4.4) rechtstreeks toepassen. Men kan ook de redenering uit bovenstaande bewijsvoering herhalen, en dit is wat meestal in de praktijk gebeurt. Omdat ∂z = (3x2 − 1)y ∂x is
Z
z=
(3x2 − 1)ydx + ϕ(y) = (x3 − x)y + ϕ(y)
Als we dit stoppen in de tweede vergelijking ∂z = x3 − x + 2y ∂y dan volgt x3 − x + ϕ0 (y) = x3 − x + 2y zodat ϕ0 (y) = 2y en ϕ(y) = y2 + c Tenslotte is z = (x3 − x)y + y2 + c
57
Juiste differentiaalvergelijkingen De differentiaalvergelijking p(x, y)dx + q(x, y)dy = 0 wordt een juiste differentiaalvergelijking genoemd indien p(x, y)dx + q(x, y)dy een totale differentiaal is. Indien dg = p(x, y)dx + q(x, y)dy, dan wordt de oplossing van de differentiaalvergelijking in impliciete vorm gegeven door de formule g(x, y) = c Voorbeelden 4.1.3 1) De vergelijking x ln ydx + dy = 0 y is een juiste differentiaalvergelijking. Immers x d(xln y) = ln ydx + dy y en de oplossing van de differentiaalvergelijking is dus xln y = c of y = ec/x 2) We beschouwen de vergelijking x − yy0 x2 + y2 = x + yy0 a2 Als we deze vergelijking oplossen naar y0 , dan vinden we x(x2 + y2 − a2 )dx + y(x2 + y2 + a2 )dy = 0 Het is gemakkelijk in te zien dat dit een onmiddellijk integreerbare vergelijking is. Immers, ∂ ∂ (x(x2 + y2 − a2 )) = (y(x2 + y2 + a2 )) = 2xy ∂y ∂x We zoeken een functie z = g(x, y) zodat dz = x(x2 + y2 − a2 )dx + y(x2 + y2 + a2 )dy Uit ∂z = x(x2 + y2 − a2 ) ∂x volgt dat z=
x4 x2 y2 x2 a2 + − + ϕ(y) 4 2 2
58
Hieruit vinden we dat ∂z = x2 y + ϕ0 (y) = y(x2 + y2 + a2 ) ∂y zodat ϕ0 (y) = y3 + ya2 en ϕ(y) =
y4 y2 a2 + +c 4 2
Tenslotte is x4 x2 y2 x2 a2 y4 y2 a2 + − + + +c 4 2 2 4 2 1 2 a2 = (x + y2 )2 + (y2 − x2 ) + c 4 2
z =
en de oplossing van de differentiaalvergelijking is (x2 + y2 )2 + 2a2 (y2 − x2 ) = d Differentiaalvergelijkingen met gescheiden veranderlijken Onderstel dat we de differentiaalvergelijking y0 = f (x, y) kunnen herschrijven onder de vorm p(x)dx = q(y)dy Deze vergelijking is onmiddellijk integreerbaar, en de algemene integraal is Z
Z
p(x)dx =
q(y)dy
Voorbeelden 4.1.4 1) y0 = y of
dy = dx y
of ln y = x + ln c of y = cex Merk op : a) door te delen door y hebben we de oplossing y = 0 uitgesloten. Deze wordt echter opnieuw ingevoerd door finaal ook c = 0 toe te laten. b) Als y negatief is, dan moeten we schrijven ln |y| = x + ln c, of −y = cex , of y = −cex . Deze oplossing wordt opnieuw ingevoerd door finaal ook c < 0 toe te laten.
59
2) y0 = − of
y x
dy dx =− y x
of ln y = −ln x + ln c of xy = c Opmerking 4.1.5 Om de vergelijking p(x, y)dx + q(x, y)dy = 0 op te lossen tracht men soms een functie µ(x, y) te vinden zodanig dat µ(x, y)p(x, y)dx + µ(x, y)q(x, y)dy = 0 een juiste differentiaalvergelijking is. Men integreert dan deze vergelijking en gaat na of er soms geen vreemde oplossingen ingevoerd werden (de oplossingen van µ(x, y) = 0) of weggelaten werden (de oplossingen van 1/µ(x, y) = 0). Het probleem is op deze manier herleid tot het vinden van µ, de integrerende factor. µ moet voldoen aan de betrekking ∂(µp) ∂(µq) = ∂y ∂x Dit is een parti¨ele differentiaalvergelijking die moeilijker op te lossen is dan de oorspronkelijke vergelijking.
4.2
Homogene differentiaalvergelijkingen
Herhaal dat een veelterm P(x, y) homogeen van graad m genoemd wordt, indien de som van de exponenten van x en y in elk van de termen van de veelterm juist m is, m.a.w. m
P(x, y) = ∑ ai xi ym−i i=0
Als p(x, y) en q(x, y) homogene veeltermen van graad m zijn, dan noemen we de differentiaalvergelijking p(x, y)dx + q(x, y)dy = 0 een homogene differentiaalvergelijking van graad m. We herschrijven de differentiaalvergelijking onder de vorm p(x, y) dx dy = q(x, y)
60
Deel teller en noemer door xm . We krijgen dan dy =
y p(1, y/x) dx = f dx q(1, y/x) x
We kiezen nu een nieuwe veranderlijke z=
y x
Dan is y = zx en dy = zdx + xdz en de vergelijking wordt zdx + xdz = f (z)dx of
dz dx = f (z) − z x
en dit is een differentiaalvergelijking met gescheiden veranderlijken. Voor de algmene integraal kunnen we volgende formule opschrijven. Z dz x = c exp f (z) − z Voorbeeld 4.2.1 (x2 − y2 )dy = 2xydx Dit is een homogene differentiaalvergelijking van graad 2. Stel y = zx, dan wordt de vergelijking (1 − z2 )(zdx + xdz) = 2zdx of
dx 1 − z2 dz = z + z3 x
of 1 − z2 dz z + z3 Z Z dz 2z = dz − z 1 + z2 = ln z − ln (1 + z2 )
ln x − ln c =
Z
Hierbij splitsten we 1 − z2 /z + z3 in parti¨ele breuken. Door de exponent van beide leden te nemen vinden we cz cyx x= = 1 + z2 x2 + y2 of x2 + y2 − cy = 0 De integraalkrommen zijn dus cirkels die de x-as raken in de oorsprong.
61
Vergelijkingen die te herleiden zijn tot homogene vergelijkingen We bekijken differentiaalvergelijkingen van de vorm y0 = f
ax + by + c a0 x + b0 y + c0
Deze vergelijking is te herleiden tot een homogene vergelijking. We beschouwen de rechten l : ax + by + c = 0 en l 0 : a0 x + b0 y + c0 = 0 eerste geval: l en l 0 snijden elkaar in het punt (x0 , y0 ). We voeren de volgende co¨ordinatentransformatie uit: x = t + x0 y = u + y0 We zien gemakkelijk in dat ax + by + c = ax + by + c − (ax0 + by0 + c) = at + bu en a0 x + b0 y + c0 = a0t + b0 u en de differentiaalvergelijking wordt at + bu du =f 0 dt a t + b0 u en deze vergelijking is homogeen. tweede geval: l en l 0 zijn evenwijdig. Dan is a0 b0 = =λ a b en de vergelijking is van de vorm y0 = f
ax + by + c λ(ax + by) + c0
Stel u = ax + by, dan wordt de vergelijking u+c u0 − a =f b λu + c0 of
du = dx u+c a + b f λu+c 0
en dit is een differentiaalvergelijking met gescheiden veranderlijken.
62
Voorbeeld 4.2.2 Beschouw de differentiaalvergelijking y0 =
x−y−4 x+y−2
Het snijpunt van de rechten x − y = 4 en x + y = 2 is (3, −1). We stellen daarom x = t + 3 en y = u − 1. De differentiaalvergelijking wordt nu du t − u = dt t +u Zoals verwacht is deze homogeen van graad 1. Stel u = tz. We krijgen dan tdz + zdt = of
1−z dt 1+z
dt 1+z 1 2(z + 1) = dz dz = − 2 t 1 − 2z − z 2 (z + 1)2 − 2
of
1 1 lnt = − ln (z + 1)2 − 2 + ln c 2 2
of t 2 (z + 1)2 − 2 = c of u2 + 2ut − t 2 = c of (y + 1)2 + 2(y + 1)(x − 3) − (x − 3)2 = c of, tenslotte y2 + 2xy − x2 − 4y + 8x − 14 = c
4.3
Lineaire differentiaalvergelijkingen van eerste orde
Een lineaire differentiaalvergelijking van orde 1 is een differentiaalvergelijking van de vorm a(x)y0 + b(x)y = d(x)
(4.5)
De vergelijking is dus een veelterm van graad 1 in y en y0 . Als d(x) = 0, dan zeggen we dat (4.5) een homogene differentiaalvergelijking is. De differentiaalvergelijking a(x)y0 + b(x)y = 0
(4.6)
die ontstaat door in (4.5) d(x) gelijk aan 0 te stellen noemen we de geassocieerde homogene differentiaalvergelijking. In de volgende stelling zullen we zien dat er een verband bestaat tussen de oplossing van een lineaire differentiaalvergelijking en de oplossing van de geassocieerde homogene vergelijking.
63
Stelling 4.3.1 Als yh de algemene integraal is van de homogene vergelijking (4.6), en y p een particuliere integraal is van de volledige vergelijking (4.5), dan wordt de algemene integraal van de volledige vergelijking (4.5) gegeven door de formule y = y p + yh
(4.7)
Bewijs. We tonen eerst aan dat (4.7) een oplossing is van (4.5): a(x)(y p + yh )0 + b(x)(y p + yh ) = a(x)y0p + b(x)y p + a(x)y0h + b(x)yh = d(x) + 0 = d(x) Omgekeerd, onderstel dat y een oplossing is van (4.5). Dan is a(x)(y − y p )0 + b(x)(y − y p ) = a(x)y0 + b(x)y − a(x)y0p − b(x)y p = d(x) − d(x) = 0 zodat y − y p een oplossing is van de homogene vergelijking (4.6). y = y p + (y − y p ) is dus van de gewenste vorm. Om een lineaire differentiaalvergelijking te integreren volstaat het dus om eerst de geassocieerde homogene vergelijking te integreren en daarna een particuliere integraal van de volledige vergelijking te vinden. Integratie van de homogene vergelijking Dit is eenvoudig: de homogene vergelijking (4.6) kan herschreven worden als een vergelijking met gescheiden veranderlijken. a(x)dy + b(x)ydx = 0 of
dy b(x) =− dx y a(x)
Integratie levert ln y = −
b(x) dx + ln c a(x)
Z
of −
R b(x)
y = ce
a(x)
dx
(4.8)
Bepalen van een particuliere oplossing : variatie van de constante Om een particuliere oplossing van de volledige vergelijking te bepalen gebruiken we een techniek die variatie van de constante genoemd wordt. Uit (4.8) weten we dat de integraal van de homogene vergelijking van de vorm yh = c f (x) is. We zoeken een particuliere integraal van (4.5) van de vorm y p = c(x) f (x)
64
(4.9)
Substitutie in (4.5) levert a(x)c0 (x) f (x) + a(x)c(x) f 0 (x) + b(x)c(x) f (x) = d(x) of, omdat a(x) f 0 (x) + b(x) f (x) = 0 ( f is een oplossing van de homogene vergelijking): a(x)c0 (x) f (x) = d(x) Hieruit vinden we gemakkelijk dat d(x) a(x) f (x)
c0 (x) = en
d(x) dx a(x) f (x) Als we (4.10) substitueren in (4.9), dan vinden we een particuliere integraal van (4.5). Z
c(x) =
(4.10)
Voorbeeld 4.3.2 Beschouw de differentiaalvergelijking sin xy0 − y = 3
(4.11)
We integreren eerst de geassocieerde homogene vergelijking sin xy0 − y = 0 We vinden
dy dx = y sin x
of
x dx = ln tg + ln c sin x 2
Z
ln y = of
x 2 We zoeken nu een particuliere oplossing door variatie van de constante: yh = ctg
y p = c(x)tg Substitutie in (4.11) levert sin xc0 (x)tg zodat c0 (x) = en
y = tg
x =3 2
3 tg
3 c(x) = 2
Tenslotte is
x 2
Z
x x x 2 2 sin 2 cos 2
=
3 2 sin2 2x
dx x = −3cotg + c 2x 2 sin 2
x x x −3cotg + c = −3 + ctg 2 2 2
65
De vergelijking van Bernoulli Een differentiaalvergelijking van de vorm a(x)y0 + b(x)y = d(x)ym
(4.12)
noemen we een vergelijking van Bernoulli. Als m = 0 of m = 1, dan is de vergelijking lineair. Voor andere waarden van m kunnen we de vergelijking van Bernoulli tot een lineaire vergelijking herleiden bij middel van de substitutie z= Immers,
1 ym−1
y0 z = (1 − m) m y 0
en (4.12)wordt a(x)
z0 + b(x)z = d(x) 1−m
(4.13)
en dit is een lineaire vergelijking. Voorbeeld 4.3.3 We beschouwen de differentiaalvergelijking xy0 + y = y2 ln x Dit is een vergelijking van Bernoulli met m = 2. Substitueer z = 1/y. Dan is z0 = −y0 /y2 en de differentiaalvergelijking wordt xz0 − z = −ln x Dit is een lineaire vergelijking. We integreren eerst de homogene vergelijking xz0 − z = 0. Dit geeft zh = cx Stel nu (variatie van de constante) z p = c(x)x Subsitutie in de vergelijking levert x2 c0 (x) = −ln x en c(x) = −
Z
ln xdx ln x = − x2 x
Z
dx ln x 1 = + +c x2 x x
We hebben dus tenslotte dat z = 1 + ln x + cx of y=
1 1 + ln x + cx
66
Hoofdstuk 5 Normale differentiaalstelsels en vergelijkingen 5.1
De methode van afleiding en eliminatie
Een normale differentiaalvergelijking van orde n is een differentiaalvergelijking van de vorm y(n) = f (x, y, y0 , y00 , · · · , y(n−1) )
(5.1)
Merk op dat een normale differentiaalvergelijking van orde 1 niets anders is dan een differentiaalvergelijking van eerste orde en eerste graad zoals besproken in het voorgaande hoofdstuk. Een stelsel differentiaalvergelijkingen van orde n is een stel vergelijkingen van de vorm y01 = f1 (x, y1 , y2 , · · · , yn ) y0 = f2 (x, y1 , y2 , · · · , yn ) 2 (5.2) .. . y0 = f (x, y , y , · · · , y ) n
n
1
2
n
waarbij de onbekenden y1 , y2 , · · · , yn functies zijn van de veranderlijke x. Merk op dat een normaal differentiaalstelsel van orde 1 niets anders is dan een differentiaalvergelijking van eerste orde en eerste graad. We kunnen (5.2) herschrijven in vectorvorm: schrijf Y = (y1 , y2 , · · · , yn ) en F = ( f1 , f2 , · · · , fn ) Dan wordt (5.2) Y 0 = F(x,Y )
(5.3)
We zullen nu zien dat er een verband is tussen normale differentiaalvergelijkingen en differentiaalstelsels: als we normale differentiaalvergelijkingen kunnen oplossen, dan kunnen we ook differentiaalstelsels oplossen, en omgekeerd. Deze eigenschap zullen we in het vervolg in beide richtingen gebruiken.
67
Afleiding en eliminatie We vertrekken van het differentiaalstelsel (5.2). Afleiden van de eerste vergelijking geeft y001 =
n ∂ f1 ∂ f1 0 +∑ y ∂x i=1 ∂yi i
Als we hierin de laatste n − 1 vergelijkingen van (5.2) invullen, dan vinden we y001 (x) = g2 (x, y1 , y01 , y2 , · · · , yn ) Nogmaals afleiden geeft, na invullen van de laatste n − 1 vergelijkingen van (5.2): 0 00 y000 1 (x) = g3 (x, y1 , y1 , y1 , y2 , · · · , yn )
Verder afleiden geeft, samengevat 0 y1 = f1 (x, y1 , y2 , · · · , yn ) 00 (x) = g (x, y , y0 , y , · · · , y ) n 2 1 1 2 1 y000 y1 (x) = g3 (x, y1 , y01 , y001 , y2 , · · · , yn ) .. . (n) (n−1) y1 = gn (x, y1 , y01 , y001 , · · · , y1 , y2 , · · · , yn )
(5.4)
Als we uit (5.4) y2 , y3 , · · · , yn elimineren, dan vinden we een betrekking van de vorm (n)
(n−1)
y1 = g(x, y1 , y01 , y001 , · · · , y1
)
(5.5)
Onderstel dat we (5.5) kunnen oplossen. Dan kennen we y1 . Bij de eliminatie van y2 , y3 , · · · , yn (n−1) uit (5.4) berekenden we y2 , y3 , · · · , yn in functie van y1 , y01 , · · · , y1 , zodat y2 , y3 , · · · , yn kunnen berekend worden. We kunnen dus het differentiaalstelsel (5.2) oplossen als we de differentiaalvergelijking (5.5) kunnen oplossen. Deze methode om een differentiaalstelsel op te lossen noemt men de methode van afleiding en eliminatie. Voorbeeld 5.1.1 Beschouw het differentiaalstelsel y01 = y2 − 2y1 y02 = 2y2 − 3y1
(5.6) (5.7)
Afleiden van (5.6) geeft, gebruik makend van (5.7) y001 = y02 − 2y01 = 2y2 − 3y1 − 2y01 Uit (5.6) volgt y2 = y01 + 2y1 en dus is y001 = 2y01 + 4y1 − 3y1 − 2y01 = y1 en y1 = Aex + Be−x Substitutie in (5.8) geeft y2 = Aex − Be−x + 2Aex + 2Be−x = 3Aex + Be−x
68
(5.8)
Via afleiding en eliminatie kan men een normaal differentiaalstelsel herleiden tot een normale differentiaalvergelijking. Omgekeerd kan men een normale differentiaalvergelijking herleiden tot een normaal differentiaalstelsel. Bekijk de vergelijking y(n) = f (x, y, y0 , y00 , · · · , y(n−1) )
(5.9)
Stel y1 = y, y2 = y0 , · · · , yi = y(i−1) , · · · , yn = y(n−1) Als y een oplossing is van (5.9), dan geldt y01 = y2 0 y2 = y3 .. . y0n−1 = yn 0 yn = y(n) = f (x, y1 , y2 , y3 , · · · , yn )
(5.10)
en Y = (y1 , y2 , · · · , yn ) is dus een oplossing van het normale differentiaalstelsel (5.10). Als we de methode van afleiding en eliminatie toepassen op (5.10), dan vinden we (5.9) terug. Voor een normale differentiaalvergelijking van orde n heeft de bestaansstelling (onder bepaalde voorwaarden) een unieke oplossing. Een analoog resultaat geldt voor normale differentiaalstelsels. Stelling 5.1.2 (Bestaansstelling) Onderstel dat G ⊂ Rn+1 een gebied is, en dat (x0 ,Y0 ) een inwendig punt van G is. Als F : G ⊂ Rn+1 → Rn continu is over G en differentieerbaar over het inwendige van G , dan bestaat er een omgeving V = (x0 − c, x0 + c) van x0 en een unieke functie G : V → Rn zodat • G is een oplossing van het differentiaalstelsel Y 0 (x) = F(x,Y (x)) voor elke x ∈ V ; • G voldoet aan de beginvoorwaarde: G(x0 ) = Y0
5.2
Lineaire differentiaalstelsels
Een lineair differentiaalstelsel is een differentiaalstelsel van de vorm y01 = a11 (x)y1 + a12 (x)y2 + · · · + a1n (x)yn + b1 (x) y0 = a21 (x)y1 + a22 (x)y2 + · · · + a2n (x)yn + b2 (x) 2 .. . y0 = a (x)y + a (x)y + · · · + a (x)y + b (x) n
n1
1
2
n2
69
nn
n
n
(5.11)
Indien de functies b1 , b2 , · · · , bn identiek gelijk zijn aan 0, dan noemen we het stelsel (5.11) een homogeen lineair differentiaalstelsel. We kunnen (5.11) herschrijven in matrixvorm. Schrijf y1 (x) b1 (x) a11 (x) a12 (x) · · · a1n (x) y2 (x) b2 (x) a21 (x) a22 (x) · · · a2n (x) Y (x) = . , B(x) = . , A(x) = . .. .. . . . . . . . . yn (x)
an1 (x) an2 (x) · · · ann (x)
bn (x)
Dan is (5.11) equivalent met Y 0 (x) = A(x)Y (x) + B(x)
(5.12)
De bestaansstelling is steeds van toepassing op een lineair differentiaalstelsel. Stelling 5.2.1 Onderstel dat A : (a, b) → Mnn (R) en B : (a, b) → Rn differentieerbaar zijn over (a, b). Voor x0 ∈ (a, b) en Y0 ∈ Rn bestaat er een unieke oplossing van (5.12) gedefinieerd op een omgeving van x0 die voldoet aan de beginvoorwaarde Y (x0 ) = Y0 Bewijs. Dit volgt onmiddellijk uit de bestaansstelling: de functie F(x,Y ) = A(x)Y +B(x) is continu over een gebied G waarvan (x0 ,Y0 ) een inwendig punt is, en differentieerbaar over het inwendige van G . Eigenschappen van het homogeen stelsel We beschouwen een homogeen lineair stelsel Y 0 = A(x)Y
(5.13)
We zullen hierbij onderstellen dat A : (a, b) → Mnn (R) differentieerbaar is. Stelling 5.2.2 De oplossingen van het homogeen stelsel (5.13) vormen een vectorruimte. Bewijs. Y = 0 is steeds een oplossing, zodat de oplossingsverzameling niet leeg is. Om aan te tonen dat de oplossingen van (5.13) een deelvectorruimte van de vectorruimte van de functies (a, b) → Rn vormen, volstaat het aan te tonen dat een lineaire combinatie van twee oplossingen opnieuw een oplossing is. Neem twee oplossingen Y1 en Y2 van (5.13), en α, β ∈ R. Dan is (αY1 + βY2 )0 = αY10 + βY20 = αA(x)Y1 + βA(x)Y2 = A(x)(αY1 + βY2 ) en dus is αY1 + βY2 ook een oplossing van (5.13).
Lemma 5.2.3 Als Y een oplossing is van (5.13), en Y (x0 ) = 0 voor een zekere x0 ∈ (a, b), dan is Y (x) = 0 voor elke x ∈ (a, b).
70
Bewijs. Uit de existentiestelling stelling 5.2.1 volgt dat de vergelijking (5.13) met beginvoorwaarde Y (x0 ) = 0 een unieke oplossing heeft. Y (x) = 0 is een oplossing die aan de beginvoorwaarde voldoet, en dus is het de enige mogelijke oplossing. In stelling 5.2.2 hebben we gezien dat de oplossingsverzameling Opl(Y 0 = A(x)Y ) een vectorruimte is. We kunnen dus spreken van lineaire (on)afhankelijke oplossingen. Een stel oplossingen Y1 ,Y2 , · · · ,Ym is lineair afhankelijk als er constanten c1 , c2 , · · · , cm bestaan waarvan er tenminste e´ e´ n verschillend van nul is zodat m
∑ ciYi(x) = 0
(5.14)
i=1
voor elke x ∈ (a, b). Het stel Y1 ,Y2 , · · · ,Ym is lineair onafhankelijk als (5.14) zich slechts kan voordoen indien c1 = c2 = · · · = cm = 0 Lemma 5.2.4 Zij x0 ∈ (a, b). Een stel oplossingen Y1 ,Y2 , · · · ,Ym van (5.13) is lineair (on)afhankelijk als en alleen als Y1 (x0 ),Y2 (x0 ), · · · ,Ym (x0 ) lineair (on)afhankelijk is in Rn . Bewijs. Onderstel eerst dat het stel oplossingen Y1 ,Y2 , · · · ,Ym lineair afhankelijk is. Er bestaan constanten c1 , c2 , · · · , cm , niet allen nul, zodat (5.14) geldt voor elke x ∈ (a, b). In het bijzonder geldt (5.14) voor x = x0 , en dus is m
∑ ciYi(x0) = 0
i=1
en Y1 (x0 ),Y2 (x0 ), · · · ,Ym (x0 ) is lineair afhankelijk in Rn . Omgekeerd, onderstel dat Y1 (x0 ),Y2 (x0 ), · · · ,Ym (x0 ) lineair afhankelijk is in Rn . Er bestaan constanten c1 , c2 , · · · , cm , niet allen nul, zodat m
∑ ciYi(x0) = 0
i=1
Uit lemma 5.2.3 volgt dat nu dat m
∑ ciYi = 0
i=1
en dus zijn Y1 ,Y2 , · · · ,Ym lineair afhankelijk.
Gevolg 5.2.5 De dimensie van de oplossingsruimte Opl(Y 0 = A(x)Y ) is ten hoogste n. Bewijs. Omdat de dimensie van Rn n is, zijn er in Rn geen stellen van meer dan n lineair onafhankelijke vectoren. Uit lemma 5.2.4 volgt dat dezelfde eigenschap geldt voor Opl(Y 0 = A(x)Y ), en dus is de dimensie van deze ruimte ten hoogste n. Stelling 5.2.6 Zij x0 ∈ (a, b). Er bestaat een omgeving I = (x0 − δ, x0 + δ) van x0 zodanig dat de vectorruimte van de oplossingen van Y 0 = A(x)Y op de omgeving I dimensie n heeft.
71
Bewijs. Neem een basis E1 , E2 , · · · , En van Rn (bijvoorbeeld de standaardbasis). Voor elke i ∈ {1, 2, · · · , n} bestaat er een omgeving van x0 waarop het differentiaalstelsel Y 0 = A(x)Y met beginvoorwaarde Y (x0 ) = Ei een unieke oplossing Yi heeft (stelling 5.2.1). Zij I de kleinste van deze n omgevingen. Uit lemma 5.2.4 volgt dat de oplossingen Y1 ,Y2 , · · · ,Yn lineair onafhankelijk zijn. Op I zijn er dus n lineair onafhankelijke oplossingen van (5.13), en dus is de dimensie van de oplossingsruimte op I tenminste n. In gevolg 5.2.5 hebben we gezien dat deze dimensie tenhoogste n is, en het resultaat volgt. Het volledig stelsel We beschouwen nu het volledige stelsel Y 0 = A(x)Y + B(x)
(5.15)
waarbij we onderstellen dat A : (a, b) → Mnn (R) en B : (a, b) → Rn continu zijn. Onderstel ook dat we het homogene stelsel (5.13) kunnen we oplossen. Om (5.15) te integreren volstaat het om een particuliere oplossing te vinden: Stelling 5.2.7 Als Yh de algemene integraal is van (5.13), en Yp een particuliere integraal van (5.15), dan is Y = Yh +Yp de algemene integraal van (5.15). Bewijs. Doe dit zelf als oefening. Het bewijs is volkomen analoog aan dat van stelling 4.3.1.
Om een particuliere integraal van (5.15) te bepalen passen we de methode van de variatie van de constante toe. Onderstel dat Y1 ,Y2 , · · · ,Yn n lineair onafhankelijke oplossingen van het homogeen stelsel (5.13) zijn. De algemene integraal is dan n
Y (x) = ∑ ciYi (x) i=1
We zoeken nu een particuliere integraal van (5.15) van de vorm n
Yp (x) = ∑ ci (x)Yi (x) i=1
Subsitutie in (5.15) levert n
n
n
i=1
i=1
i=1
∑ ci(x)Yi0(x) + ∑ c0i(x)Yi(x) = ∑ ci(x)A(x)Yi(x) + B(x)
72
Aangezien Yi0 = A(x)Yi , krijgen we n
∑ c0i(x)Yi(x) = B(x)
(5.16)
i=1
(5.16) is een stelsel lineaire vergelijkingen met onbekenden c01 , c02 , · · · , c0n . De matrix van dit stelsel is de matrix met kolommen Y1 ,Y2 , · · · ,Yn . Omdat Y1 ,Y2 , · · · ,Yn lineair onafhankelijk zijn is de matrix van het stelsel regulier, en dit voor elke x. (5.16) heeft dus een unieke oplossing. Oplossen van (5.16) geeft c01 , c02 , · · · , c0n . Integreren geeft c1 , c2 , · · · , cn en de particuliere oplossing Yp .
5.3
Lineaire differentiaalvergelijkingen van orde n
Een differentiaalvergelijking van orde n noemen we lineair als ze van de volgende vorm is. y(n) + a1 (x)y(n−1) + · · · + an−1 (x)y0 + an (x)y = b(x)
(5.17)
Hierbij wordt ondersteld dat de functies a1 , · · · , an en b continu zijn over een interval (a, b). Als b(x) = 0, dan noemen we (5.17) een homogene lineaire differentiaalvergelijking. (5.17) is een normale vergelijking. In § 5.1 hebben we gezien dat een normale vergelijking steeds te herleiden is tot een normaal differentiaalstelsel. Hiertoe stellen we y1 y 0 y2 y Y = .. = .. . . y(n−1)
yn
Y voldoet aan het differentiaalstelsel y01 = y2 0 y2 = y3 .. . y0 = yn n−1 y0n = y(n) = −an (x)y1 − an−1 (x)y2 − · · · − a1 (x)yn + b(x)
(5.18)
of Y 0 = A(x)Y + B(x)
(5.19)
met A(x) =
0
1
0
···
0
0
0 .. .
0 .. .
1 .. .
···
0 .. .
0 .. .
0
0
0
···
1
0
0
0
···
0
−an
−an−1
−an−2
· · · −a2
73
en B(x) = 0 1 −a1
0
0 b(x) 0 .. .
(5.20)
Als de differentiaalvergelijking (5.17) homogeen is, dan is het overeenstemmende differentiaalstelsel (5.19) ook homogeen. We kunnen de resultaten uit § 5.2 toepassen op (5.19) en vertalen naar (5.17). Meer bepaald geldt de existentiestelling voor lineaire differentiaalvergelijkingen van orde n. Stelling 5.3.1 Onderstel dat a1 , · · · , an , b : (a, b) → R continu zijn. (n−1) ∈ R bestaat er een unieke oplossing van (5.17) gedefinieerd op Voor x0 ∈ (a, b) en y0 , y00 , · · · , y0 een omgeving van x0 die voldoet aan de beginvoorwaarden y(x0 ) = y0 y0 (x0 ) = y00 .. . (n−1) (n−1) y (x0 ) = y0 Eigenschappen van de homogene vergelijking y(x) is een oplossing is van de homogene vergelijking y(n) + a1 (x)y(n−1) + · · · + an−1 (x)y0 + an (x)y = 0
(5.21)
als en alleen als Y (x) = (y, y0 , · · · , y(n−1) ) een oplossing is van het homogene stelsel Y 0 = A(x)Y
(5.22)
Uit stelling 5.2.2, gevolg 5.2.5, stelling 5.2.6 volgt nu onmiddellijk Stelling 5.3.2 De oplossingen van (5.21) vormen een vectorruimte van dimensie ten hoogste n. Voor elke x0 ∈ (a, b) bestaat er een omgeving (x0 − δ, x0 + δ) waarop de oplossingen van (5.21) een vectorruimte van dimensie n vormen. Wanneer zijn n oplossingen y1 , · · · , yn van (5.21) lineair onafhankelijk? Een nodige en voldoende voorwaarde is dat de n overeenstemmende oplossingen Y1 , · · · ,Yn van (5.22) lineair onafhankelijk zijn. Uit lemma 5.2.4 volgt dat het voldoende is dat Y1 (x0 ), · · · ,Yn (x0 ) lineair onafhankelijk zijn in Rn , voor een zekere x0 ∈ (a, b). Een nodige en voldoende voorwaarde hiervoor is dat de determinant y1 y2 ··· yn 0 y02 ··· y0n y1 6= 0 W (y1 , · · · , yn ) = . .. . . . (n−1) (n−1) (n−1) y y · · · y n 1 2 Deze determinant wordt de determinant van Wronski of de Wronskiaanse determinant genoemd. Stelling 5.3.3 n oplossingen y1 , y2 , · · · , yn van (5.21) zijn lineair afhankelijk als de Wronskiaan W (y1 , · · · , yn ) identiek gelijk aan 0 is over het interval (a, b). Als de Wronskiaan in 1 punt x0 gelijk aan 0 is, dan is hij overal 0.
74
De volledige lineaire vergelijking Om de algemene integraal van de volledige differentiaalvergelijking (5.17) te bepalen volstaat het om de algemene integraal van de homogene vergelijking (5.21) te bepalen, en een particuliere oplossing van de volledige vergelijking (5.17). Dit volgt onmiddellijk uit stelling 5.2.7. Stelling 5.3.4 Als yh de algemene integraal is van (5.21), en y p een particuliere integraal van (5.17), dan is y = yh + y p de algemene integraal van (5.17). Onderstel dat we n lineair onafhankelijke oplossingen y1 , y2 , · · · , yn van de homogene vergelijking (5.21) kunnen bepalen. Om een particuliere integraal van (5.17) te bepalen kunnen we de methode (n−1) ). De van de variatie van de constante toepassen op het stelsel (5.19). Stel Yi = (yi , y0i , · · · , yi algemene integraal van (5.19) is dan n
Yh = ∑ ciYi i=1
en we zoeken een particuliere integraal van de vorm n
Yp = ∑ ci (x)Yi i=1
Het lineair stelsel (5.16) neemt nu de volgende vorm aan: n c0i (x)yi (x) = 0 ∑ i=1 n ∑ c0 (x)y0 (x) = 0 i i i=1
.. . n (n−1) 0 (x) = b(x) ∑ ci (x)yi
(5.23)
i=1
De determinant van (5.23) is juist de Wronskiaanse determinant, en we weten dat die verschillend van nul is, omdat y1 , y2 , · · · , yn lineair onafhankelijk zijn. We kunnen (5.23) dus oplossen, en we vinden c0i (x), en na integratie ci (x). De gezochte particuliere integraal is dan n
y p (x) = ∑ ci (x)yi (x) i=1
Stelling 5.3.5 Als y1 (x) een oplossing is van de differentiaalvergelijking y(n) + a1 (x)y(n−1) + · · · + an−1 (x)y0 + an (x)y = b1 (x) en y2 (x) een oplossing van y(n) + a1 (x)y(n−1) + · · · + an−1 (x)y0 + an (x)y = b2 (x)
75
dan is y1 (x) + y2 (x) een oplossing van y(n) + a1 (x)y(n−1) + · · · + an−1 (x)y0 + an (x)y = b1 (x) + b2 (x)
Bewijs. Oefening
5.4
Constante co¨effici¨enten
In deze paragraaf bestuderen we de differentiaalvergelijking (5.17) in het geval waarin de co¨effici¨enten ai onafhankelijk zijn van x. y(n) + a1 y(n−1) + · · · + an−1 y0 + an y = b(x)
(5.24)
In principe kan men van (5.24) altijd de algemene integraal bepalen. De homogene vergelijking Beschouw eerst het geval n = 1. De homogene vergelijking y0 + a1 y = 0 kan dan ge¨ıntegreerd worden met behulp van de methode der scheiding van veranderlijken. Men vindt gemakkelijk dat y = ce−a1 x Hierdoor ge¨ınspireerd stellen we voor de algemene homogene vergelijking y(n) + a1 y(n−1) + · · · + an−1 y0 + an y = 0
(5.25)
oplossing voor van de vorm y = eλx waarbij λ een te bepalen parameter is. Afleiden geeft y(i) = λi eλx en na substitutie in (5.25) vinden we λx n n−1 e λ + a1 λ + · · · + an = 0 of P(λ) = λn + a1 λn−1 + · · · + an = 0
(5.26)
y = eλx is dus een oplossing van (5.25) indien λ een wortel is van de n-de graadsvergelijking (5.26). Men noemt de veelterm P de karakteristieke veelterm van de differentiaalvergelijking (5.25). De
76
vergelijking (5.26) noemen we de karakteristieke vergelijking. De karakteristieke vergelijking ontstaat door in de differentiaalvergelijking de i-de afgeleide van y te vervangen door de i-de macht van de veranderlijke λ. In principe kunnen we nu de algemene integraal van de homogene vergelijking bepalen. De karakteristieke vergelijking is een n-de graadsvergelijking en heeft dus n wortels. Elke wortel levert een oplossing, en zo verwachten we n lineair onafhankelijke oplossingen, die een basis vormen voor de vectorruimte der oplossingen. Er zijn echter twee moeilijkheden. 1. Het kan zijn dat de karakteristieke veelterm meervoudige wortels heeft; 2. het kan zijn dat de karakteristieke veelterm complexe wortels heeft. Als f een differentieerbare functie is, dan noteert men D f voor de afgeleide functie. De tweede afgeleide kunnen we noteren als D2 f , de derde als D3 f , enzovoort. We kunnen D beschouwen als de differentiaaloperator die een functie omzet in haar afgeleide. Dan is D2 = D ◦ D, D3 = D ◦ D ◦ D, enz. Nog algemener kunnen we schrijven: (D − λ) f = f 0 − λ f Merk op dat (D − λ)(xn eλx ) = nxn−1 eλx
(5.27)
Immers (D − λ)(xn eλx ) = nxn−1 eλx + λxn eλx − λxn eλx = nxn−1 eλx Verder gebruik makend van onze nieuwe notatie kunnen we de differentiaalvergelijking (5.25) herschrijven als volgt P(D)( f ) = 0 (5.28) waarbij P(λ) = λn + a1 λn−1 + · · · + an de karakteristieke veelterm is. Onderstel nu dat λ1 , λ2 , · · · , λr de - eventueel complexe - wortels zijn van de karakteristieke vergelijking, en dat de multipliciteit van λi het getal mi is. Dan kan de karakteristieke veelterm P ontbonden worden tot P(λ) = (λ − λ1 )m1 (λ − λ2 )m2 · · · (λ − λr )mr (5.28) kan nu herschreven worden onder de vorm P(D)( f ) = (D − λ1 )m1 (D − λ2 )m2 · · · (D − λr )mr f = 0 We beweren nu dat, voor m < mi ,
y = xm eλi x
77
(5.29)
een oplossing is van (5.28). Immers, Uit (5.27) volgt dat (D − λi )m (xm eλi x ) = m!eλi x en dus is (D − λi )mi (xm eλi x ) = 0 en P(D)(xm eλi x ) = 0 We hebben dus in het totaal m1 + m2 + · · · mr = n oplossingen van (5.28), namelijk eλ1 x , xeλ1 x , x2 eλ1 x , · · · , xm1 −1 eλ1 x eλ2 x , xeλ2 x , x2 eλ2 x , · · · , xm2 −1 eλ2 x .. .
(5.30)
eλr x , xeλr x , x2 eλr x , · · · , xmr −1 eλr x We gaan nu aantonen dat deze n oplossingen lineair onafhankelijk zijn. Lemma 5.4.1 Als P een veelterm is van graad n, en λ 6= 0, dan is D(P(x)eλx ) = Q(x)eλx waarbij Q ook een veelterm van graad n is. Bewijs. D(P(x)eλx ) = (λP(x) + P0 (x))eλx en de graad van λP + P0 is dezelfde als die van P.
Lemma 5.4.2 Onderstel dat λ1 , λ2 , · · · , λr niet aan elkaar gelijk zijn en dat P1 , P2 , · · · , Pr veeltermen zijn. Als P1 (x)eλ1 x + P2 (x)eλ2 x + · · · + Pr (x)eλr x = 0 (5.31) voor elke x, dan zijn de veeltermen P1 , P2 , · · · , Pr allemaal nul. Bewijs. Onderstel dat niet alle veeltermen nul zijn, bijvoorbeeld P1 6= 0. We vermenigvuldigen (5.31) met e−λr x . P1 (x)e(λ1 −λr )x + P2 (x)e(λ2 −λr )x + · · · + Pr (x) = 0 (5.32) We leiden deze betrekking gr(Pr ) + 1 maal af (bij conventie onderstellen we dat de graad van de nulveelterm −1 is). We vinden dan een betrekking van de vorm Q1 (x)e(λ1 −λr )x + Q2 (x)e(λ2 −λr )x + · · · + Qr−1 (x)e(λr−1 −λr )x = 0
(5.33)
Uit lemma 5.4.1 volgt dat gr(Qi ) = gr(Pi ), en in het bijzonder is Q1 6= 0. Als we dezelde procedure nogmaals toepassen, dan vinden we R1 (x)e(λ1 −λr−1 )x + R2 (x)e(λ2 −λr−1 )x + · · · + Rr−2 (x)e(λr−2 −λr−1 )x = 0
78
(5.34)
waarbij Ri een veelterm is van dezelde graad als Qi en Pi . Als we de procedure r − 1 keer herhalen vinden we tenslotte S(x)e(λ1 −λ2 )x = 0 (5.35) waarbij S een veelterm is met dezelfde graad als P1 (en dus S 6= 0). Uit (5.35) volgt echter ook dat S = 0, en dit is een contradictie. Gevolg 5.4.3 De n oplossingen (5.30) van (5.28) zijn lineair onafhankelijk. Voorbeelden 5.4.4 1) y00 − y = 0. De karakteristieke vergelijking is λ2 − 1 = 0 en de wortels zijn λ = ±1. y = ex en y = e−x zijn dus oplossingen, en uit gevolg 5.4.3 weten we dat deze lineair onafhankelijk zijn. De algemene integraal is dus y = c1 ex + c2 e−x 2) y00 + 2y0 + y = 0. De karakteristieke vergelijking is nu λ2 + 2λ + 1 = 0 of (λ + 1)2 = 0 −1 is een dubbele wortel van de karakteristieke vergelijking, en y = e−x en y = xe−x zijn twee lineair onafhankelijke oplossingen van de differentiaalvergelijking. De algemene integraal is y = (c1 x + c2 )e−x 3) y(4) − 6y000 + 13y00 − 12y0 + 4y = 0. De karakteristieke vergelijking is nu λ4 − 6λ3 + 13λ2 − 12λ + 4 = 0 of (λ − 1)2 (λ − 2)2 = 0 De algemene integraal is y = c1 ex + c2 xex + c3 e2x + c4 xe2x 4) y(n) = 0. De karakteristieke vergelijking is λn = 0. λ = 0 is de enige wortel, met multipliciteit n. De algemene integraal is y = c1 xn−1 + c2 xn−2 + · · · + cn De hierboven geschetste methode laat toe om van een homogene differentiaalvergelijking met constante co¨effici¨enten de algemene integraal te bepalen. De enige beperking tot nu toe is dat de wortels van de karakteristieke vergelijking re¨eel moeten zijn. Wat gebeurt er als de karakteristieke vergelijking complexe wortels heeft? Formeel is er geen enkel probleem, als we gebruik maken van de complexe exponenti¨ele functie (zie volume 1). Herinner dat ea+ib = ea (cos b + i sin b)
79
Herhaal ook dat de formule
d λx e = λeλx dx blijft gelden als λ een complex getal is. De procedure hierboven geschetst levert dus, ook in het geval dat de karakteristieke vergelijking complexe wortels heeft, n lineair onafhankelijke oplossingen van de differentiaalvergelijking (5.25). Alleen staan deze oplossingen genoteerd in complexe vorm, en voor vele toepassingen hebben we liever de re¨ele vorm. We kunnen echter makkelijk overgaan naar de re¨ele vorm als volgt. Als λ = α + iβ een complexe wortel is van de karakteristieke vergelijking, dan is ook de complex toegevoegde λ = α − iβ een wortel, en wel met dezelfde multipliciteit. Als i kleiner is dan deze multipliciteit, dan zijn xi eλx = xi eαx (cos βx + i sin βx) en xi eλx = xi eαx (cos βx − i sin βx) oplossingen. Lineaire combinaties hiervan (eventueel met complexe co¨effici¨enten) zijn ook oplossingen, en dus zijn xi eαx cos βx en xi eαx sin βx oplossingen van (5.25). Deze oplossingen staan in re¨ele vorm, en op deze manier verkrijgen we n lineair onafhankelijke oplossingen van (5.25). Voorbeelden 5.4.5 1) y00 + y = 0. De karakteristieke vergelijking is λ2 + 1 = 0 De wortels zijn i en −i. De complexe oplossingen zijn cos x + i sin x en cos x − i sin x of, in re¨ele vorm, cos x en sin x. De algemene integraal is dus y = c1 cos x + c2 sin x 2) y00 + 2y0 + 2y = 0. Karakteristieke vergelijking: λ2 + 2λ + 2 = 0 De wortels zijn λ = −1 ± i, en de algemene integraal is y = e−x (c1 cos x + c2 sin x) 3) y000 − y = 0. Karakteristieke vergelijking: λ3 − 1 = 0 De wortels zijn
√ −1 ± i 3 λ = 1 en λ = 2
80
en de algemene integraal is x
− 2x
y = c1 e + e
√ √ 3 3 x + c3 sin x) (c2 cos 2 2
4) yiv + 2y00 + y = 0. Karakteristieke vergelijking: λ4 + 2λ2 + 1 = 0 of (λ2 + 1)2 = 0 We hebben twee dubbele wortels λ = ±i, en de algemene integraal is y = (c1 + c2 x) cos x + (c3 + c4 x) sin x De volledige vergelijking We bekijken nu de volledige vergelijking y(n) + a1 y(n−1) + · · · + an−1 y0 + an y = b(x)
(5.36)
Om een particuliere integraal te bepalen kunnen we altijd de methode van de variatie van de constante toepassen. Voorbeeld 5.4.6 y00 + y = ex . De oplossing van de homogene vergelijking is yh = A cos x + B sin x Als particuliere oplossing stellen we voor: y p = A(x) cos x + B(x) sin x A0 en B0 zijn dan de oplossingen van het lineair stelsel (zie (5.23)) 0 A cos x + B0 sin x = 0 −A0 sin x + B0 cos x = ex Oplossen geeft 0 sin x A0 (x) = x = −ex sin x e cos x cos x 0 0 B (x) = = ex cos x x − sin x e Integratie geeft 1 1 A(x) = (cos x − sin x)ex en B(x) = (cos x + sin x)ex 2 2 Voor de particuliere oplossing vinden we dus yp =
1 1 ex cos x(cos x − sin x)ex + sin x(cos x + sin x)ex = 2 2 2
81
Zoals uit dit eenvoudig voorbeeld blijkt is de methode nogal omslachtig. Als b(x) van een speciale vorm is, dan bestaan er meer eenvoudige technieken. Vooraleer we het meest algemene geval bespreken, geven we eerst enkele eenvoudige voorbeelden. Voorbeelden 5.4.7 1) We hernemen het voorgaande voorbeeld y00 + y = ex . De techniek bestaat er in om een particuliere oplossing voor te stellen die van ”dezelfde vorm¨ıs als het rechterlid. We proberen y p = Aex met A een te bepalen constante. Substitutie in de vergelijking geeft Aex + Aex = ex waaruit volgt dat A = 1/2. De algemene integraal van de differentiaalvergelijking is 1 y = ex + A cos x + B sin x 2 2) y00 + y = x2 + x + 1. Het rechterlid is een veelterm van graad 2 in x, en als particuliere oplossing stellen we een veelterm van graad 2 voor: y p = Ax2 + Bx +C Invullen in de vergelijking levert 2A + Ax2 + Bx +C = x2 + x + 1 waaruit A = 1, B = 1, 2A +C = 1 of C = −1 De particuliere integraal is dus y p = x2 + x − 1 en de algemene integraal van de differentiaalvergelijking is y = x2 + x − 1 + A cos x + B sin x 3) y00 − y = ex . Zoals in voorbeeld 1) proberen we y p = Aex Substitutie in de vergelijking levert Aex − Aex = ex of 0 = 1. In dit geval werkt onze methode dus niet! Eigenlijk is dit niet verwonderlijk, aangezien yh = Cex + De−x Onze voorgestelde particuliere oplossing was dus een oplossing van de homogene vergelijking, en a fortiori geen oplossing van de volledige vergelijking. Om toch een particuliere oplossing
82
te vinden gebruiken we een truucje dat we verderop zullen veralgemenen. Stel als particuliere oplossing voor: y p = Axex Invullen in de vergelijking levert nu A(x + 2)ex − Axex = ex waaruit volgt dat A = 1/2. De algemene integraal van de differentiaalvergelijking is dus x y = (C + )ex + De−x 2 Onderstel dat het rechterlid van de vergelijking (5.36) van de vorm b(x) = V (x)eµx is, waarbij V een veelterm van graad p is. De karakteristieke vergelijking van de homogene vergelijking is λn + a1 λn−1 + · · · + an−1 λ + an = 0 (5.37) Onderstel dat µ, λ1 , · · · , λr de wortels van (5.37) zijn, met multipliciteiten respectievelijk m, m1 , · · · , mr . Als µ geen wortel is van de karakteristieke vergelijking, dan stellen we gewoon m = 0, en al hetgeen volgt blijft geldig. (5.36) kan dan herschreven worden als volgt: (D − µ)m (D − λ1 )m1 · · · (D − λr )mr y = V (x)eµx
(5.38)
Onderstel dat y p een oplossing is van (5.36). We passen op beide leden van (5.38) de operator (D − µ) p+1 . We vinden, gebruik makend van (5.27), dat (D − µ)m+p+1 (D − λ1 )m1 · · · (D − λr )mr y p = (D − µ) p+1 (V (x)eµx ) = 0 y p is dus de oplossing van een homogene differentiaalvergelijking met constante co¨effici¨enten. We kunnen besluiten dat y p kan geschreven worden onder de vorm r y p = c1 xm+p + c2 xm+p−1 + · · · + cm+p x + cm+p+1 eµx + ∑ Qi (x)eλi x i=1
waarbij de ci constanten zijn, en Qi veeltermen van graad kleiner dan mi . Nu is r c p+2 xm−1 + c p+3 xm−2 + · · · + cm+p x + cm+p+1 eµx + ∑ Qi (x)eλi x i=1
een oplossing van de homogene vergelijking. Dus is (c1 xm+p + c2 xm+p−1 + · · · + c p+1 xm )eµx = xm (c1 x p + c2 x p−1 + · · · + c p+1 )eµx een particuliere oplossing van de volledige vergelijking (5.24). We kunnen onze resultaten samenvatten als volgt:
83
Stelling 5.4.8 Beschouw de differentiaalvergelijking y(n) + a1 y(n−1) + · · · + an−1 y0 + an y = V (x)eµx
(5.39)
waarbij V (x) een veelterm van graad p is. Onderstel dat de multipliciteit van µ als wortel van de karakteristieke vergelijking λn + a1 λn−1 + · · · + an−1 λ + an = 0 m is. Dan bestaat er een particuliere oplossing van (5.39) van de vorm y p = xmV1 (x)eµx waarbij V1 een veelterm van dezelfde graad als V is. Voorbeelden 5.4.9 1) y000 + y00 = x2 . De oplossing van de homogene vergelijking is yh = c1 e−x + c2 + c3 x 0 is dus een dubbele wortel van de karakteristieke vergelijking. Het rechterlid bestaat uit e0x vermenigvuldigd met een veelterm van graad 2. Daarom stellen we een particuliere oplossing voor van de vorm y p = x2 (Ax2 + Bx +C) Invullen in de vergelijking geeft y0p = 4Ax3 + 3Bx2 + 2Cx y00p = 12Ax2 + 6Bx + 2C y000 p = 24Ax + 6B Invullen in de vergelijking levert 12Ax2 + (6B + 24A)x + 6B + 2C = x2 waaruit A=
1 1 ; B=− ; C=1 12 3
en de algemene integraal is y = c1 e−x + c2 + c3 x + x2 −
x3 x4 + 3 12
2) y00 −2y0 +y = ex . λ = 1 is een dubbele wortel van de karakteristieke vergelijking, en de oplossing van de homogene vergelijking is yh = (c1 x + c2 )ex Als particuliere oplossing stellen we voor y p = Ax2 ex
84
Afleiden geeft y0p = A(x2 + 2x)ex y00p = A(x2 + 4x + 2)ex Invullen in de vergelijking levert 2Aex = ex , en A = 1/2. De algemene integraal is dus y=(
x2 + c1 x + c2 )ex 2
3) y00 − y = cos x. De algemene integraal van de homogene vergelijking is yh = c1 ex + c2 e−x Het rechterlid is cos x = 21 (eix + e−ix ). Als particuliere oplossing stellen we daarom een lineaire combinatie van eix en e−ix voor, of, hetgeen op hetzelfde neerkomt, een lineaire combinatie van cos x en sin x: y p = A cos x + B sin x Afleiden geeft y00p = −A cos x − B sin x Invullen in de vergelijking levert −2A = 1 en 2B = 0. 1 y p = − cos x 2 is dus een particuliere oplossing van de vergelijking. 4) y00 + y = sin x. De oplossing van de homogene vergelijking is nu yh = c1 cos x + c2 sin x i en −i zijn enkelvoudige wortels van de karakteristieke vergelijking. Als particuliere oplossing stellen we nu een lineaire combinatie van xeix en xe−ix voor, of, hetgeen op hetzelfde neerkomt, een lineaire combinatie van x cos x en x sin x: y p = Ax cos x + Bx sin x Afleiden geeft y0p = A cos x + B sin x − Ax sin x + Bx cos x y00p = −2A sin x + 2B cos x − Ax cos x − Bx sin x Invullen in de vergelijking levert −2A sin x + 2B cos x = sin x zodat A = −1/2 en B = 0. De algemene integraal van de vergelijking is x yh = (c1 − ) cos x + c2 sin x 2
85
5.5
De methode der eigenwaarden en eigenvectoren
We beschouwen een homogeen differentiaalstelsel met constante co¨effici¨enten. y01 = a11 y1 + a12 y2 + · · · + a1n yn y0 = a21 y1 + a22 y2 + · · · + a2n yn 2 .. . y0 = a y + a y + · · · + a y n
n1 1
n2 2
(5.40)
nn n
In matrixvorm: Y 0 = AY
(5.41)
Om dit stelsel op te lossen kunnen we altijd de methode van afleiding en eliminatie toepassen. Dit herleidt het stelsel tot een homogene differentiaalvergelijking met constante co¨effici¨enten. Het volstaat dan om de resultaten uit § 5.4 toepassen. In deze paragraaf be bespreken we in het kort een alternatieve methode. Herhaal dat V een eigenvector is van de matrix A met bijhorende eigenwaarde λ als AV = λV In de cursus lineaire algebra hebben we gezien hoe we de eigenwaarden en eigenvectoren van een matrix A kunnen bepalen. Lemma 5.5.1 Als V een eigenvector is van A met bijhorende eigenwaarde λ, dan is Y = Veλx een oplossing van (5.41) Bewijs. Y 0 = λVeλx = AVeλx = AY Onderstel nu dat er n lineair onafhankelijke eigenvectoren V1 ,V2 , · · · ,Vn van A bestaan, met bijhorende eigenwaarden λ1 , λ2 , · · · , λn . Dit is steeds het geval als A een symmetrische matrix is (zie cursus lineaire algebra). In dit geval is de algemene integraal van (5.41) n
Y = ∑ ciVi eλi x i=1
De oplossing van (5.41) is dus herleid tot het bepalen van de eigenwaarden en eigenvectoren van de matrix van het stelsel.
86
Hoofdstuk 6 Oplossing van differentiaalvergelijkingen door reeksontwikkeling 6.1
Inleiding
We beschouwen opnieuw de homogene lineaire differentiaalvergelijking a0 (x)y(n) + a1 (x)y(n−1) + · · · + an−1 (x)y0 + an (x)y = 0
(6.1)
In het vorige hoofdstuk hebben we gezien dat we deze altijd kunnen oplossen als de co¨effici¨enten ai constanten zijn. In het algemeen weten we dat de oplossingen een vectorruimte vormen. In de omgeving van een gegeven punt a is de dimensie van de oplossingsruimte juist n, de orde van de differentiaalvergelijking. In het algemeen kan men de oplossingen van (6.1) niet exact bepalen. Daarom zoekt men naar benaderde oplossingen. In dit hoofdstuk is het idee om oplossingen te zoeken in de omgeving van een punt a, in de vorm van een machtreeks ∞
y(x) =
∑ cn(x − a)n
(6.2)
n=0
We zullen zien hoe men de co¨effici¨enten cn recursief kan bepalen. Om de berekeningen te vereenvoudigen zullen we ons beperken tot homogene lineaire vergelijkingen van orde 2. Onze methoden zijn ook voor hogere orden van toepassing. De vergelijking die we bestuderen is dus van de vorm y00 + a1 (x)y0 + a2 (x)y = 0
(6.3)
De methode heeft natuurlijk alleen maar zin als er oplossingen van de vergelijking bestaan waarvan de Taylorreeks in het punt a bestaat en convergeert in een omgeving van een punt a. Daarom voeren we eerst het volgende begrip in. Definitie 6.1.1 Een numerieke functie f : R → R noemen we analytisch in het punt a als de Taylorreeks in het punt a bestaat en convergeert naar f (x) voor x in een omgeving van a.
87
Merk op dat dit een zware eis is voor de functie f . Een functie die analytisch is in het punt a is een oneindig aantal keer differentieerbaar in een omgeving van a. Definitie 6.1.2 We noemen a ∈ R een gewoon punt van de differentiaalvergelijking (6.3) als de functies a1 en a2 analytisch zijn in het punt a. Een punt dat geen gewoon punt is noemt men een singulier punt. Definitie 6.1.3 We noemen a ∈ R een regelmatig singulier punt van de differentiaalvergelijking (6.3) indien de vergelijking kan herschreven worden onder de vorm (x − a)2 y00 + (x − a)p(x)y0 + q(x)y = 0
(6.4)
waarbij de functies p en q analytisch zijn in het punt a. In de volgende paragrafen zullen we zien dat de hierboven geschetste methode kan gebruikt worden in de omgeving van een gewoon punt of een regelmatig singulier punt.
6.2
Oplossing in een omgeving van een gewoon punt
Zonder bewijs geven we de volgende algemene stelling. Stelling 6.2.1 Onderstel dat a een gewoon punt is van (6.3), en dat r het minimum is van de convergentiestralen van de Taylorreeksen van a1 (x) en a2 (x). Dan is elke oplossing van (6.3) analytisch in a, en de Taylorreeks ervan convergeert zeker voor |x − a| < r. Hoe bepalen we deze oplossing? We schrijven ∞
∑ cn(x − a)n
y(x) =
n=0
en vullen dit in in (6.3). Binnen het convergentieinterval geldt y0 (x) = y00 (x) =
∞
∞
n=0 ∞
n=0
∑ ncn(x − a)n−1 = ∑ (n + 1)cn+1(x − a)n n−2
∑ n(n − 1)cn(n − a)
∞
=
n=0
∑ (n + 2)(n + 1)cn+2(x − a)n
n=0
Voor |x − a| < r geldt tevens ∞
a1 (x) =
∑ αk (x − a)k k=0 ∞
a2 (x) =
∑ βk (x − a)k
k=0
88
(6.5)
Als we dit allemaal in (6.3) invullen vinden we ∞
∑ (n + 2)(n + 1)cn+2(x − a)n
n=0
"
#"
∞
#
∞
+ ∑ αk (x − a)k ∑ (n + 1)cn+1 (x − a)n #" n=0 # "k=0 ∞
+
∑ βk (x − a)k k=0
∞
∑ cn(x − a)n
= 0
n=0
Het linkerlid is een som van producten van machtreeksen. Zoals we gezien hebben, mogen we machtreeksen optellen en vermenigvuldigen zoals veeltermen, zolang x in het inwendige van het convergentiegebied van alle machtreeksen blijft. Het linkerlid is dus nul als de co¨effici¨ent van xn nul is voor elke n, dus als n
n
(n + 2)(n + 1)cn+2 + ∑ αk (n − k + 1)cn−k+1 + ∑ βk cn−k = 0 k=0
(6.6)
k=0
(6.6) bepaalt cn+2 recursief in functie van c0 , c1 , c2 , · · · , cn+1 , en dit voor n = 0, 1, 2, · · ·. De eerste vergelijking bepaalt dus c2 in functie van c0 en c1 , de tweede c3 in functie van c0 , c1 en c2 , enz. c0 en c1 kunnen we beschouwen als integratieconstanten (de vergelijking is van orde 2, en er zijn dus twee integratiecostanten). Indien men de oplossing wenst die voldoet aan de beginvoorwaarden y(a) = y0 en y0 (a) = y1 dan volstaat het c0 = y0 en c1 = y1 te stellen. De vergelijking van Legendre Dit is een vergelijking van de vorm (1 − x2 )y00 − 2xy0 + α(α + 1)y = 0
(6.7)
Het punt x = 0 is een gewoon punt van deze vergelijking. Dit ziet men door (6.7) te delen door 1 − x2 . Men krijgt x α(α + 1) y00 − 2 y0 + y=0 2 1−x 1 − x2 De Taylorreeksen a1 (x) = −2 a2 (x) =
∞ x = −2 ∑ x2n+1 1 − x2 n=0
∞ α(α + 1) = α(α + 1) ∑ x2n 1 − x2 n=0
89
convergeren voor |x| < 1. Uit stelling 6.2.1 volgt dat de reeksontwikkeling van de oplossing ook convergeert voor |x| < 1. Stel opnieuw ∞
y=
∑ cnxn
n=0
We subsitueren deze reeks in (6.7). Dit geeft ∞ ∞ ∞ (1 − x2 ) ∑ n(n − 1)cn xn−2 − 2x ∑ ncn xn−1 + α(α + 1) ∑ cn xn = 0 of
n=0
n=1
n=2
∞
∞
∞
∞
n=2
n=1
n=0
∑ (n + 2)(n + 1)cn+2xn − ∑ n(n − 1)cnxn − 2 ∑ ncnxn + α(α + 1) ∑ cnxn = 0
n=0
De co¨effici¨ent van
xn
moet 0 zijn. Dit geeft
(n + 2)(n + 1)cn+2 − n(n − 1)cn − 2ncn + α(α + 1)cn = 0 of cn+2 =
(n − α)(n + α + 1) cn (n + 2)(n + 1)
(6.8)
De recursiebetrekking (6.8) bepaalt de co¨effici¨enten c3 , c4 , · · · in functie van c0 en c1 . Als we c0 = 1 en c1 = 0 stellen verkrijgen we een eerste particuliere oplossing y1 (x) = 1 −
α(α + 1) 2 α(α − 2)(α + 1)(α + 3) 4 x + x +··· 2! 4!
Deze oplossing convergeert voor |x| < 1 en voldoet aan de beginvoorwaarden y1 (0) = c0 = 1 en y01 (0) = c1 = 0 Als we c0 = 0 en c1 = 1 stellen verkrijgen we een tweede particuliere oplossing y2 (x) = x −
(α − 1)(α + 2) 3 (α − 1)(α − 3)(α + 2)(α + 4) 5 x + x +··· 3! 5!
Deze oplossing convergeert eveneens voor |x| < 1 en voldoet aan de beginvoorwaarden y1 (0) = c0 = 0 en y01 (0) = c1 = 1 De algemene integraal van de vergelijking van Legendre is nu y(x) = c0 y1 (x) + c1 y2 (x)
90
De veeltermen van Legendre Als α = n een natuurlijk getal is, dan is een van de twee oplossingen van de vergelijking van Legendre een veelterm. Als α = n een even getal is, dan volgt uit de recursiebetrekking dat cn+2 = 0. Dan volgt ook dat cn+4 = cn+6 = · · · = 0, zodat de oplossing y1 een veelterm van graad n is. Als α = n oneven is, dan volgt uit de recursiebetrekking nog steeds dat cn+2 = cn+4 = cn+6 = · · · = 0, zodat de oplossing y2 een veelterm van graad n is. Voor elke n ∈ N vinden we dus een veelterm van graad n. Deze is bepaald op een veelvoud na. Als we de hoogstegraadsco¨effici¨ent van deze veelterm gelijk stellen aan 1 (2n)! 2n (n!)2 dan verkrijgen we de veeltermen van Legendre: Definitie 6.2.2 Voor n ∈ N noemen we ! 1 (2n)! n n(n − 1) n−2 n(n − 1)(n − 2)(n − 3) n−4 Pn (x) = n x + x −··· x − 2 (n!)2 2(2n − 1) 2.4(2n − 1)(2n − 3)
(6.9)
de n-de veelterm van Legendre. Voor lage waarden van n kunnen we de veeltermen van Legendre gemakkelijk berekenen. P0 (x) = 1 P1 (x) = x 3x2 − 1 P2 (x) = 2 5x3 − 3x P3 (x) = 2 35x4 − 30x2 + 3 P4 (x) = 8 De veeltermen van Legendre hebben vele toepassingen. We formuleren hier alvast enkele eigenschappen, we zullen in een later hoofdstuk op de veeltermen van Legendre terugkomen. Stelling 6.2.3 Als n 6= m, dan is Z 1 −1
Pn (x)Pm (x)dx = 0
Bewijs. Pn is een oplossing van de differentiaalvergelijking van Legendre (1 − x2 )Pn00 (x) − 2xPn0 (x) + n(n + 1)Pn (x) = 0 of
i dh 2 0 (1 − x )Pn (x) + n(n + 1)Pn (x) = 0 dx
91
(6.10)
Op dezelfde manier geldt i dh (1 − x2 )Pm0 (x) + m(m + 1)Pm (x) = 0 dx
(6.11)
We vermenigvuldigen (6.10) met Pm (x) en (6.11) met Pn (x) en trekken af. i i dh dh 2 0 2 0 (1 − x )Pm (x) − Pm (x) (1 − x )Pn (x) n(n + 1) − m(m + 1) Pn (x)Pm (x) = Pn (x) dx dx of
i dh (1 − x2 )(Pn (x)Pm0 (x) − Pm (x)Pn0 (x)) dx We integreren beide leden van −1 tot 1. Dan volgt (n − m)(n + m + 1)Pn (x)Pm (x) =
(n − m)(n + m + 1)
Z 1 −1
=
Z 1 dh −1
=
dx
(1 − x
2
Pn (x)Pm (x)dx i dx
)(Pn (x)Pm0 (x) − Pm (x)Pn0 (x))
h i1 (1 − x2 )(Pn (x)Pm0 (x) − Pm (x)Pn0 (x)) =0 −1
Voor n 6= m volgt hieruit het gestelde.
Gevolg 6.2.4 Als p(x) een veelterm is van graad kleiner dan n, dan is Z 1 −1
p(x)Pn (x)dx = 0
Bewijs. De veeltermen {P0 , P1 , · · · , Pn−1 } vormen een basis van de vectorruimte der veeltermen van graad kleiner dan n. Immers, de overgangsmatrix naar de standaardbasis {1, x, x2 , · · · , xn−1 } is een driehoeksmatrix, en is dus inverteerbaar. Een veelterm p(x) van graad kleiner dan n kan dus geschreven worden als een lineaire combinatie n−1
p(x) =
∑ aiPi(x)
i=0
en het gestelde volgt onmiddellijk uit stelling 6.2.3. Stelling 6.2.5 Pn (x) =
1 dn 2 (x − 1)n 2n n! dxn
Bewijs. We rekenen het rechterlid uit. dn 2 (x − 1)n dxn
92
! dn n n = x2n − x2n−2 + x2n−4 − · · · dxn 1 2 n (2n − 2)! n−2 n (2n − 4)! n−4 (2n)! n x − x + x −··· = 1 (n − 2)! 2 (n − 4)! n! ! (2n)! n n(n − 1) n−2 n(n − 1)(n − 2)(n − 3) n−4 = x − x + x −··· n! 2(2n − 1) 2.4(2n − 1)(2n − 3) Als we dit vergelijken met definitie 6.2.2 dan volgt het resultaat.
Op analoge manier kan men de veeltermen van Laguerre en Hermite invoeren. Voor verdere details verwijzen we naar de oefeningen.
6.3
Oplossing in een omgeving van een regelmatig singulier punt
Het idee is om de voorgestelde oplossing ∞
y=
∑ cn(x − a)n
n=0
in het geval van een gewoon punt te vervangen door ∞
y=
∑ cn(x − a)n+ρ
(6.12)
n=0
waarbij ρ een getal is dat re¨eel en zelfs complex kan zijn. Voor x < a kunnen er problemen zijn (bijvoorbeeld als ρ = 1/2), en daarom gebruiken we (6.12) enkel voor x > a. Voor x < a stellen we een oplossing voor van de vorm ∞
y=
∑ cn(a − x)n+ρ
(6.13)
n=0
In het vervolg zullen we ons beperken tot het geval x > a. We zoeken dus oplossing op de rechterhelft van een omgeving van a. Dat de hierboven geschetste methode werkt, blijkt uit de volgende stelling die we hier zonder bewijs geven. Stelling 6.3.1 Onderstel dat a een regelmatig singulier punt is van de vergelijking (x − a)2 y00 + (x − a)p(x)y0 + q(x)y = 0
(6.14)
en dat r het minimum is van de convergentiestralen van de Taylorreeksen van p(x) en q(x). Dan bestaat er minstens een oplossing van (6.14) van de vorm (6.12), en deze reeks convergeert zeker voor a < x < a + r.
93
Om de co¨effici¨enten cn praktisch te berekenen gaan we tewerk zoals in de vorige paragraaf. Stel ∞
∑ pn(x − a)n
p(x) =
∞
en q(x) =
n=0
∑ qn(x − a)n
n=0
Omdat we machtreeksen term per termen mogen afleiden vinden we ook ∞
y0 =
∑ (n + ρ)cn(x − a)n+ρ−1
n=0 ∞
y00 =
∑ (n + ρ)(n + ρ − 1)cn(x − a)n+ρ−2
n=0
Als we dit allemaal vervangen in (6.14), dan vinden we ∞
∑ (n + ρ)(n + ρ − 1)cn(x − a)n+ρ
n=0
"
#"
∞
#
∞
+ ∑ pn (x − a)n ∑ (n + ρ)cn (x − a)n+ρ "n=0 #" n=0 # ∞
∞
∑ qn(x − a)n
+
∑ cn(x − a)n+ρ
n=0
= 0
n=0
We delen alles door (x − a)ρ en vinden ∞
∑ (n + ρ)(n + ρ − 1)cn(x − a)n
n=0
"
#"
∞
#
∞
+ ∑ pn (x − a)n ∑ (n + ρ)cn (x − a)n "n=0 #" n=0 # ∞
+
∑ qn(x − a)n
n=0
(6.15)
∞
∑ cn(x − a)n
= 0
n=0
Dit is weer een machtreeks in x − a die identiek gelijk is aan 0. Bijgevolg moeten alle co¨effici¨enten gelijk zijn aan 0. Als we kijken naar de constante term van de reeks, dan vinden we c0 ρ(ρ − 1) + p0 ρc0 + q0 c0 = 0 We kunnen onderstellen dat de co¨effici¨ent c0 niet nul is. Immers, als c0 = c1 = · · · = ci−1 = 0, dan volstaat het om de co¨effici¨enten te hernummeren (stel c0m = cm+i ) en ρ te vervangen door ρ + i. We vinden dus ρ2 + (p0 − 1)ρ + q0 = 0 (6.16) Deze vergelijking noemt men de karakteristieke vergelijking. Oplossing van (6.16) geeft de mogelijke waarden van ρ. Als we de co¨effici¨ent van (x − a)n in (6.15) gelijk aan 0 stellen, dan vinden we n
n
(n + ρ)(n + ρ − 1)cn + ∑ pn−s (s + ρ)cs + ∑ qn−s cs = 0 s=0
s=0
94
of
n−1 (n + ρ)(n + ρ − 1) + p0 (n + ρ) + q0 cn + ∑ pn−s (s + ρ) + qn−s cs = 0
(6.17)
s=0
Neem een oplossing ρ van de karakteristieke vergelijking (6.16). (6.17) laat toe om de co¨effici¨enten cn recursief te bepalen in functie van c0 . Dit levert een oplossing van de differentiaalvergelijking. Omdat de karakteristieke vergelijking een kwadratische vergelijking is, zijn er in het algemeen twee oplossingen, ρ1 en ρ2 . Deze geven aanleiding tot twee oplossingen y1 en y2 van de vergelijking. Als deze oplossingen lineair onafhankelijk zijn, dan kennen we dus de algemene integraal van de differentiaalvergelijking. Verschillende gevallen zijn mogelijk: 1) ρ1 en ρ2 zijn verschillend en re¨eel. Als ρ1 − ρ2 ∈ / Z, dan kan men aantonen dat y1 en y2 lineair onafhankelijk zijn. Indien ρ1 − ρ2 ∈ Z, dan kan het zijn dat y1 een veelvoud is van y2 . In dit geval levert onze methode dus slechts een lineair onafhankelijke oplossing. 2) ρ1 en ρ2 zijn toegevoegd complex. Dan zijn y1 en y2 lineair onafhankelijk, maar wel in complexe vorm geschreven. Door geschikte lineaire combinaties van y1 en y2 te nemen vindt men twee lineair onafhankelijke re¨ele oplossingen ; de details worden uitgewerkt in het oefeningenboek. In de praktijk blijkt dit geval niet zoveel voor te komen, zodat we er hier niet verder op ingaan. 3) ρ1 = ρ2 . Onze methode levert dan uiteraard maar e´ e´ n lineair onafhankelijke oplossing. In dit geval bestaan er methodes om een tweede lineair onafhankelijke oplossing in reeksvorm op te schrijven. We verwijzen hier weer naar het oefeningenboek. Voorbeeld: de vergelijking van Bessel We beschouwen de vergelijking x2 y00 + xy0 + (x2 − ν2 )y = 0
(6.18)
Hierin is ν een niet-negatieve parameter. x = 0 is een regelmatig singulier punt van deze vergelijking. We zoeken dus een oplossing van de vorm ∞
y=
∑ cnxn+ρ
n=0
Term per term afleiden geeft y0 =
∞
∑ (n + ρ)cnxn+ρ−1
n=0
en y00 =
∞
∑ (n + ρ)(n + ρ − 1)cnxn+ρ−2
n=0
Als we dit substitueren in (6.18), dan vinden we ∞
∞
∞
n=0
n=0
n=0
∑ (n + ρ)(n + ρ − 1)cnxn+ρ + ∑ (n + ρ)cnxn+ρ + ∑ cnxn+ρ+2 − ν2
95
∞
∑ cnxn+ρ = 0
n=0
of
∞
∞
∞
∞
n=0
n=0
n=2
n=0
∑ (n + ρ)(n + ρ − 1)cnxn + ∑ (n + ρ)cnxn + ∑ cn−2xn − ν2 ∑ cnxn = 0
(6.19)
Als we kijken naar de constante term van deze machtreeks, dan vinden we (ρ(ρ − 1) + ρ − ν2 )c0 = 0 of ρ2 = ν2
(6.20)
De oplossingen van de karakteristieke vergelijking zijn dus ρ = ±ν. We bekijken nu de lineaire term in (6.19). Dit geeft ((1 + ρ)(1 + ρ − 1) + (1 + ρ) − ν2 )c1 = 0 of ((1 + ρ)2 − ν2 )c1 Behalve in het geval waarin ρ = −ν = −1/2 volgt hieruit dat c1 = 0. Als we de co¨effici¨ent van xn in (6.19) gelijk aan nul stellen, dan vinden we ((n + ρ)(n + ρ − 1) + (n + ρ) − ν2 )cn + cn−2 = 0 of (n + ρ)2 − ν2 cn + cn−2 = 0
(6.21)
Dit laat toe om c2 , c4 , c6 , · · · recursief te berkenen uit c0 . Bovendien volgt hieruit dat c1 = c3 = c5 = · · · = 0. Eerste geval Beschouw nu eerst het bijzonder geval waarin ρ1 − ρ2 = 2ν ∈ / N. Uit de algemene theorie volgt dan dat we twee lineair onafhankelijke oplossingen moeten vinden. Beschouw eerst het geval waarin ρ = ν. De recursiebetrekking (6.21) neemt nu de volgende vorm aan: cn = −
cn−2 n(n + 2ν)
Herhaaldelijk toepassen van deze betrekking levert c0 c0 =− 2 2(2 + 2ν) 2 (ν + 1) c2 c2 c0 = − =− 2 = 4 4(4 + 2ν) 2 .2(2 + ν) 2 .2(ν + 1)(ν + 2) .. . c0 = (−1)n 2n 2 n!(ν + 1)(ν + 2) · · · (ν + n)
c2 = − c4
c2n
Hierbij is c0 een willekeurige constante. Als we c0 gepast kiezen, dan kunnen we voor c2n een elegante formule opschrijven. Herhaal dat de gammafunctie gedefinieerd wordt door de formule Z ∞
Γ(p) =
t p−1 e−t dt
0
96
(p > 0)
(6.22)
Een belangrijke eigenschap van de gammafunctie is dat Γ(p) = (p − 1)Γ(p − 1)
(6.23)
Voor details verwijzen we naar volume 2. Als we c0 = stellen, dan vinden we dat c2n =
1 2ν Γ(ν + 1)
(−1)n n!Γ(ν + n + 1)22n+ν
De gevonden oplossing is dus (−1)n x Jν (x) = ∑ n=0 n!Γ(ν + n + 1) 2 ∞
!2n+ν (6.24)
We noemen Jν (x) de Besselfunctie van de eerste soort en orde ν. de Besselfuncties J0, J1, J2 en J3
de Besselfuncties J4, J5, J6 en J7
1
0.8
0.8
0.6
0.6
0.4
0.4
y-as
y-as
1
0.2
0.2
0
0
-0.2
-0.2
-0.4
-0.4 0
2
4
6
8
10
12
14
16
18
0
20
2
x-as
4
6
8
10
12
14
16
18
20
x-as
Figuur 6.1: De Besselfuncties J0 , J1 , J2 , J3 , J4 , J5 , J6 , J7 De oplossing behorend bij ρ = −ν wordt op analoge wijze berekend: het volstaat om overal ν door −ν te vervangen. Het enige probleem is dat Γ(−ν+n+1) niet gedefinieerd is voor −ν+n+1 < 0. De integraal (6.22) is immers divergent voor p ≤ 0. Men lost dit probleem als volgt op: voor p ∈ (−1, 0) definieert men Γ(p) door de formule (6.23): Γ(p) =
Γ(p + 1) p
Vervolgens gebruiken we dezelfde formule om Γ(p) te defini¨eren voor p ∈ (−2, −1). Als we zo doorgaan, dan kunnen we de gammafunctie defini¨eren voor p ∈ R \ {0, −1, −2, −3, · · ·}. We vinden dan een tweede oplossing voor de vergelijking van Bessel: !2n−ν ∞ (−1)n x J−ν (x) = ∑ (6.25) n=0 n!Γ(−ν + n + 1) 2
97
De algemene integraal van de vergelijking van Bessel is dus y = AJν (x) + BJ−ν (x) y=gamma(x)
y=1/gamma(x)
4
6 3
4
2
y-as
y-as
2 0 -2
1
0
-4
-1
-6 -5
-4
-3
-2
-1
0
1
2
3
4
-2
5
-4
-2
x-as
0
2
4
6
x-as
Figuur 6.2: De Gamma functie Tweede geval Indien 2ν ∈ N een geheel getal is, dan kunnen we de bovenstaande redenering herhalen, en we vinden nog steeds de oplossing Jν . Als 2ν een even getal is, dan zijn er problemen bij het bepalen van de tweede oplossing: uit de recursiebetrekking −n(n − 2ν)cn = cn−2 volgt dat c2ν−2 = 0. Maar dan is ook c2ν−4 = c2ν−6 = · · · = c0 = 0. Dit is in strijd met onze hypothese dat c0 6= 0. Als we c0 = 0 toelaten, dan vinden we dat de eerste co¨effici¨ent die verschilt van 0, c2ν is. De gevonden oplossing is dan van de vorm y = c2ν x2ν−ν + · · · = c2ν xν + · · · en dit is juist c2ν Jν (x). Onze methode levert dus slechts e´ e´ n lineair onafhankelijke oplossing. Als 2ν een oneven getal is, dan levert de methode een tweede lineair onafhankelijke oplossing J−ν . Speciale aandacht verdient het geval ν = 1/2. In dit geval is ρ = −1/2 een oplossing van de karakteristieke vergelijking. Zoals we hierboven gezien hebben kunnen we in dit geval niet automatisch besluiten dat c1 = 0. Het is echter eenvoudig om in te zien dat de oplossing die we krijgen door c0 = 0 en c1 = 1 te stellen juist de Besselfunctie J1/2 (x) is.
98
Hoofdstuk 7 Eerste integralen 7.1
Eerste integralen
Beschouw een normaal stelsel differentiaalvergelijkingen van n vergelijkingen met n onbekende functies: y01 = f1 (x, y1 , y2 , · · · , yn ) y0 = f2 (x, y1 , y2 , · · · , yn ) 2 (7.1) .. . y0 = f (x, y , y , · · · , y ) n
n
1
2
n
Een eerste integraal van het stelsel (7.1) is een functie φ(x, y1 , y2 , · · · , yn ) die constant is (dit wil zeggen: onafhankelijk van x) als men y1 , y2 , · · · , yn vervangt door een oplossing van het stelsel. Een eerste integraal kunnen we dus schrijven als φ(x, y1 , y2 , · · · , yn ) = c of, korter, φ=c Het is evident dat als φ = c een eerste integraal is van (7.1), en f : R → R een willekeurige functie is, dat dan f ◦ φ = c0 ook een eerste integraal is. Immers, op een oplossing geldt dat f (φ(x, y1 , y2 , · · · , yn )) = f (c) = c0 Deze nieuwe eerste integraal levert echter geen nieuwe informatie op over de oplossingen van het differentiaalstelsel. Daarom voeren we de volgende definitie in: twee eerste integralen φ1 = c1 en φ2 = c2 worden verschillend genoemd als er geen enkele functie g bestaat zodanig dat φ1 = g ◦ φ2 of φ2 = g ◦ φ1 . Meer algemeen zeggen we dat m eerste integralen φ1 = c1 , · · · , φm = cm verschillend zijn als er er geen enkele functie in m − 1 veranderlijken bestaat zodat φi = g(φ1 , · · · , φi−1 , φi+1 , · · · , φm ). Stelling 7.1.1 Als men een eerste integraal φ(x, y1 , · · · , yn ) = c van het stelsel (7.1) kent, dan kan men het stelsel vervangen door een equivalent stelsel van orde n − 1, waarin ook de constante c voorkomt.
99
Bewijs. Los een van de onbekende functies, bijvoorbeeld y1 op uit de eerste integraal. Men verkrijgt y1 = ψ(x, y2 , · · · , yn ) Substitutie hiervan in n − 1 vergelijkingen van het stelsel (7.1) geeft een stelsel in de n − 1 onbekenden y2 , · · · , yn , dat bovendien afhangt van c. Stelling 7.1.2 Als de algemene integraal van het stelsel (7.1) afhangt van n constanten, dan bezit het stelsel juist n verschillende eerste integralen. Omgekeerd, indien men n verschillende eerste integralen van het stelsel kent, dan kan men hieruit ook de algemene integraal van het stelsel bepalen. Bewijs. Bij onderstelling is de algemene integraal van het stelsel van de vorm y1 = y1 (x, c1 , · · · , cn ) y2 = y2 (x, c1 , · · · , cn ) .. . y = y (x, c , · · · , c ) n
1
n
(7.2)
n
Als we dit stelsel oplossen naar de constanten c1 , · · · , cn dan verkrijgen we (met behulp van de stelling van de impliciete functies): c1 = φ1 (x, y1 , · · · , yn ) c2 = φ2 (x, y1 , · · · , yn ) (7.3) .. . c = φ (x, y , · · · , y ) n
1
n
n
We hebben dus n eerste integralen. Deze n eerste integralen zijn verschillend, want anders bestond er een betrekking van de vorm φ j = ψ(φ1 , · · · , φ j−1 , φ j+1 , · · · , φn ) met andere woorden c j = ψ(c1 , · · · , c j−1 , c j+1 , · · · , cn ) maar dit zou betekenen dat de constanten c1 , · · · , cn niet willekeurig zijn. We moeten nog aantonen dat er geen (n + 1)-de eerste integraal is die verschillend is van de vorige. Onderstel dat cn+1 = φn+1 (x, y1 , · · · , yn ) (7.4) zulk een eerste integraal. (7.3) en (7.4) laten dan toe om x, y1 , · · · , yn te berekenen in functie van de constanten c1 , · · · , cn+1 . Maar dan is x geen onafhankelijke veranderlijke meer. Veronderstel tenslotte dat we n verschillende eerste integralen (7.3) kennen. Hieruit kunnen we y1 , y2 , · · · , yn oplossen, zodat we een stelsel van de vorm (7.2) krijgen. Dit is een algemene integraal. Men kan stelling 7.1.2 als volgt interpreteren: de algemene integraal van (7.1) is in feite een verzameling krommen in Rn die afhangt van n parameters c1 , · · · , cn . Als we deze algemene integraal
100
schrijven onder de vorm (7.2), dan betekent dit dat we van elke oplossing een stel parametervergelijkingen kennen. Als we n verschillende eerste integralen (7.3) kennen, dan betekent dit dat we elke oplossing kennen als de doorsnede van n hyperoppervlakken in Rn+1 .
7.2
Het bepalen van eerste integralen
Om de eerste integralen van een differentiaalstelsel te bepalen volstaat het uiteraard om het stelsel op te lossen en dan stelling 7.1.2 toe te passen. Zoals we reeds ondervonden is dit niet altijd mogelijk, en daarom is het nuttig om alternatieve technieken te bestuderen om eerste integralen te bepalen. Ook hier bestaan er geen algemene methoden; in deze paragraaf geven we enkele elementaire methoden. Algemene natuurkundige wetten (bijvoorbeeld de wet van het behoud van energie of hoeveelheid beweging) leveren soms eerste integralen van stelsels differentiaalvergelijkingen die afkomstig zijn van fysische problemen. Soms kan men ook op andere manieren eerste integralen vinden. We herschrijven het stelsel (7.1) onder de vorm dyn dx dy1 = ··· = = (= k) (7.5) g g1 gn waarbij g, g1 , · · · , gn functies zijn van x, y1 , · · · , yn . Als een van de vergelijkingen van het stelsel slechts x, yi en y0i bevat, dan levert de algemene integraal van deze vergelijking een eerste integraal. Voorbeeld 7.2.1 We beschouwen het stelsel dy dz dx = = x−y x+y z De eerste vergelijking is een homogene differentiaalvergelijking: (x + y)dx = (x − y) dy We lossen deze op via de substitutie y = ux ; dy = udx + xdu We vinden achtereenvolgens (1 + u)dx = (1 − u)(udx + x du)
en tenslotte vinden we dat
dx (1 − u)du = x 1 + u2 1 ln |x| = bgtg u − ln (1 + u2 ) + ln c 2 p x2 + y2 =c exp bgtg xy
een eerste integraal is van (7.6).
101
(7.6)
Soms kan men ook op de volgende manier een eerste integraal vinden. Onderstel dat we functies α, α1 , · · · , αn van x, y1 , · · · , yn kan vinden zodat αg + α1 g1 + · · · + αn gn = 0 en αdx + α1 dy1 + · · · + αn dyn de totale differentiaal is van een functie h(x, y1 , · · · , yn ). Als (y1 , · · · , yn ) een oplossing is van het stelsel, dan geldt dus dh = αdx + α1 dy1 + · · · + αn dyn = k(αg + α1 g1 + · · · + αn gn ) = 0 zodat h(x, y1 , · · · , yn ) = c een eerste integraal is. We noemen deze methode ook wel de methode der multiplicatoren. De moeilijkheid is natuurlijk om de functies α, α1 , · · · , αn te vinden. We bekijken enkele voorbeelden. Voorbeeld 7.2.2 We beschouwen het stelsel dy dz dx = = y−z z−x x−y
(7.7)
Stel eerst α = α1 = α2 = 1. Dan is α(y − z) + α1 (z − x) + α2 (x − y) = 0 en bovendien is dx + dy + dz = d(x + y + z) een totale differentiaal. We hebben dus een eerste integraal x + y + z = c1 Stel nu α = x, α1 = y, α2 = z. Weer is α(y − z) + α1 (z − x) + α2 (x − y) = 0 en
1 xdx + ydy + zdz = d(x2 + y2 + z2 ) 2 is een totale differentiaal. Een tweede eerste integraal is dus x2 + y2 + z2 = c2 De oplossingen van (7.7) zijn dus de krommen x + y + z = c1 x2 + y2 + z2 = c2
(7.8)
De oplossingen zijn dus cirkels met middelpunt op de as R(1, 1, 1) en gelegen in een vlak loodrecht op deze as.
102
Voorbeeld 7.2.3 We beschouwen weer het stelsel (7.6). We proberen α = −xz, α1 = −yz, α2 = x2 + y2 Inderdaad is −xz(x − y) − yz(x + y) + z(x2 + y2 ) = 0 Jammer genoeg is z −xzdx − yzdy + (x2 + y2 )dz = − d(x2 + y2 ) + (x2 + y2 )dz 2 geen totale differentiaal. Alles zou ok zijn als de factor −1/2 er niet zou gestaan hebben. Merk op dat de factor 1 − d(x2 + y2 ) 2 voorkomt bij de differenti¨ering van (x2 + y2 )−1/2 , immers 1 d(x2 + y2 )−1/2 = − (x2 + y2 )−3/2 d(x2 + y2 ) 2 Daarom delen we α, α1 en α2 door (x2 + y2 )3/2 . We vinden dan (x2 + y2 )dz zd(x2 + y2 ) + αdx + α2 dy + α3 dz = − 2 2(x + y2 )3/2 (x2 + y2 )3/2 z = d p x2 + y2 zodat
z p = c2 2 x + y2
een tweede eerste integraal van het stelsel (7.6) levert. We kunnen dus concluderen dat de algemene integraal van het stelsel bestaat uit de krommen van de vorm p x2 + y2 = c1 exp bgtg xy z2 = c2 x2 + y2
103
Hoofdstuk 8 Lineaire parti¨ele differentiaalvergelijkingen van orde 1 8.1
Inleiding
Een parti¨ele differentiaalvergelijking is een betrekking tussen een functie y van n veranderlijken x1 , x2 , · · · , xn en de parti¨ele afgeleiden van y naar de xi . De orde van de parti¨ele differentiaalvergelijking is de hoogste orde van afleiden die voorkomt in de vergelijking. Een oplossing van de differentiaalvergelijking is een functie y = f (x1 , · · · , xn ) die identiek voldoet aan de vergelijking (tenminste op een zeker domein). Een probleem dat men zich kan stellen is het bepalen van de algemene integraal van een parti¨ele differentiaalvergelijking. Uit de volgende voorbeelden zal blijken dat dit niet altijd een erg relevant probleem is. Voorbeeld 8.1.1 z is een functie van de veranderlijken x en y. De oplossingen van de parti¨ele differentiaalvergelijking ∂z =0 ∂x wordt gegeven door de formule z = ψ(y) waarbij ψ een willekeurige functie van e´ e´ n veranderlijke is. Immers, als we z afleiden naar de veranderlijke x, dan vinden we 0. Men kan de willekeurige functie ψ vergelijken met de integratieconstante c bij een gewone differentiaalvergelijking. De algemene integraal van een parti¨ele differentiaalvergelijking van orde 2 hangt nu af van twee willekeurige functies. Dit blijkt uit het volgend voorbeeld. Voorbeeld 8.1.2 We nemen nu de volgende parti¨ele differentiaalvergelijking van orde 2: a2
∂2 z ∂2 z − =0 ∂x2 ∂y2
104
(8.1)
Deze parti¨ele differentiaalvergelijking wordt ook wel de e´ e´ ndimensionale golfvergelijking genoemd. We beweren nu dat z = φ(x − ay) + ψ(x + ay) een oplossing is van (8.1), ongeacht wat de (tweemaal differentieerbare) functies φ en ψ zijn. Immers, ∂z ∂x ∂z ∂y ∂2 z ∂x2 ∂2 z ∂y2
= φ0 (x − ay) + ψ0 (x + ay) = −aφ0 (x − ay) + aψ0 (x + ay) = φ00 (x − ay) + ψ00 (x + ay) = a2 φ00 (x − ay) + a2 ψ00 (x + ay)
en hieruit volgt onmiddellijk dat z een oplossing is van (8.1). Voor een differentiaalvergelijking van drie veranderlijken wordt de situatie nog ingewikkelder: de algemene integraal hangt nu af van een willekeurige functie van twee veranderlijken. Voorbeeld 8.1.3 u is een functie van drie veranderlijken x, y en z. De algemene integraal van de differentiaalvergelijking ∂u =0 ∂x is nu u = ψ(y, z) waarbij ψ een willekeurige functie van twee veranderlijken is. Uit deze voorbeelden blijkt dat een fysisch probleem dat zich laat vertalen in een parti¨ele differentiaalvergelijking nog niet opgelost is als we de algemene integraal van de parti¨ele differentiaalvergelijking kunnen bepalen: we moeten ook nog de onbekende functies bepalen die in de algemene integraal optreden (gebruik makende van rand- of beginvoorwaarden), en dit probleem is dikwijls moeilijker dan het bepalen van de algemene integraal. Daarom tracht men dikwijls de oplossing van het probleem rechtstreeks te vinden, zonder eerst de algemene integraal uit te rekenen. In dit hoofdstuk zullen we zien hoe men de algemene integraal van een lineaire parti¨ele differentiaalvergelijking van orde e´ e´ n kan bepalen. Notatie Zij z = f (x, y) een functie van twee veranderlijken. Men gebruikt dikwijls de volgende notaties: ∂z ∂z p= q= ∂x ∂y r=
∂2 z ∂x2
s=
∂2 z ∂x∂y
105
t=
∂2 z ∂y2
8.2
Homogene lineaire parti¨ele differentiaalvergelijking van orde e´ e´ n
Een parti¨ele differentiaalvergelijking van de vorm a1 (x1 , · · · , xn )
∂y ∂y + · · · + an (x1 , · · · , xn ) =0 ∂x1 ∂xn
(8.2)
noemt men homogene lineaire parti¨ele differentiaalvergelijking van orde e´ e´ n . Om deze te integreren beschouwt men het geassocieerd differentiaalstelsel. dx2 dxn dx1 = = ··· = a1 (x1 , · · · , xn ) a2 (x1 , · · · , xn ) an (x1 , · · · , xn )
(8.3)
Zoals we zullen zien is er een verband tussen de oplossingen van (8.2) (meetkundig gezien hyperoppervlakken in Rn+1 ) en de eerste integralen van (8.3). Onderstel dat f (x1 , x2 , · · · , xn ) = c een eerste integraal is van het geassocieerd stelsel (8.3). Op een oplossing van (8.3) is dus df =
∂f ∂f ∂f dxn = 0 dx1 + dx2 + · · · ∂x1 ∂x2 ∂xn
en, vanwege (8.3) ∂f ∂f + · · · + an (x1 , · · · , xn ) =0 (8.4) ∂x1 ∂xn Waar (8.3) voldoet aan de voorwaarden van de bestaansstelling gaat er door elk punt een oplossing van (8.3). In een dergelijk punt geldt dus (8.4). (8.4) drukt dus uit dat y = f (x1 , · · · , xn ) een oplossing is van (8.2). Elke eerste integraal van (8.3) levert dus een oplossing van (8.2). Omgekeerd, onderstel dat y = f (x1 , · · · , xn ) een oplossing is van (8.2). We zullen aantonen dat f (x1 , · · · , xn ) = c dan een eerste integraal is van (8.3). Uit het voorgaande hoofdstuk weten we dat we n − 1 verschillende eerste integralen van (8.3) kunnen bepalen. Beschouw zulk een stel verschillende eerste integralen. f1 (x1 , · · · , xn ) = c1 .. . f (x , · · · , x ) = c n n−1 1 n−1 a1 (x1 , · · · , xn )
We weten dan ook dat y = f1 (x1 , · · · , xn ), · · · , y = fn−1 (x1 , · · · , xn ) oplossingen zijn van (8.2). We hebben dus ∂ f1 ∂ f1 + · · · + an (x1 , · · · , xn ) =0 a1 (x1 , · · · , xn ) ∂x1 ∂xn .. . ∂ fn−1 ∂ fn−1 + · · · + a (x , · · · , x ) =0 a (x , · · · , x ) n n n 1 1 1 ∂x1 ∂xn ∂f ∂f a1 (x1 , · · · , xn ) =0 + · · · + an (x1 , · · · , xn ) ∂x1 ∂xn
106
We kunnen dit beschouwen als een lineair stelsel met a1 , · · · , an als onbekenden. Omdat de ai niet allen identiek gelijk aan nul zijn (anders is de parti¨ele differentiaalvergelijking gewoon 0 = 0) bezit het stelsel een van nul verschillende oplossing. Dus is de determinant van het stelsel gelijk aan nul. Deze determinant is de Jacobiaanse determinant van f1 , · · · , fn−1 , f ten opzichte van x1 , · · · , xn : ∂( f1 , · · · , fn−1 , f ) =0 ∂(x1 , x2 , · · · , xn ) Uit Stelling 7.4.2, Analyse I, volgt nu dat er een betrekking bestaat tussen f1 , · · · , fn−1 en f . Omdat de eerste integralen f1 = c1 , · · · , fn−1 = cn−1 verschillend zijn kan er geen betrekking bestaan tussen f1 , · · · , fn−1 . Daarom is noodzakelijk f = φ( f1 · · · , fn−1 ) en dus is f = c een eerste integraal van (8.3). Hiermee hebben we de volgende stelling bewezen: Stelling 8.2.1 y = f (x1 , · · · , xn ) is een oplossing van (8.2) als en alleen als f (x1 , · · · , xn ) = c een eerste integraal is van (8.3). Om de algemene integraal van (8.2) te vinden volstaat het om n − 1 verschillende eerste integralen f1 = c1 , · · · , fn−1 = cn−1 van (8.3) te bepalen. De algemene integraal is dan y = φ f1 (x1 , · · · , xn ), · · · , fn−1 (x1 , · · · , xn ) Voorbeeld 8.2.2 We beschouwen de parti¨ele differentiaalvergelijking y
∂u ∂u ∂u +x +z = 0 ∂x ∂y ∂z
(8.5)
Het geassocieerd stelsel is dx dy dz = = y x z We bepalen hiervan twee eerste integralen. De eerste vergelijking kan herschreven worden als xdx = ydy en levert de eerste integraal f1 (x, y, z) = x2 − y2 = c1 Om een tweede eerste integraal te vinden gebruiken we de methode der multiplicatoren. We zoeken functies α, β, γ zodat αy + βx + γz = 0 Stel α = β = z en γ = −(x + y) We zouden willen dat zd(x + y) − (x + y)dz een totale differentiaal is. Dit is het geval als we alles delen door (x + y)2 , immers zd(x + y) − (x + y)dz z = d − (x + y)2 x+y
107
De tweede eerste integraal is dus f2 (x, y, z) = −
z = c2 x+y
De algemene integraal van (8.5) is dus z 2 2 u = φ x −y , x+y waarbij φ een willekeurige functie in twee veranderlijken is.
8.3
De volledige lineaire vergelijking van orde e´ e´ n
Een parti¨ele differentiaalvergelijking van de vorm a1 (x1 , · · · , xn , y)
∂y ∂y + · · · + an (x1 , · · · , xn , y) = b(x1 , · · · , xn , y) ∂x1 ∂xn
(8.6)
noemen we een lineaire parti¨ele differentiaalvergelijking van orde 1. We kunnen deze herleiden tot een homogene vergelijking. We zoeken oplossingen in impliciete vorm, d.w.z. y wordt bepaald als functie van x1 , · · · , xn door een betrekking van de vorm ψ(x1 , · · · , xn , y) = 0
(8.7)
Dit betekent dat ψ nu de nieuwe onbekende functie is. Als we (8.7) afleiden naar xi , dan vinden we ∂ψ ∂ψ ∂y + =0 ∂xi ∂y ∂xi (zie ook het hoofdstuk over impliciete functies in “Analyse I”). Als we deze betrekking substitueren in (8.6), dan vinden we a1 (x1 , · · · , xn , y)
∂ψ ∂ψ ∂ψ + · · · + an (x1 , · · · , xn , y) + b(x1 , · · · , xn , y) =0 ∂x1 ∂xn ∂y
(8.8)
Dit is een homogene lineaire vergelijking in ψ. Het geassocieerd stelsel is dx1 dx2 dxn dy = = ··· = = a1 a2 an b
(8.9)
Dit is een stelsel van orde n, en het bezit n verschillende eerste integralen f1 = c1 , · · · , fn = cn . De algemene integraal van (8.8) is dus ψ = φ( f1 , f2 , · · · , fn ) waarbij φ een willekeurige functie in n veranderlijken is. De oplossing van (8.6) wordt impliciet gegeven door de betrekking φ( f1 , f2 , · · · , fn ) = 0 (8.10)
108
Voorbeeld 8.3.1 We beschouwen de parti¨ele differentiaalvergelijking yp + xq = z
(8.11)
We zoeken oplossingen van de vorm ψ(x, y, z) = 0 (8.8) wordt nu ∂ψ ∂ψ ∂ψ +x +z =0 ∂x ∂y ∂z en dit is juist (8.5). We hebben deze vergelijking reeds opgelost, (voorbeeld 8.2.2) en de algemene integraal is z ψ = φ x2 − y2 , x+y De algemene integraal van (8.11) wordt dus in impliciete vorm gegeven door de formule z 2 2 φ x −y , =0 x+y y
Met behulp van de stelling van de impliciete functies kunnen we deze oplossen naar bijvoorbeeld de tweede veranderlijke. We vinden z = ϕ(x2 − y2 ) x+y of z = (x + y)ϕ(x2 − y2 ) waarbij ϕ een willekeurige functie van e´ e´ n variabele is. Meetkundige interpretatie We beschouwen het geval n = 2. De vergelijking (8.6) kunnen we dan schrijven als a(x, y, z)p + b(x, y, z)q = c(x, y, z)
(8.12)
Elke oplossing van (8.12) stelt een oppervlak in R3 voor. We noemen zulk een oppervlak een integraaloppervlak. Neem een integraaloppervlak, en een punt (x, y, z) op dit oppervlak. De vergelijking van het raakvlak in dit punt is dan (x − x)p + (y − y)q = z − z (8.13) Hierbij is (x, y, z) een lopend punt in het raakvlak. Uit (8.12) en (8.13) volgt dat de rechte met vergelijking x−x y−y z−z = = (8.14) a b c in het raakvlak (8.13) ligt. Voor elk punt (x, y, z) ∈ R3 kunnen we de rechte met vergelijking (8.14) neerschrijven.
109
Een kromme die in elk punt raakt aan de rechte (8.14) door dat punt noemen we een karakteristieke kromme. Aangezien (dx, dy, dz) de componenten zijn van een vector aan de raaklijn van een kromme, worden de karakteristieke krommen bepaald door het stelsel dx dy dz = = a b c
(8.15)
en dit is juist het geassocieerd stelsel aan de parti¨ele differentiaalvergelijking (8.12). De karakteristieke krommen zijn dus de oplossingen van het geassocieerd differentiaalstelsel. Onderstel nu dat f1 (x, y, z) = c1 en f2 (x, y, z) = c2 twee verschillende eerste integralen zijn van het geassocieerd stelsel. Voor elke vaste waarde van c1 en c2 bepaalt het stelsel vergelijkingen f1 (x, y, z) = c1 f2 (x, y, z) = c2 een oplossing van het geassocieerd stelsel, en dus een karakteristieke kromme. Als we c1 en c2 over alle re¨ele getallen laten lopen, dan vinden we alle karakteristieke krommen. Beschouw nu de unie van alle karakteristieke krommen waarvoor de bijhorende constanten c1 en c2 voldoen aan een betrekking van de vorm φ(c1 , c2 ) = 0 De punten (x, y, z) die tot deze unie behoren voldoen aan de betrekking φ( f1 (x, y, z), f2 (x, y, z)) = 0 en dit is de vergelijking van een integraaloppervlak. We kunnen dit nog als volgt herformuleren: een integraaloppervlak ontstaat door van de verzameling karakteristieke krommen (die afhangt van de twee parameters c1 en c2 ) een deelfamilie te nemen die afhangt van slechts 1 parameter, en dan de unie van deze familie te nemen. Meer bepaald kunnen we besluiten dat elk integraaloppervlak beschreven wordt door karakteristieke krommen.
8.4
Het vraagstuk van Cauchy
Beschouw weer de parti¨ele differentiaalvergelijking (8.12) a(x, y, z)p + b(x, y, z)q = c(x, y, z) Gegeven is een kromme k in R3 . Gevraagd wordt om het integraaloppervlak van (8.12) te bepalen dat door deze kromme k gaat. Men kan dit beschouwen als het analogon van het vraagstuk bij gewone differentiaalvergelijkingen van orde 1 dat er in bestaat om een oplossing te vinden die door een gegeven punt gaat. In het geval van lineaire parti¨ele differentiaalvergelijking en van orde 1 is de situatie iets ingewikkelder: we hebben namelijk twee verschillende gevallen.
110
Eerste geval k is geen karakteristieke kromme. Men beschouwt alle karakteristieke krommen die door k gaan. Deze beschrijven een oppervlak dat door k gaat, en dat ook een integraaloppervlak is. Dit oppervlak is dus de gevraagde oplossing. Onderstel dat k gegeven wordt door de vergelijkingen f (x, y, z) = 0 g(x, y, z) = 0 De karakteristieke krommen worden gegeven door de eerste integralen f1 (x, y, z) = c1 f2 (x, y, z) = c2 Neem een punt (x, y, z) gelegen op de kromme k. Dit punt moet aan de vier vergelijkingen voldoen. Als we x, y en z uit de vier vergelijkingen elimineren, dan vinden we een betrekking tussen c1 en c2 , namelijk φ(c1 , c2 ) = 0 De c1 en c2 die behoren bij een karakteristieke vergelijking die het gevraagde oppervlak beschrijft moeten hieraan voldoen. De vergelijking van het gevraagde oppervlak is dus φ( f1 (x, y, z), f2 (x, y, z)) = 0
(8.16)
Tweede geval k is zelf een karakteristieke kromme. Iedere familie karakteristieke krommen afhankelijk van e´ e´ n parameter, die k bevat is een oplossing van het vraagstuk. In dit geval zijn er dus oneindig veel oplossingen van het vraagstuk. Om de vergelijking van een dergelijke oplossing neer te schrijven volstaat het om een willekeurige functie φ(u, v) van twee veranderlijken te beschouwen. Aangezien de kromme k een karakteristieke kromme is, kan deze geschreven worden onder de vorm f1 (x, y, z) = a f2 (x, y, z) = b voor zekere a, b ∈ R. De familie karakteristieke krommen met vergelijking φ( f1 , f2 ) − φ(a, b) = 0 bevat k en is dus een oplossing van het vraagstuk. Voorbeeld 8.4.1 We hernemen de parti¨ele differentiaalvergelijking yp + xq = z Welk integraaloppervlak gaat door de hyperbool met vergelijking y=0 xz = 1
111
We elimineren x, y en z tussen de vergelijkingen y=0 xz = 1 x2 − y2 = c1 z =c 2 x+y Dit levert
xz = 1 x2 = c1 z x = c2
of
x2 = c1 x2 = c12
en dus c1 c2 = 1 De vergelijking van het gevraagde oppervlak is dus z(x − y) = 1 Dit is een cilinder waarvan de beschrijvenden evenwijdig zijn met de rechte x=y z=0
112
Hoofdstuk 9 Variatierekening 9.1
De vergelijkingen van Euler
In volume 1 hebben we extreme waarden van functies van e´ e´ n of meer veranderlijken besproken. Een van de resultaten die we bewezen is het volgende: onderstel dat de scalaire functie f : Rn → R een (locaal) minimum of maximum bereikt in het punt ~a. Als f differentieerbaar is in ~a, dan geldt noodzakelijk dat grad f (~a) = ~0, of, equivalent ∂f ∂f ∂f (~a) = (~a) = · · · = (~a) = 0 ∂x1 ∂x2 ∂xn
(9.1)
(9.1) is dus een nodige voorwaarde opdat f een extremum zou bereiken. In dit hoofdstuk willen we dit resultaat uitbreiden tot de situatie waarbij de verzameling Rn vervangen wordt door een verzameling V die zelf een verzameling functies is. Om dit verder te verduidelijken vertrekken we van een klassiek meetkundig probleem. In R3 beschouwen we een oppervlak S met vergelijking z = z(x, y) Beschouw ook twee punten A en B op het oppervlak S. We stellen onszelf nu de vraag: langs welke (continue) kromme C gelegen op het oppervlak S die de punten A en B verbindt is de afstand van het punt A tot het punt B minimaal? Zulk een kromme noemt men een geodeet van het oppervlak S. Zoals we zien is de veranderlijke in dit probleem een kromme, en een kromme wordt in het algemeen beschreven door een stel functies. In ons geval kunnen we een kromme C gelegen op het oppervlak S die A en B verbindt beschrijven door een stel parametervergelijkingen x = x(t) y = y(t) C: z = z(x(t), y(t)) waarbij A = (x(t1 ), y(t1 ), z(x(t1 ), y(t1 )), B = (x(t2 ), y(t2 ), z(x(t2 ), y(t2 )) voor zekere waarden t1 en t2 van de parameter t. De verzameling V bestaat dus uit de functies (x, y) : [t1 ,t2 ] ⊂ R → R2
113
die continu zijn en die t1 en t2 afbeelden op respectievelijk de eerste twee componenten van A en B. Om ons probleem oplosbaar te maken zullen we bovendien onderstellen dat de functies x en y tweemaal continu differentieerbaar zijn over [t1 ,t2 ]. We zullen dit formeel noteren als volgt: x, y ∈ C2 [t1 ,t2 ] waarbij
Cm [t1 ,t2 ] = { f : [t1 ,t2 ] → R | f bezit continue afgeleiden tot op orde m op [t1 ,t2 ]} Ook de functie l : V → R kunnen we opschrijven, aangezien we formules hebben die toelaten om de booglengte te berekenen (zie volume 2). Voor de lengte l van de kromme C vinden we: Z t2 q x0 (t)2 + y0 (t)2 + z0 (t)2 dt l= t1
De lengte wordt dus gevonden door een functie van de kromme C en de afgeleiden van haar componenten te integreren van t1 tot t2 . Het algemeen probleem kunnen we formuleren als volgt: Probleem 9.1.1 (het algemeen probleem van de variatierekening) Gegeven is een functie f : R2n+1 → R, t1 ,t2 ∈ R en A = (a1 , a2 , · · · , an ), B = (b1 , b2 , · · · , bn ) ∈ Rn . Bepaal een kromme C in Rn met parametervergelijkingen x1 = x1 (t) x2 = x2 (t) C: .. . x = x (t) n
n
waarbij xi (t1 ) = ai , xi (t2 ) = bi zodanig dat de functionaal Z t2
I= t1
f (t, x1 , x2 , · · · , xn , x10 , x20 , · · · , xn0 )dt
extremaal wordt. Hiermee bedoelen we dat I voor de gezochte kromme minimaal of maximaal moet zijn onder alle naburige krommen die ook aan de hierboven opgesomde voorwaarden voldoen. We gaan een nodige voorwaarde opstellen waaraan de gezochte kromme C (of het gezochte stel functies x1 , x2 , · · · , xn ) moet voldoen. Vooraleer we dit kunnen doen moeten we eerst een lemma bewijzen. Lemma 9.1.2 Onderstel dat f : [a, b] → R continu, en dat Z b
f (x)φ(x)dx = 0 a
voor elke functie φ ∈ C2 [a, b] waarvoor φ(a) = φ(b) = 0. Dan geldt noodzakelijk dat f (x) = 0 voor elke x ∈ [a, b].
114
Bewijs. Onderstel dat er een c ∈ [a, b] bestaat zodat f (c) 6= 0. Omdat f continu is, bestaat er dus een omgeving van x waarop f (c) 6= 0. Hieruit volgt dat we c verschillend van a en b kunnen kiezen, zodat c ∈ (a, b). Er bestaat dan een interval (α, β) ⊂ (a, b) dat c bevat en waarop f hetzelfde teken heeft als f (c). We defini¨eren nu de functie φ : [a, b] → R als volgt: ( (t − α)3 (t − β)3 als t ∈ (α, β) φ(t) = 0 als t 6∈ (α, β) De functie φ is dan tweemaal continu differentieerbaar over [a, b] (bewijs dit zelf in de punten α en β). Ook is φ(a) = φ(b) = 0 zodat Z b
0=
Z β
f (x)φ(x)dx = a
f (x)(x − α)3 (x − β)3 dx
α
Deze laatste integraal kan echter niet nul zijn, aangezien de integrand een constant positief of negatief teken heeft over (α, β). Dus is er een tegenstrijdigheid, en dit bewijst dat f (x) = 0 voor elke x. We bekijken nu opnieuw het algemene probleem van de variatierekening. Om de notaties te vereenvoudigen kijken we naar het geval n = 3. We noteren x, y en z voor de veranderlijken, zodat we de extremalen moeten zoeken van de functionaal Z t2
I=
f (t, x, y, z, x0 , y0 , z0 )dt
t1
We zullen - zoals reeds eerder vermeld - enkel een nodige voorwaarde opstellen. Onderstel dus dat (x, y, z) een oplossing is van het probleem. Neem nu een (voorlopig vaste maar willekeurige) functie φ ∈ C2 [t1 ,t2 ] waarvoor geldt dat φ(t1 ) = φ(t2 ) = 0. Bekijk nu de krommen van de vorm x(t, α) = x(t) + αφ(t) y(t, β) = y(t) + βφ(t) z(t, γ) = z(t) + γφ(t) Deze verzameling krommen hangt af van de drie parameters α, β en γ. Omdat φ(t1 ) = φ(t2 ) = 0 gaan alle krommen door de punten A en B. We bekijken nu de functie Z t2
I(α, β, γ) =
f (t, x(t, α), y(t, β), z(t, γ), x0 (t, α), y0 (t, β), z0 (t, γ))dt
t1
I bereikt een extremum in (x, y, z) als (x, y, z) loopt over alle naburige krommen die A en B verbinden. Als we de klasse krommen beperken tot een deelklasse bereikt I natuurlijk nog een extremum. Daarom bereikt I(α, β, γ) een extremum voor α = β = γ = 0. Een nodige voorwaarde hiervoor is dat ∂I ∂I ∂I = = =0 ∂α ∂β ∂γ
115
in het punt (0, 0, 0). We berekenen nu deze parti¨ele afgeleiden met behulp van de regel van Leibniz. Hiervoor is het voldoende te onderstellen dat f en de parti¨ele afgeleiden van f continu zijn. We vinden Z t2 ∂f ∂I ∂f = φ(t) + 0 φ0 (t) dt ∂α ∂x ∂x t1 Z t2 Z t2 ∂f ∂f 0 = φ(t)dt + φ (t)dt 0 t1 ∂x t1 ∂x Z t2 h∂f it2 Z t2 ∂f d ∂f dt = φ(t)dt + φ(t) − φ(t) ∂x0 dt ∂x0 t1 t1 ∂x t1 (parti¨ele integratie) ! Z t2 ∂f d ∂f = φ(t) dt − ∂x dt ∂x0 t1 De ge¨ıntegreerde term viel weg omdat φ(t1 ) = φ(t2 ) = 0. Voor de extremaalkromme (x, y, z) geldt dus dat ! Z t2 d ∂f ∂I ∂f = − dt = 0 φ(t) ∂α ∂x dt ∂x0 t1 Omdat dit geldt voor elke φ ∈ C2 [t1 ,t2 ] waarvoor φ(t1 ) = φ(t2 ) = 0 vinden we, gebruik makend van lemma 9.1.2 dat d ∂f ∂f = ∂x dt ∂x0 Op analoge wijze vinden we dat d ∂f ∂f = ∂y dt ∂y0 d ∂f ∂f = ∂z dt ∂z0 Men noemt dit stelsel differentiaalvergelijkingen de vergelijkingen van Euler. De oplossingen noemt men de extremalen. De Eulervergelijkingen vormen een differentiaalstelsel van drie vergelijkingen van orde 2, en de algemene integraal hangt dus af van 6 constanten. Als we eisen dat de extremaal door de punten A en B gaat, dan krijgen we zes bijkomende voorwaarden die toelaten om de 6 constanten te bepalen. Onze resultaten kunnen veralgemeend worden tot het algemeen geval. Stelling 9.1.3 Een nodige voorwaarde opdat de functionaal Z t2
I= t1
f (t, x1 , x2 , · · · , xn , x10 , x20 , · · · , xn0 )dt
extreem wordt voor de kromme (x1 , · · · , xn ) is dat voldaan is aan het stelsel differentiaalvergelijkingen ∂f d ∂f = (9.2) ∂xi dt ∂xi0 voor i = 1, 2, · · · , n.
116
Net zoals in het driedimensionaal geval noemen we de vergelijkingen (9.2) de vergelijkingen van Euler. De oplossingen zijn de extremalen, en deze hangen af van 2n constanten. Deze constanten kan men bepalen door de eisen dat de kromme door de punten A en B gaat. Het geval n = 1 We bekijken nu het allereenvoudigste geval, namelijk dat waarbij de fucntionaal I van de volgende vorm is: Z t2 I= f (x, y, y0 )dx t1
De vergelijking van Euler wordt nu ∂f d ∂f = ∂y dx ∂y0
(9.3)
Dit is een differentiaalvergelijking van orde 2, die men in sommige gevallen kan herleiden tot een differentiaalvergelijking van orde 1. Eerste geval f hangt niet af van y, met andere woorden Z t2
I=
f (x, y0 )dx
t1
(9.3) wordt dan d ∂f =0 dx ∂y0 of
∂f =c ∂y0
Dit is een differentiaalvergelijking van orde 1. Tweede geval f hangt niet af van x, met andere woorden Z t2
I=
f (y, y0 )dx
t1
Nu vinden we
d ∂f ∂f ∂f ∂f d ∂f f − y0 0 = y0 + 0 y00 − y00 0 − y0 =0 dx ∂y ∂y ∂y ∂y dx ∂y0
en we kunnen besluiten dat f − y0
∂f =c ∂y0
en dit is weer een differentiaalvergelijking van orde 1.
117
y B A
x
0 Figuur 9.1: Geodeten in het vlak Voorbeeld 1: de geodeten in het vlak
Gegeven twee punten (x1 , y1 ) en (x2 , y2 ) in het vlak. Welk is de geodeet in het vlak (de kortste weg die de twee punten verbindt)? We parameteren de kromme in x. De functionaal is dus I=
Z x2 p
1 + y02 dx
x1
De integrand is niet afhankelijk van y en dus wordt de Eulervergelijking y0 p
1 + y02
=c
waaruit volgt dat y0 = c1 en y = c1 x + c2 Zoals verwacht is dit de vergelijking van een rechte. De gevraagde extremaalkromme is dus de rechte die A en B verbindt. Voorbeeld 2: het vraagstuk van de brachistochrone Gegeven zijn twee punten O en A in R3 , waarbij O hoger ligt dan A. Een massapunt met massa m wordt in O losgelaten, zonder beginsnelheid, en beweegt zonder wrijving langs een kromme die O en A verbindt onder invloed van de zwaartekracht. Welke vorm moeten we aan de kromme C geven opdat de tijd nodig om van O naar A te gaan minimaal is? We kiezen een rechthoekig assenstelsel met de oorsprong in O. De x-as is gericht volgens de verticale, en het xy-vlak gaat door A. Daar de beweging gebeurt zonder wrijving is op elk ogenblik
118
y
0 m A (h, a, o)
x Figuur 9.2: Het vraagstuk van de brachistochrone de kinetische energie gelijk aan de potenti¨ele energie, met andere woorden p mv2 = mgx of v = 2gx 2 ds Hierbij is g de valversnelling, en v = de snelheid van het deeltje. dt We parametreren de kromme in x. t en s worden dan functie van x en ds ds dt dt = =v dx dt dx dx s ds 1 + y02 + z02 = dt = √ dx 2gx 2gx De valtijd t is dus gelijk aan 1 ds =√ t= √ 2g C 2gx Z
Z hr 1 + y02 + z02
x
0
De Eulervergelijkingen zijn of
!
d dx
1 y0 √ p x 1 + y02 + z02
!
d dx
1 z0 √ p x 1 + y02 + z02
=0 =0
1 y0 p √ = c1 x 1 + y02 + z02 1 z0 = c2 √ p x 1 + y02 + z02
Hieruit volgt dat y0 c2 = z0 c1
119
dx
en c2 y + c3 z = c4 en dit is de vergelijking van een vlak evenwijdig met de x-as. Aangezien de gezochte kromme C in dit vlak ligt, en bovendien door de punten O en A gaat, volgt noodzakelijk dat dit vlak het xy-vlak is. Met andere woorden, z=0 We kunnen dit ook formeel aantonen door uit te drukken dat C door de punten O en A gaat: voor y = 0 is ook z = 0, en dus is c4 = 0. Voor y = a is ook z = 0, en dus is c2 a = 0 en c2 = 0. We vinden dus dat z = 0. We konden dit trouwens intu¨ıtief verwachten, maar het wordt hier dus formeel bevestigd. Er blijft nu slechts 1 Eulervergelijking over: y0 1 √ p = c1 x 1 + y02 y komt in deze differentiaalvergelijking niet voor. We lossen ze op door ze in parametervorm te herschrijven (zie Reeks 3 van de oefeningen over differentiaalvergelijkingen). We stellen y0 = tgt Dan is 1 + y02 = sec2 t en uit de vergelijking volgt dat 1 √ sint = c1 x Geparametreerd in t wordt de differentiaalvergelijking dus sin2 t x= 2 c1 0 y = tgt Nu is dy = tgt dx = tgt
dx dt dt
2 2 tgt sint costdt = 2 sin2 tdt 2 c1 c1 1 = 2 (1 − cos 2t)dt c1
=
Integratie levert y=
2t − sin 2t 2c21
120
We vervangen 2t door u en vinden volgend stel parametervergelijkingen voor de oplossing: 1 − cos u x= 2c21 u − sin u y= 2c21 Dit is de vergelijking van een cyclo¨ıde (zie volume 2).
9.2
Het vraagstuk met nevenvoorwaarden
We werken in R3 . Gevraagd wordt weer om de functionaal Z t2
I=
f (t, x, y, z, x0 y0 , z0 )dt
(9.4)
t1
extremaal te maken, maar nu eist men dat de beschouwde krommen op een gegeven oppervlak g(x, y, z) = 0
(9.5)
liggen. Een voor de hand liggende manier om dit vraagstuk aan te pakken is de volgende. Uit de stelling van de impliciete functies volgt dat, in de onderstelling dat grad g 6= 0, minstens e´ e´ n van de drie componenten x, y of z kan opgelost worden in functie van een van de twee overigen. We lossen (9.5) op naar e´ e´ n van de veranderlijken, bijvoorbeeld z. Dit levert als vergelijking voor het oppervlak z = h(x, y) (9.6) Afleiden naar t geeft, met behulp van de kettingregel: z0 =
∂h 0 ∂h 0 x+ y ∂x ∂y
Als we (9.6) en (9.7) substitueren in (9.4), dan vinden we Z t2 ∂h ∂h I= f t, x, y, h(x, y), x0 , y0 , x0 + y0 dt ∂x ∂y t1
(9.7)
(9.8)
We kunnen nu gewoon de Eulervergelijkingen oplossen voor (9.8). Het nadeel hieraan is dat de functie h soms moeilijk te bepalen is, of ingewikkeld wordt. We zullen daarom trachten om h uit de vergelijkingen van Euler te elimineren. Schrijf ∂h ∂h k(t, x, y, x0 , y0 ) = f t, x, y, h(x, y), x0 , y0 , x0 + y0 ∂x ∂y De vergelijkingen van Euler worden dan ∂k d ∂k − =0 ∂x dt ∂x0 ∂k d ∂k − =0 ∂y dt ∂y0
121
of
∂ f ∂ f ∂h ∂ f ∂z0 d ∂ f d ∂ f ∂h + + − − =0 ∂x ∂z ∂x ∂z0 ∂x0 dt ∂x0 dt ∂z0 ∂x ∂ f + ∂ f ∂h + ∂ f ∂z − d ∂ f − d ∂ f ∂h = 0 ∂y ∂z ∂y ∂z0 ∂y dt ∂y0 dt ∂z0 ∂y
(9.9)
We bekijken de eerste vergelijking van (9.9) van naderbij. Aangezien ∂ ∂h 0 ∂h 0 ∂2 h 0 ∂2 h 0 d ∂h ∂z0 = x + y = 2x + y = ∂x ∂x ∂x ∂y ∂x ∂x∂y dt ∂x wordt deze vergelijking ! ∂f d ∂ f ∂h ∂ f d ∂f + =0 − − ∂x dt ∂x0 ∂x ∂z dt ∂z0
(9.10)
Afleiden van de betrekking g(x, y, h(x, y)) = 0 naar x levert ∂g
∂h ∂x = − ∂g ∂x
(9.11)
∂z
Substitutie van (9.11) in (9.10) geeft nu ∂f d ∂f d ∂f ∂f − − ∂x dt ∂x0 = ∂z dt ∂z0 ∂g ∂g ∂x ∂z Merk op dat, zoals we wensten, de functie h(x, y) uit de vergelijking is verdwenen. Op volkomen analoge manier kunnen we uit de tweede vergelijking van (9.9) halen dat ∂f d ∂f ∂f d ∂f − − 0 ∂y dt ∂y0 = ∂z dt ∂z ∂g ∂g ∂y ∂z Alles samengevat vinden we ∂f d ∂f ∂f d ∂f d ∂f ∂f − − − ∂x dt ∂x0 = ∂y dt ∂y0 = ∂z dt ∂z0 ∂g ∂g ∂g ∂x ∂y ∂z
(9.12)
Dit zijn de vergelijkingen van Euler voor het variatieprobleem met nevenvoorwaarden. Samen met de vergelijking g(x, y, z) = 0 laten deze ons toe om de extremalen (x(t), y(t), z(t) te bepalen.
122
Voorbeeld: de geodeten op een oppervlak Gevraagd wordt de kromme (x(t), y(t), z(t)) op het oppervlak met vergelijking g(x, y, z) = 0 zo te bepalen dat de lengte Z q t2
I=
x0 (t)2 + y0 (t)2 + z0 (t)2 dt
t1
langs deze kromme minimaal wordt. De Eulervergelijkingen (9.12) worden in dit geval d x0 y0 d p p dt dt x0 (t)2 + y0 (t)2 + z0 (t)2 x0 (t)2 + y0 (t)2 + z0 (t)2 = − − ∂g ∂g ∂x ∂y d z0 p dt x0 (t)2 + y0 (t)2 + z0 (t)2 = − ∂g ∂z We kunnen dit differentiaalstelsel op een meer elegante wijze herschrijven als we de parameter t vervangen door de hoofdparameter s langs de kromme. Herhaal dat de hoofdparameter s de booglengte langs de kromme gemeten van een zeker vast punt is. We hebben dan dat ds2 = dx2 + dy2 + dz2 zodat
q dx dx ds dx = = x0 (t)2 + y0 (t)2 + z0 (t)2 dt ds dt ds en het differentiaalstelsel wordt d2 y d2 z d2 x ds2 = ds2 = ds2 ∂g ∂g ∂g ∂x ∂y ∂z We kunnen (9.13) ook nog als volgt herformuleren: de vector ∂g ∂g ∂g grad g = , , ∂x ∂y ∂z x0 =
is evenwijdig met de vector d2 x d2 y d2 z , , ds2 ds2 ds2 Wat is de betekenis van deze laatste vector? We weten dat de vector ~T = dx , dy , dz ds ds ds een eenheidsvector is die raakt aan de kromme. We zien ook dat d~T d2 x d2 y d2 z = , , ds ds2 ds2 ds2
123
(9.13)
Afleiden van k~T k = ~T · ~T = 1 naar s geeft ons d~T ~ ~ d~T ·T +T · =0 ds ds of
zodat
d~T ~ ·T = 0 ds d~T ~ ⊥T ds
d~T staat dus loodrecht op de kromme. De eenheidsvector evenwijdig met deze vector ds noteert men ~N en noemt men de hoofdnormaal langs de kromme. De lengte van deze vector noteert men κ en noemt men de kromming. Men heeft dus
De vector
d~T = κ~N ds Hoe groter κ, hoe sneller de kromme van richting verandert, en dit is de reden waarom we κ de kromming noemen. Het vlak door een punt van de kromme en met richtingsvectoren ~T en ~N is het vlak dat het dichtste bij de kromme ligt in een omgeving van het beschouwde punt. We noemen dit vlak het osculatievlak. We kunnen alles als volgt samenvatten: een geodeet op een oppervlak is een kromme waarvoor de hoofdnormaal ~N in elk punt evenwijdig is met de normaal op het oppervlak, grad g. Anders gezegd: een geodeet op een oppervlak is een kromme waarvoor het osculatievlak in elk punt loodrecht staat op het raakvlak aan het oppervlak. In sommige gevallen kan men de geodeten expliciet uitrekenen. De geodeten op de bol blijken de grote cirkels te zijn, terwijl de geodeten op de cilinder de schroeflijnen zijn. Voor details verwijzen we naar de oefeningen.
124
Deel III Numerieke Analyse
125
Hoofdstuk 10 Het oplossen van vergelijkingen en stelsels 10.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
(10.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 (10.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. (10.1) is dan equivalent met x = x + k(x) f (x) (10.2) De oplossingen van van (10.2) zijn de dekpunten van de operator O : [a, b] → R : x 7→ x + k(x) f (x)
126
We stellen als algoritme voor: kies x0 ∈ [a, b], en bereken xn+1 uit xn met behulp van de recursieformule xn+1 = O(xn ) = xn + k(xn ) f (xn ) (10.3) Als de rij (xn ) convergeert, dan is de limiet x zeker een oplossing van (10.2) (en (10.1)). Het volstaat in (10.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 (10.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 (10.1) het snijpunt van deze raaklijn en de x-as (zie Figuur 10.1). y
y = f(x) f(x)n
xn+1
xn
x
x
Figuur 10.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 )
(10.4)
(10.4) is het speciaal geval van (10.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. Zonder bewijs geven we volgende eigenschap. Stelling 10.1.1 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. Het nadeel is dat men bij elke iteratieslag de afgeleide f 0 (x) moet berekenen. Daarom gebruikt men soms de volgende vereenvoudiging.
127
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 )
(10.5)
(10.5) is weer een speciaal geval van (10.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 10.2). y y = f(x)
x xn–1
xn+1
x xn
Figuur 10.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 )
(10.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.
128
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 10.1.2 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 xn2 − a 1 a xn+1 = xn − = xn + 2xn 2 xn √ Men kan aantonen dat lim xn = a als x0 > 0. n→∞ √ We berekenen als voorbeeld 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
129
(10.7)
Startwaarde was x(1) = 2; na drie iteraties vinden we reeds x(4) = 1.73205080756888; alle decimalen zijn correct. Voorbeeld 10.1.3 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 1−
a =0 x2
De methode van Newton levert het volgende algoritme 1 (3axn − xn3 ) (10.8) 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 (10.8) √ 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 xn+1 =
1 (axn − xn3 ) 2a √ xn √ = ( a − xn )( a + xn ) 2a
xn+1 − xn =
√ 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 dus convergent. De limiet kan enkel a zijn. √ √en begrensd, Voor x ∈ ( a, 3a), is f 0 (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
130
√ % 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).
10.2
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 Definitie 10.2.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 10.2.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 m1 = kmα0 = kmα−1 0 m0 = cm0 en 1+α m2 = kmα1 = kcα mα0 = kcα mα−1 m0 0 m0 = c
en
2
2
m3 = kmα2 = kcα+α mα0 = c1+α+α m0
131
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 10.2.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 −
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.
132
√ 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 =c n→∞ mα n lim
Het bewijs maakt gebruik van de getallen van Fibonacci. 5) De volgende iteratieve methode is van orde 3: xn+1 = g(xn ) = xn −
f (xn ) f 00 (xn ) f (xn )2 − f 0 (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.
10.3
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
(10.9)
We zullen dit soms verkort opschrijven in vectorvorm. Stellen we f1 x1 f2 x2 . F = en X = .. .. . xn
fn dan kunnen we (10.9) herschrijven als F(X) = 0
(10.10)
We gaan tewerk als in § 10.1, en trachten een rij vectoren (Xm ) te construeren die naar de oplossing convergeert. (10.10) is equivalent met X = X + M(X)F(X)
(10.11)
indien M(X) een reguliere n × n-matrix is, voor elke waarde van X. Naar analogie met (10.2) stellen we volgende iteratieve formule voor Xm+1 = Xm + M(Xm )F(Xm )
133
(10.12)
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 (10.10) gevonden hebben. We benaderen nu de functies fi door hun Taylor ontwikkeling tot op orde 1: fi (x1 , . . . , xn ) ≈
(m) (m) fi (x1 , . . . , xn ) +
n
∂ fi
(m)
∑ ∂x j (x1
(m)
(m)
, . . . , xn )(x j − x j )
j=1
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 )
(10.13)
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 ). (10.13) wordt ∂ f2 ∂ f1 f1 (xm , ym ) (xm , ym ) − f2 (xm , ym ) (xm , ym ) ∂y ∂y xm+1 = xm − ∂ f1 ∂ f2 ∂ f1 ∂ f2 (xm , ym ) (xm , ym ) − (xm , ym ) (xm , ym ) ∂x ∂y ∂y ∂x (10.14) ∂ f ∂ f 1 2 f2 (xm , ym ) (xm , ym ) − f1 (xm , ym ) (xm , ym ) ∂x ∂x ym+1 = ym − ∂ f2 ∂ f1 ∂ f2 ∂ f1 (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 )
134
(10.15)
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 ∂ f2 0 (X) · · · 0 ∂x2 D(X) = .. .. .. . . . ∂ fn 0 0 ··· (X) ∂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
(10.16)
(10.16) is natuurlijk eenvoudiger dan (10.13), maar zoals we kunnen verwachten, is de convergentie van (10.16) minder goed dan die van (10.13). Voorbeeld 10.3.1 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+1 =
yn 2xn2 2 − + 2 9yn yn
135
We nemen als startwaarden bijvoorbeeld x0 = y0 = 2. Een eenvoudig MATLAB programma x(1) = 2; y(1) = 2; f or i = 1 : 20 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 10.3.2 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.
136
10.4
Stelsels lineaire vergelijkingen
In de cursus “Algebra en beschrijvende meetkunde”hebben we gezien hoe ineaire stelsels kunnen opgelost worden, bijvoorbeeld met de eliminatiemethode van Gauss of Gauss-Jordan, of, in het geval van een regulier stelsel, met de methode van Cramer. Deze methoden worden minder handig als de orde van het stelsel groot wordt (in de praktijk treden gemakkelijk 100 bij 100 of 1000 bij 1000 stelsels op). Zeker de methode van Cramer wordt dan irrelevant. Daarom ontwikkelt men iteratieve methoden om lineaire stelsels op te lossen. We bespreken slechts enkele eenvoudige voorbeelden. 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 (10.17) 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 c2 d2 a2 − xm − zm = b2 b2 b2 b3 d3 a3 = − xm − ym c3 c3 c3
xm+1 =
(10.18)
ym+1
(10.19)
zm+1
(10.20)
Merk op dat we de formules gemakkelijk terugvinden door (10.17) op te lossen naar x, y en z. 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 (10.19) xm beter vervangen door xm+1 . Om dezelfde reden vervangen we in (10.20) xm en ym door xm+1 en ym+1 . Dit levert de formules van GaussSeidel .
d1 b1 c1 − ym − zm a1 a1 a1 d2 a2 c2 = − xm+1 − zm b2 b2 b2 d3 a3 b3 = − xm+1 − ym+1 c3 c3 c3
xm+1 =
(10.21)
ym+1
(10.22)
zm+1
137
(10.23)
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 10.4.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. 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 =
(10.24)
αi
(10.25)
138
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 Opmerkingen 10.4.2 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 (10.24)(10.25) moet men dan vermenigvuldigen met de inverse van de matrix αi .
139
Hoofdstuk 11 Afleiden en integreren 11.1
De interpolatieformule van Lagrange
De oplossingen van de problemen in hoofdstuk 10 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 (11.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: 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
140
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 (11.2) Stel L j (x) = k(x − x0 )(x − x1 ) · · · (x − x j−1 )(x − x j+1 ) · · · (x − xn ) Dan is voldaan aan (11.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
(11.3)
De interpolatieveelterm van Lagrange We keren terug naar het algemeen probleem (11.1). Stel pn (x) = f0 L0 (x) + f1 L1 (x) + · · · + fn Ln (x)
(11.4)
pn voldoet aan (11.1) en wordt de interpolerende veelterm van Lagrange genoemd. In de volgende stelling bepalen we de restterm. Stelling 11.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)!
(11.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. 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 (11.5).
141
Voorbeeld 11.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.2) (1.15 − 1.0)(1.15 − 1.1) +0.88021 + 0.91031 (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)|
11.2
1.6 ≈ 0.000098 6
Numeriek Afleiden
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.
142
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 11.1.1) opschrijven: voor elke x ∈ I hebben we n
n
f (x) = ∑ fi Li (x) + ∏ (x − x j ) j=0
i=0
f (n+1) (ξ(x)) (n + 1)!
(11.6)
We leiden (11.6) af en vinden n
f 0 (x) = ∑ 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 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 0 (xk ) = ∑ fi Li0 (xk ) + (−1)n−k hn i=0
k!(n − k)! (n+1) f (ξk ) (n + 1)!
(11.7)
Praktische formules Na enig rekenwerk kunnen we (11.7) 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 (x − x0 )(x − x2 ) 1 L1 (x) = = − 2 (x − x0 )(x − x0 − 2h) (x1 − x0 )(x1 − x2 ) h (x − x0 )(x − x1 ) 1 L2 (x) = = 2 (x − x0 )(x − x0 − h) (x2 − x0 )(x2 − x1 ) 2h
L0 (x) =
143
Afleiden geeft 1 (2x − 2x0 − 3h) 2h2 1 L10 (x) = − 2 (2x − 2x0 − 2h) h 1 (2x − 2x0 − h) L20 (x) = 2h2 L00 (x) =
Substitutie in (11.7) geeft, na enig gereken 1 h2 (−3 f0 + 4 f1 − f2 ) + f (3) (ξ0 ) 2h 3 2 h 1 (− f0 + f2 ) − f (3) (ξ1 ) f 0 (x1 ) = 2h 6 1 h2 f 0 (x2 ) = ( f0 − 4 f1 + 3 f2 ) + f (3) (ξ2 ) 2h 3
f 0 (x0 ) =
(11.8) (11.9) (11.10)
Opmerkingen 11.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 (11.9), 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 (11.9) 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 (11.7), 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 (11.9) 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
1 h3 (−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
2
144
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
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
11.3
h2 3 3D
3
h3 4 12 D
3
4
4
4
h4 5 5D
Numeriek integreren
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
145
in de spilpunten x0 , x1 , . . . , xn . Net als in het § 11.2 vertrekken we van de interpolatieformule van Lagrange. n
f (x) = ∑ fi Li (x) + rn (x)
(11.11)
i=0
met rn (x) =
f (n+1) (ξ(x)) n (x − x j ) (n + 1)! ∏ j=0
(11.11) integreren geeft n
Z b
f (x)dx = ∑ fi
a
We stellen ci =
Rb a
Z b
i=0
a
Z b
Li (x)dx +
a
rn (x)dx
Li (x)dx, en vinden als integratieformule n
Z b a
f (x)dx = ∑ ci fi + E( f )
(11.12)
i=0
waarbij de gemaakte fout E( f ) gegeven wordt door de formule Z b
E( f ) = a
rn (x)dx
(11.13)
Opmerkingen 11.3.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. h
a = x0
x1
b x2
x n–1
xn
Figuur 11.1: De gesloten Newton-Cotes formules De integratieformules die 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 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
We zullen nu de co¨effici¨enten ci berekenen in een aantal eenvoudige gevallen.
146
(11.14)
a
x1
x0
xn
b
Figuur 11.2: De open Newton-Cotes formules Voorbeeld 11.3.2 De open formule met 1 spilpunt In dit geval is x0 = (a + b)/2, en L0 = 1. (11.14) 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 11.3.3 De trapeziumregel We schrijven de gesloten formule met twee spilpunten neer. Gebruik makend van (11.14) vinden we Z Z h 1 0 1 x1 tdt = (x − x1 )dx = − c0 = − h x0 h −h 2 en Z 1 x1 h c1 = (x − x0 )dx = h x0 2 zodat
Z x1 x0
h f (x)dx ≈ ( f0 + f1 ) 2
(11.15)
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 11.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 ) (11.16) 2 x0 De trapeziumregel wordt meestal in deze vorm gebruikt.
147
Voorbeeld 11.3.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 (11.14) volgt (−1)i ci = 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 We krijgen Z 2
Z x2
f (x)dx = h 0
x0
≈ 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
(11.17)
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
(11.18)
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
f (x)dx = x0
h h3 ( f0 + f1 ) − f 00 (ξ) 2 12
148
(11.19)
h h5 ( f0 + 4 f1 + f2 ) − f (4) (ξ) 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
Z x2
f (x)dx =
x0 Z x3 x0 Z x4 x0
(11.20) (11.21) (11.22)
Open formules Z x1 x−1 Z x2 x−1 Z x3 x−1 Z x4 x−1
h3 00 f (ξ) 3 3h3 00 3h ( f0 + f1 ) + f (ξ) f (x)dx = 2 4 4h 28h5 (4) f (x)dx = (2 f0 − f1 + 2 f2 ) + f (ξ) 3 90 95h5 (4) 5h (11 f0 + f1 + f2 + 11 f3 ) + f (ξ) f (x)dx = 24 144 f (x)dx = 2h f0 +
(11.23) (11.24) (11.25) (11.26)
(11.25) staat bekend als de formule van Milne. Voorbeeld 11.3.5 We berekenen het getal π met de formule π=
Z 2p
4 − x2 dx
0
1, 6
A B 0
C 1, 2
2
Figuur 11.4: Berekening van π We gebruiken de formule van Simpson met h = 0.2, en vinden π≈
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
149
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. 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
f (a + th)dt 0
a n
=
∑ ci fi + hchk+1
i=0
f (k+1) (ξ) (k + 1)!
en vandaar de volgende definitie: Definitie 11.3.6 Een integratieformule wordt een simplexformule genoemd indien men de rest kan schrijven onder de vorm f (k+1) (ξ) E( f ) = chk+2 (11.27) (k + 1)!
150
Zonder bewijs geven we volgende stelling Stelling 11.3.7 De Newton-Cotesformules zijn simplexformules. Voor een gegeven Newton-Cotes formule kunnen we nu gemakkelijk h en c bepalen. We geven twee voorbeelden. Voorbeeld 11.3.8 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 2 E(x ) = 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 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
(11.28)
Men ziet dit in door de substitutie x = a + ht toe te passen. Voorbeeld 11.3.9 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 4 1 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
(11.29)
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.
151
11.4
Differentiaalvergelijkingen
We beperken ons tot normale differentiaalvergelijkingen van orde 1 y0 = f (x, y)
(11.30)
y(x0 ) = y0
(11.31)
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. De methode van Euler-Heun We zullen (11.30) 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 . (11.30) en (11.31) zijn equivalent met Z x
y(x) = y0 +
f (x, y)dx
(11.32)
x0
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
(11.33)
(11.33) 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 y∗1 = y0 + h f (x0 , y0 ) (11.34) Dit betekent dat we de kromme y = y(x) vervangen door de raaklijn in het punt (x0 , y0 ). (11.34) noemt men de voorspellingsformule Men vervangt y∗1 in het rechterlid van (11.33) en voert een eerste iteratie uit h y1 = y0 + ( f (x0 , y0 ) + f (x1 , y∗1 )) 2
152
(11.35)
y1* y1
y0 x0 0
x1
1
Figuur 11.5: De methode van Euler-Heun (11.35) noemt men de verbeteringsformule. Meestal stopt men na e´ e´ n iteratie. Het heeft geen zin om de wortels van (11.33) nauwkeuriger te bepalen, aangezien (11.33) zelf reeds een benadering is, omdat men de rest in de trapeziumformule niet kan meenemen. Als we ons tevreden stellen met de voorspellingsformule (11.34), dan noemen we onze methode de methode van Euler. (11.34) gecombineerd met (11.35) 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. 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 baseren we ons op de methode van Simpson in plaats van de trapeziumregel. We beperken ons hier tot het algoritme, en laten de afleiding ervan achterwege. 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
(11.36)
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.
153
Runge-Kutta II
k1 = h f (x0 , y0 ) h k1 k2 = h f x0 + , y0 + 2 2 k2 h k3 = h f x0 + , y0 + 2 2 k4 = h f x0 + h, y0 + k3 1 y1 = y0 + (k1 + 2k2 + 2k3 + k4 ) 6
154
(11.37)
Index analytische functie 25, 87 Besselfunctie 97 binomiaalreeks 27 Cauchyrij 7 convergente rij 5 convergentiecriterium van Cauchy 8 convergentieinterval 18 convergentiestraal 18 convolutieproduct 44 criterium van Weierstrass 16 cyclo¨ıde 121 differentiaaloperator 77 differentiaalstelsel 67 differentiaalvergelijking 49 discretisatie 152 distributie 47 divergente rij 5 eigenwaarden en eigenvectoren 86 Eulervergelijkingen 116 exponenti¨ele groei 51 formule van Milne 149 formule van Newton 127 formule van Simpson 148 Fourierreeks 33 Fouriertransformatie 38 frequentie 31 geodeet 113 gesloten Newton-Cotes formule 146 gewoon punt 88 Gibbsverschijnsel 35 goniometrische reeks 30 homogeen lineair differentiaalstelsel 70 homogene differentiaalvergelijking 60 inwendig product 45 juiste differentiaalvergelijking 58 karakteristieke veelterm 76 karakteristieke vergelijking 77
Legendre, veeltermen van 91 Legendre, vergelijking van 89 lineair differentiaalstelsel 69 lineaire differentiaalvergelijking 63 lineaire differentiaalvergelijking van orde n 73 lineaire parti¨ele differentiaalvergelijking 106 linkervoorwaarde van Lipschitz 34 machtreeks 16 methode van de vaste richting 128 methode van Euler 153 methode van Gauss-Seidel 137 methode van Heun 153 methode van Runge-Kutta 153 Morrey 134 Newton-Raphson 134 norm van een functie 14 normale differentiaalvergelijking 49, 67 numerieke reeks 7 open Newton-Cotes formule 146 parti¨ele differentiaalvergelijking 104 positieve reeks 8 product van reeksen 22 pulsatie 31 puntsgewijze convergentie 14 puntsgewijze convergentie 9 rechtervoorwaarde van Lipschitz 34 regelmatig singulier punt 88 regula falsi 129 Runge-Kutta I 153 Runge-Kutta II 153 simplexformule 150 singulier 88 stelling van Abel 20 Taylorreeks 25 trapeziumregel 147 tridiagonale matrix 138 uniforme convergentie 10
155
uniforme convergentie 15 veralgemeende functie 47 verbeteringsformule 153 vergelijking van Bernoulli 66 voorspellingsformule 152 zifteigenschap 47
156
