8 Oefeningen Numerieke Analyse
8.a Fouriertransformatie
1. Een complexwaardige functie f op de re¨ele rechte en zijn Fouriergetransformeerde g voldoen aan de volgende relaties:
g(t) = Z ∞
−∞
f (x) eixt dx en f (x) = 1 2π
Z ∞
−∞
g(t) e−ixt dt . (8.1) Laat zien, dat de volgende symmetrie¨en gelden (f∗ is de complex toegevoegde van f ) :
f (x) = f∗(x) (f is re¨eel) g(−t) = g∗(t) f (x) = −f∗(x) (f is imaginair) g(−t) = −g∗(t)
f (x) = f (−x) (f is even) g(t) = g(−t) (g is even) f (x) = −f(−x) (f is oneven) g(t) = −g(−t) (g is oneven) f is re¨eel en even g is re¨eel en even
f is re¨eel en oneven g is imaginair en oneven f is imaginair en even g is imaginair en even f is imaginair en oneven g is re¨eel en oneven
Wat zijn de overeenkomstige symmetrie¨en bij de Discrete Fourier Transformatie (DFT) voor de vectoren x en y ∈ CN, x = (x0, x1, · · · xN −1)T en identificeer xN met x0 (idem voor y):
yk =
N −1X
j=0
xj e2πikj/N en xj = 1 N
N −1X
k=0
yk e−2πikj/N
voorbeeld: x = x∗ (x is re¨eel) impliceert yN −j = y∗j.
2. De klassieke FFT werkt voor re¨ele en complexe vectoren. In het geval van re¨ele vectoren kunnen zuiniger algorithmen gevonden worden.
a. Schrijf een algorithme die de Fouriergetransformeerde van twee re¨ele vectoren x en y ∈ IRN tegelijk berekent door gebruik te maken van de complexe vector x + iy .
b. Schrijf een algorithme die de Fouriergetransformeerde berekent van een enkele re¨ele vector.
Hint : Ga van een re¨ele transformatie op N punten over naar een complexe op N/2 punten.
Plaats de data met even index in het re¨ele, en de data met oneven index in het imaginaire gedeelte. Anders gezegd: doe de eerste reduktieslag N → N/2 apart in re¨ele arithmetiek.
3. Naast de (complexe) Fouriertransformatie
Hn =
N −1X
k=0
hk e2πikn/N,
bestaan er ook de (re¨ele) Fourier sinus- en cosinustransformaties, resp.
Fk =
N −1X
j=1
fjsin(πjk/N ) en Fk=
N −1X
j=0
fjcos(πjk/N ).
a. Maak van de re¨ele sinustransformatie op N − 1 punten een ( complexe ) FFT op 2N punten door bij een vector x ∈ IRN een geschikte vector y ∈ IR2N te defini¨eren
( gebruik sin x = 2i1(eix− e−ix) ).
b. Maak gebruik van de symmetrie y2N −j = −yj en van het feit dat yj re¨eel is om het proces effici¨enter te maken.
8.b Approximatie
a. Zij I een interval dat de k steunpunten t1, t2, · · · , tk−1, tk bevat. Zij f ∈ C2k(I) een functie op I. Zij πf(x) het Hermite-interpolatiepolynoom van graad 2k − 1, dat f(x) interpoleert in de gegeven steunpunten, d.w.z.
f (ti) = πf(ti) en f′(ti) = πf′(ti) voor i = 1, 2, · · · , k − 1, k . Bewijs dat er een ξx ∈ I bestaat zo dat
f (x) − πf(x) = Yk i=0
(x − ti)2 f2k(ξx) (2k)!
b. Benader f (x) := tg(x) door de lineaire interpolant ϕ(x) zo dat ϕ(π
4) = f (π
4) en ϕ(1.47112767) = f (1.47112767) = 10.000 .
Bepaal een (goede) bovengrens voor de afwijking |f(x) − ϕ(x)| op het interval [π4, 1.47112767].
Bepaal een benadering voor tg(π3) door ϕ(π3) te berekenen en bepaal ook
| tg(π
3) − ϕ(π 3) | .
c. Benader f (x) := cotg(x) door de lineaire interpolant ψ(x) zo dat ψ(π
4) = f (π
4) en ψ(1.47112767) = f (1.47112767) .
Bepaal een (goede) bovengrens voor de afwijking |f(x) − ψ(x)| op het interval [π4, 1.47112767]
bepaal een benadering voor cotg(π3) door ψ(π3) te berekenen.
Bepaal | f(π3) − ψ(π3) |.
Bepaal een benadering voor tg(π3) door ψ(1π
3) te berekenen en bepaal dan ook
| tg(π
3) − 1 ψ(π3) |
Waarom is deze laatste benadering van tg(3) veel beter dan die van 2) ?
d. Als de functie ϕ een compositie is van twee functies f en g, zodat ϕ(x) = g(f (x)) en als we over numeriek stabiele algoritmen voor de berekening van f (x) uit x en g(y) uit y, dan behoeft de compositie van beide algoritmen niet numeriek stabiel te zijn:
Laten Cf en Cg de conditiegetallen zijn van f resp. g (voor een relevante waarde van hun parameter), bewijs dan, dat Cg(1 + Cf)η een bovengrens is voor de absolute waarde van de relatieve afrondfout in de berekende waarde van g(f (x)).
Pas dit toe op de berekening van arctan(x) via de identiteit arctan(x) = arcsin( x
√1 + x2).
8.c Orthogonale polynomen
1. De rij Legendrepolynomen {Pn| n ∈ IN } kunnen we defini¨eren als de rij polynomen die ontstaat als we de rij functies 1, α1x, α2x2, · · · , αnxn, · · · met het Gram-Schmidt proces orthogonaliseren;
de co¨effici¨enten αn worden zo gekozen dat
Pn(1) = 1 ∀n (P.1)
De orthogonaliteit moet gelden t.o.v. het inprodukt (f, g) :=
Z 1
−1 f (t)g(t)dt (P.2)
a. Bewijs dat Pn loodrecht staat op alle polynomen van graad kleiner dan n.
b. Bewijs dat er co¨effici¨enten βn en γn bestaan zodat de polynomen voldoen aan de drieterms- recursierelatie
Pn+1(x) = βnxPn(x) + γnPn−1(x) . (P.3) c. Bewijs dat de functie
Qn(x) := dn
dxn(1 − x2)n
een polynoom is van graad n, dat loodrecht op alle polynomen van graad < n staat. Laat zien dat hieruit volgt
Pn(x) = (−1)n 2nn!
dn
dxn(1 − x2)n (dit is de z.g. Rodrigues relatie voor Pn ). (P.4) d. Bewijs dat
k Pnk2 = 2
2n + 1, βn= 2n + 1
n + 1 , γn= − n
n + 1 en Pn(x) = (2n)!
2n(n!)2xn+ l.o.t . (P.5) e. Als | x | ≤ 1 en als z voldoend klein is, dan kan de functie
V (x, z) := 1
√1 − 2xz + z2 (P.6)
ontwikkeld worden in een machtreeks rond z = 0, V (x, z) =
X∞ k=0
Rk(x)zk Bewijs dat Rk = Pk ∀k.
aanwijzing: Differentieer V en haar machtreeksontwikkeling partieel naar z en leidt hieruit een recurrente betrekking af voor Rk . V heet de ”voortbrengende functie” voor de Legendrepoly- nomen .
f. Laat zien dat V voldoet aan de parti¨ele differentiaalvergelijking
∂
∂x(1 − x2) ∂
∂xV + ∂
∂zz2 ∂
∂zV = 0 (P.7)
en dat hieruit volgt dat Pn voldoet aan de gewone differentiaalvergelijking
(1 − x2)y′′ − 2xy′ + n(n + 1)y = 0 (P.8)
g. Laat zien dat uit (P.8) opnieuw de orthogonaliteit van de Legendrepolynomen volgt.
h. Laat zien dat Pnprecies n onderling verschillende re¨ele nulpunten heeft in het open interval (-1, 1).
aanwijzing: stel dat Pn precies k tekenwisselingen in de punten x1, x2, · · · , xk in (-1, 1) heeft.
Bekijk dan het inprodukt ( Pn , Πki=1(x − xi) ) en leid hieruit een tegenspraak af voor k < n.
i. Laat zien dat tussen ieder tweetal nulpunten van Pn+1 een nulpunt van Pn ligt.
aanwijzing: Leid uit (8) een inhomogene eerste orde gewone differentiaalvergelijking af voor de
“Wronski-determinant” W ,
W := Pn Pn+1′ − Pn′ Pn+1,
en integreer deze (met variatie van constante ). Stel nu Pn+1(α) = Pn+1(β) = 0 en laat zien dat uit de gevonden uitdrukking voor W volgt dat Pn tussen α en β een tekenwisseling moet hebben.
j. Bewijs (uit de recurrente betrekking voor Pn) de somformule van Christoffel-Darboux:
Xn k=0
(2k + 1) Pk(x) Pk(y) = (n + 1)Pn+1(x) Pn(y) − Pn(x) Pn+1(y)
x − y (P.9)
Door de limiet y → x te nemen vinden we dat de Wronskiaan W strikt positief is (ga na).
2. Men definieert de Chebyshev-polynomen van de eerste soort Tn(x) op het interval [-1, 1] door Tn(cos(t)) = cos(nt) t ∈ [0, π], n ∈ IN (T.1) a. Toon aan dat Tn(x) een veelterm is van graad n, meerbepaald:
Tn(x) = 2n−1xn + termen van lagere orde ( n 6= 0 ), T0 = 1 . (T.2) b. Toon aan dat er co¨effici¨enten βn en γn bestaan zodat de polynomen voldoen aan de drieterms-
recursierelatie
Tn+1(x) = βn x Tn(x) + γn Tn−1(x) . (T.3) Bepaal βn en γn .
c. Toon aan dat Tn(x) en Tk(x) (n 6= k) orthogonaal zijn voor het inprodukt ( f, g) :=
Z 1
−1 f (t) g(t) 1
√1 − t2 dt . (T.4)
Aanwijzing. Transformeer in de integraal t = cos(τ ).
d. Bewijs dat Tn(x) loodrecht staat op alle veeltermen met graad strikt kleiner dan n.
e. Bepaal de norm k Tn(x) k (in de norm behorend bij inprodukt (T.4)).
f. Bewijs dat de functie
Qn(x) = p1 − x2 dn
dxn (1 − x2)n−12 (T.5)
een polynoom van graad n is, dat loodrecht op alle polynomen van graad strikt kleiner dan n staat, zodat Tn(x) = cn Qn(x) (d.w.z. Tn(x) is een veelvoud van Qn(x)).
g. Zoek de voortbrengende functie
V (x, z) :=
X∞ n=0
Tn(x) zn
aanwijzing: substitueer x = cos t = (eit + e−it)/2.
h. Zij Rn(x) = 2n−1xn + lagere orde termen een willekeurig polynoom van graad n met gelijke kopco¨effici¨ent als Tn, toon dan aan dat
x∈]−1,1[max | Tn(x) | = 1 en dat max
x∈]−1,1[ | Pn(x) | ≥ 1 . (T.6) Hieruit volgt, dat Tn van alle n-de graads polynomen met kopco¨effici¨ent 2n−1 het kleinste max- imum heeft.
3. Men definieert de Chebyshev-polynomen van de tweede soort Un(x) op het interval [-1, 1] door Un(cos(t)) := sin(n + 1)t
sin(t) t ∈ ]0, π], n ∈ IN . (U.1)
a. Toon aan dat Un(x) een veelterm is van graad n van de vorm
Un(x) = 2nxn + termen van lagere orde. (U.2) b. Toon aan dat er co¨effici¨enten βn en γn bestaan zodat de polynomen voldoen aan de drieterms-
recursierelatie
Un+1(x) = βn x Un(x) + γn Un−1(x) . (U.3) Bepaal βn en γn.
c. Toon aan dat Un(x) en Uk(x) voor n 6= k orthogonaal zijn voor het inprodukt ( f, g ) :=
Z 1
−1
f (t) g(t) p1 − t2 dt (U.4)
d. Bewijs dat Un(x) loodrecht staat op alle veeltermen met graad strikt kleiner dan n.
e. Bepaal k Un(x) k (in de norm behorend bij inprodukt (U.4)).
f. Bewijs dat de functie
Qn(x) := 1
√1 − x2 dn
dxn(1 − x2)n+12 (U.5) een polynoom van graad n is, dat loodrecht op alle polynomen van graad strikt kleiner dan n staat, zodat Un(x) = cnQn(x) en bepaal cn .
g. Zoek een voortbrengende functie
V (x, z) :=
X∞ n=0
Un(x) zn
8.d Een spline-benadering
Laat n ∈ IN een natuurlijk getal zijn en h := 1/n . Laat ∆ := {0, h, 2h, · · · , nh = 1} een verdeling zijn van [0,1] in n deelintervallen van gelijke lengte. Van een functie f ∈ C3([0, 1]) zijn de functiewaarden fk:= f (xk) (k = 0 · · · n + 1) gegeven in de punten
x0 := 0, xk:= kh −12h (k = 1 · · · n) en xn+1:= 1 .
We willen f in de punten {xk, k = 0 · · · n+1} interpoleren met een kwadratische spline π ∈ M2,1(∆) (stuksgewijs tweedegraads polynoom met continu¨ıteit van functie en afgeleide) waarvoor geldt:
π(xk) = fk. We kiezen voor de restriktie πkvan π tot het deelinterval [kh − h, kh] de representatie πk(x) := 1
2hαk (x − kh + h)2 − 1
2hαk−1 (kh − x)2 + βk . (1) a. Bepaal de dimensie van M2,1(∆) en laat zien dat representatie (1) de continu¨ıteit van π′ im-
pliceert.
b. Elimineer uit (1) de parameters βk (k = 1 · · · n) met de vergelijkingen π(xk) = fk (k = 1 · · · n) en leid vervolgens uit de continu¨ıteit van π in de punten kh (k = 1 · · · n − 1) vergelijkingen af voor de onbekende α’s van de vorm
αk−1+ 6αk+ αk+1 = 8
h(fk+1− fk) (k = 1 · · · n − 1). (2a) Leid uit π(0) = f0 en π(1) = fn+1 ook de sluitvergelijkingen af:
3α0+ α1 = 8
h(f1− f0) en αn−1+ 3αn= 8
h(fn+1− fn) (2b) c. Herschrijf het stelsel van n + 1 vergelijkingen (2) in matrix-vector vorm M a = r. Bewijs vervolgens, dat kMak∞≥ kak∞, ∀a ∈ IRn+1en dat dit stelsel vergelijkingen een unieke oplossing heeft.
d. We kunnen het stelsel M a = r oplossen met Gausseliminatie. Schrijf de algoritme voor deze speciale tridiagonale matrix op en laat zien dat er bij iedere slag van deze eliminatie de spil (het grootste element van de kolom) reeds op de hoofddiagonaal staat en dat er dus geen rijverwis- selingen nodig zijn. Hoeveel ”flops” heeft uw algoritme nodig? Is deze algoritme vectoriseerbaar?
Zo ja, leg uit hoe, zo nee, waarom niet?
e. Aangezien αk = π′(kh) zal deze f′(kh) moeten benaderen; bewijs dat deze benadering van orde O(h2) is. Ga hierbij alsvolgt te werk: definieer de vector d ∈ IRn+1 met componenten dk:= f′(kh), (k = 0 · · · n) en laat (met Taylorontwikkeling) zien,
kM(d − a)k∞= O(h2) . (3)
Laat zien, dat met (c) hieruit volgt kd − ak∞= O(h2). (Behandel de nulde en n-de component apart.)
f. Laat (met lineaire interpolatie in een deelinterval) zien, dat uit (e) volgt: kπ′ − f′k∞ = O(h2) (max-norm voor functies).
g. Laat (door integratie in een deelinterval) zien, dat uit (f) volgt: kπ − fk∞= O(h3).
h. Generalisatie. Laten de steunpunten {x0 := 0 < x1 < · · · < xn+1 := 1} willekeurig zijn en laat de verdeling ∆ := {t0:= 0, t1, · · · , tn:= 1} gegeven zijn met tk:= 12(xk+1+xk) (k = 1 · · · n−1), hk:= tk− tk−1 (k = 1 · · · n), h := max{|hk| : 1 ≤ k ≤ n} en
πk(x) := 1
2hkαk (x − tk−1)2 − 1
2hkαk−1 (tk− x)2 + βk . (4) Stel de matrix M op en laat zien, dat de afleiding van de conclusie van (g.) volledig analoog is aan het bewijs in het equidistante geval.
8.e Numerieke integratie
1. Laten { t1 , t2 , · · · , tn } de steunpunten en { w1 , w2 , · · · , wn } de gewichten zijn van de n-punts Gaussformule op [ −1 , 1 ] ; laten verder { Pn| n = 0, 1, 2, · · · , } de Legendrepolynomen zijn, waarvoor de drietermsrecursierelatie geldt:
( n + 1 ) Pn+1 (x) − ( 2n + 1 ) x Pn (x) + n Pn−1 (x) = 0 (n ≥ 0) (1) a. Laat zien dat ti een eigenwaarde is van de tridiagonale matrix
0 1 0 · · · 0 0 0
1
3 0 23 · · · 0 0 0
0 25 0 · · · 0 0 0
0 0 37 · · · 0 0 0
· · · · 0 0 0 · · · 0 2n−5n−2 0 0 0 0 · · · 2n−3n−2 0 2n−3n−1 0 0 0 · · · 0 2n−1n−1 0
en dat de vector ( P0(ti) , P1(ti) , · · · , Pn−1(ti) )T de bijbehorende eigenvector is.
b. Zij A = ( aij ) de matrix met elementen
aij :=
s2i − 1
2 wj Pi−1(tj) , 1 ≤ i, j ≤ n . Laat zien dat A orthogonaal is en dat hieruit volgt:
2 wj
= Xn i=1
(2i − 1) Pi−12 (tj) . (2)
c. Ga na dat de integrand in formule (3.12a) gelijk is aan een veelvoud van (x − xi)−2Pn2(x). Laat zien, dat de formule (3.13) voor het gewicht
wi= 2
(1 − t2i) [Pn′(ti)]2 (3) hieruit kan worden afgeleid door parti¨ele integratie, omdat (x − xi)−1Pn(x) Pn′(x) weer een polynoom van graad 2n − 2 is waarvoor de Gaussformule exact is.
d. Bewijs de relatie
(1 − x2) Pn′(x) = n Pn−1(x) − n x Pn(x) (4) en laat met de somformule (P.9) zien, dat hieruit de identiteit van de formules voor wi volgt.
2. Laten { t0 := −1 , t1 , · · · , tn := 1 } de steunpunten en { w0 , w1 , · · · , wn} de gewichten van een (n + 1)-punts integratieformule zijn, die exact is voor polynomen van graad ≤ n. Definieer voor n ≥ 1 het polynoom Qn van graad n − 1 door
Qn(x) := (−1)n+1 n(n + 1) 2n n! (1 − x2)
µ d dx
¶n−1
( 1 − x2 )n (1)
a. Toon aan dat
Qn(x) = −n(n + 1) 1 − x2
Z x
−1
Pn(t) dt (2)
waarbij Pn de Legendrepolynomen zijn.
b. Bewijs de orthogonaliteitsrelatie Z 1
−1
Qn(x) Qm(x) (1 − x2) dx = 0 als m 6= n (3)
c. Bewijs dat uit de differentiaalvergelijking van Legendre volgt Qn(x) = Pn′(x) = n(2n)!
2n(n!)2 xn−1 + lagere orde termen (4) d. Bewijs met behulp van (4)
Z 1
−1 (1 − x2) Q2n(x) dx = n(n + 1) Z 1
−1 Pn2(x) dx = 2n(n + 1)
2n + 1 (5)
e. Laat zien dat er voor ieder polynoom f ∈ IP2n−1 polynomen g en r ∈ IPn−2 bestaan zo dat f (x) = 1
2 f (−1) (1 − x) + 1
2 f (1) (1 + x) + (1 − x2) { Qn(x) g(x) + r(x) } (6) en leid hieruit af dat de steunpunten van de (n+1)-punts Lobattoformule gelijk zijn aan de nulpunten van (1 − x2) Pn′(x)
f. Laat zien dat deze steunpunten in het interval [ -1 , 1 ] liggen en dat de gewichten positief zijn.
g. Laat zien dat voor de restterm van de Lobattoformule geldt Z 1
−1 f (x) dx = Xn i=0
wi f (ti) + 2n + 2 n(2n + 1)
22n(n!)4
((2n)!)3 f(2n)(ξ) (7) h. Laat zien dat, analoog aan oefening 1.c de gewichten w0 en wn bepaald kunnen worden met de
polynomen (1 ± x) (Pn′(x))2 en de andere met (1 − x2) (Pn′(x))2/(x − tk)2.
3. Een integratieformule heet een Newton-Cotes formule als de steunpunten equidistant zijn.
Bepaal de gewichten van de Newton-Cotes formules op 4 en 5 punten en bepaal de resttermen ervan met behulp van de stelling van Peano. Merk op dat de restterm van de 4-punts formule groter is dan die van de driepuntsformule !
4. Voor functies met singulariteiten kunnen we speciale formules maken die de betreffende singu- lariteit exact integreren. Bepaal de gewichten w0 en w1 en het steunpunt t0 zo dat de formule
Z 1 0
√x f (x) dx = w0 f (t0) + w1 f (1) + rest
exact is voor polynomen van graad 0, 1 en 2 en bepaal de restterm (in de vorm C f(p)(ξ)).
5. Laat M (a; h)f een integratieformule op [ a , a + h ] zijn waarvoor geldt Z a+h
a
f (x) dx = M (a; h)f + d hp+1 f(p)(ξ), (ξ ∈ ( a , a + h ) )
Laat zien dat voor de gerepeteerde formule op [ a , b ] met h := b−aj geldt
Rh f :=
j−1X
i=0
M (a + ih; h)f = Z b
a f (x) dx + C hp + O(kp+1)
als f ∈ Cp+1 [ a , b ] waarbij C niet afhangt van h. Bewijs ook dat hiervoor de algemene extrapolatieformule geldt ( met stapgrootten h < k )
Rhf + hp
kp − hp(Rhf − Rkf ) = Z b
a f (x) dx + O(hp+1) 6. Bepaal de formules voor een Rombergschema op basis van de rij stapgrootten
h , 2h 3 , h
2 , h 3 , h
4 , h 6 , h
8 , h 12 , · · ·
7. Laat Th de gerepeteerde trapeziumregel op [a, b] zijn met stapgrootte h := (b − a)/n. Laat zien dat de formule die ontstaat door h2-extrapolatie op Th en Th/2 gelijk is aan de gerepeteerde Simpsonregel op [a, b] met stapgrootte h .
8. Laat zien dat voor de integraalbenaderingen Th(i) in het Rombergschema geldt : Th(i) =
Xp l=0
wl f (tl) met wl≥ 0 ∀ l (a)
Th(i) is exact voor polynomen van graad ≤ 2i + 1 (b) N.B. (a) impliceert dat het Rombergschema numeriek stabiel is.
9. Orthogonale polynomen en Gauss-integratie met een gewichtsfunctie: Laat w een niet-negatieve gewichtsfunctie zijn op (a, b) met eindige momenten (i.e. Rab xjw(x) dx < ∞ , ∀ j ) en defini¨eer de rij orthogonale polynomen Q0 , Q1 , · · · door
Q−1 := 0 , Q0 := 1 ,
Qk := xk +
k−1X
j=0
qkj xj zo, dat Z b
a Qk(x) xj w(x) dx = 0 als j < k a. Bewijs dat er getallen βk en γk (6= 0) bestaan zo dat
Qk+1 = (x + βk ) Qk − γk Qk−1
en dat de nulpunten van Qk gelijk zijn aan de eigenwaarden van de tridiagonale matrix
−β0
γ1 0
· · · 0 0
1
−β1
γ2
· · · 0 0
0 1
−β2
· · · 0 0
0 0 1
· · · 0 0
· · ·
· · ·
· · ·
· · ·
· · ·
· · · 0 0 0
· · · γk−2
0 0 0 0
· · ·
−βk−2
γk−1 0 0 0
· · · 1
−βk−1
Wat zijn de eigenvectoren ?
b. Bewijs dat Qk precies k enkelvoudige nulpunten heeft in (a , b).
c. Laten t1 , · · · , tk de nulpunten van Qk zijn en w1 , · · · , wk gewichten zo, dat Z b
a
xj w(x) dx = Xk i=1
wi tji j = 0 , 1 , · · · , k − 1 Bewijs dan dat voor alle polynomen f van graad ≤ 2k − 1 geldt
Z b
a
f (x) w(x) dx = Xk i=1
wi f (ti)
d. Bewijs dat er voor iedere f ∈ C2k [a, b] een ξ ∈ [ a , b ] bestaat zo dat Z b
a f (x) w(x) dx = Xk i=1
wi f (ti) + f(2k)(ξ) (2k)!
Z b
a Q2k(x) w(x) dx
e. Definieer di := Rab Q2i(x) w(x) dx en bewijs dat de vectoren vj ∈ IRk met componenten vij :=
s wi
dj−1 Qj−1 (ti) i , j = 1 , 2 , · · · k een orthonormaal stel vormen in IRk en dat hieruit volgt
1 wi =
k−1X
j=0
Q2j(ti) dj
10. Bepaal met de methode van de onbepaalde co¨effici¨enten de steunpunten, de gewichten en de co¨effici¨ent van de afbreekfout voor de volgende formules op het interval [−1 , 1]:
• de 3- en 4-punts-Gauss formules,
• de 3-, 4- en 5-punts-Lobatto formules.
Gebruik dat deze integratieformules veeltermen tot een zekere graad correct integreren. Maak ook gebruik van de symmetrie van het probleem.
11.
a. Bepaal de parameters w1, w2, t en de foutconstante C voor een tweepunts Radau integratiefor- mule van orde 3 van de vorm
Z a+h
a
f (x) dx = w1h f (a + th) + w2 h f (a + h) + Ch4f(3)(ξ)
en bewijs dat er voor iedere driemaal differentieerbare functie f een ξ ∈ (a, a + h) bestaat waarvoor deze formule korrekt is. (Aanwijzing: voor het bewijs kun je bijvoorbeeld gebruik maken van een geschikt interpolerend polynoom)
b. Laat ook zien dat voor de gerepeteerde tweepunts Radau formule geldt : Z b
a f (x)dx = h
n−1X
j=0
{w1f (a + jh + th) + w2 f (a + jh + h)} + Ch3(f′′(b) − f′′(a)) + o(h3) (h → 0)
n In 2log(In− In/2) 1 0.6699006
2 0.6815621 -6.42 4 0.6827792 -9.68 8 0.6829221 -12.77 16 0.6829395 -15.81
c. Voor de functie f (x) := (1 + x2) sin x werden bij gebruik van de gerepeteerde tweepunts Radau formule voor de integraal over het interval [0, 1] de volgende waarden verkregen: Komen deze uitkomsten overeen met wat U verwacht en kunt U op grond hiervan een betere benadering voor de integraal vinden? Motiveer Uw antwoord en berekeningen.
12. De benadering vanπ door Archimedes. Het oppervlak π van de eenheidsscirkel ligt ingeklemd tussen de oppervlakken In van de ingeschreven en Un van de omgeschreven regelmatige n-hoeken.
a. Ga na dat In en Ungegeven worden door In= 12n sin 2π
n en Un= n tan π n en dat verdubbeling van het aantal hoeken de volgende recursie geeft:
I2n=pInUn en 2 U2n = 1
I2n + 1 Un.
b. Ga na dat het oppervlak van de ingeschreven cirkel feitelijk een integraalbenadering geeft d.m.v. de (gerepeteerde) trapeziumregel en de omgeschreven cirkel een benadering door de mid- delpuntsregel, zodat Dn:= 13In+23Uneen benadering door de (veel nauwkeuriger) Simpsonregel is.
c. Bereken met de bovenstaande iteratie de waarde van Dk voor k = 4j, j = 1 · · · 6 uitgaande van I4 = 2 en U4 = 4 en verifieer, dat de relatieve fout in D128 al kleiner dan 10−7 is.
d. Verifieer, dat de fout in D2j ongeveer 161 maal de fout in Dj is en dat we dus het berekende resultaat kunnen verbeteren met h4-extrapolatie.
e. Verifieer ook, dat we met het Rombergschema uit de boven berekende zes waarden van Dj een benadering van π kunnen vinden met een relatieve fout kleiner dan de gangbare machineprecisie 2 × 10−16.
Opmerking: Ludolf van Ceulen (Antwerpen - Leiden 1540 – 1610) berekende 35 decimalen met de methode van Archimedes en VELE JAREN REKENWERK; ze zijn gepubliceerd op zijn grafsteen in de Leidse Pieterskerk, (http://www.math.rug.nl/˜top/pi-dag/5juli.html) een facsimile van deze steen is onlangs teruggeplaatst.
8.f Lineaire Algebra
1) Zij A ∈ IRn×n een reguliere matrix en u, v ∈ IRn kolomvectoren. Veronderstel vTA−1u 6= −1 . Bewijs de Sherman-Morrison formule :
(A + uvT)−1 = A−1−A−1uvTA−1 1 + vTA−1u .
2) Zij k · k een vectornorm op IRn en A ∈ IRn×n een matrix. Toon aan dat de afbeelding k · k : IRn×n 7→ IR+ : A → k A k met k A k := max
x6=0
k Ax k k x k , een matrixnorm is. Een dergelijke matrixnorm noemen we een geassocieerde norm.
3) Beschouw voor p ∈ [1, ∞) de vectornormen op IRn kxkpp :=
Xn i=1
|xi|p en kxk∞:= max
i |xi| . Bewijs voor de geassocieerde matrixnormen :
kAk1= max
j
Xn i=1
|aij| , kAk∞= max
i
Xn j=1
|aij| en kAk2 = max
x6=0,y6=0
|(Ax, y)|
kxk2kyk2
met (x, y) :=Pni=1xiyi.
4) Laat zien dat voor een matrixnorm k · k, geassocieerd aan een vectornorm, geldt dat
kABk ≤ kAkkBk ( multiplicatieve eigenschap) en kIk = 1 (I : eenheidsmatrix) . 5) Bewijs voor x ∈ IRn en A ∈ IRn×n de ongelijkheden:
1) kxk2 ≤ kxk1≤√
nkxk2 2) kxk∞≤ kxk2 ≤√ nkxk∞
3) 1
√nkAk2≤ kAk1 ≤√
nkAk2 4) 1
√nkAk∞≤ kAk2 ≤√
nkAk∞
Toon met voorbeelden aan dat de ongelijkheden scherp zijn, d.w.z. dat er gelijkheid kan optreden.
6) Toon aan dat de “Frobenius”-norm
k · kF : IRn×n → IR+ : A 7→
vu ut
Xn i,j=1
| aij |2
een matrixnorm is, maar dat deze niet geassocieerd is aan een vectornorm. Laat ook zien dat deze norm multiplicatief is en invariant onder unitaire transformaties (kU AkF = kAkF voor iedere unitaire afbeelding U ).
7) Toon voor de Frobeniusnorm aan:
k A k2F = spoor ( AT A ) = som der eigenwaarden van AT A k A k22 = grootste eigenwaarde van AT A
√1
n k A kF ≤ k A k2 ≤ k A kF
N.B. De wortels van de eigenwaarden van ATA zijn de singuliere waarden van A.
8) Het conditiegetal van de matrix A is κ(A) := kAkkA−1k. Bewijs:
i) κ(A) ≥ 1 voor alle multiplicatieve matrixnormen, ii) κ(αA) = κ(A) ∀α 6= 0,
iii) κ(A) = 1 als A unitair is en k · k geassocieerd is met de Euclidische vectornorm,
iv) κ(A) is gelijk aan de wortel uit het quotient van de grootste en de kleinste eigenwaarde van ATA , als k · k geassocieerd is met de Euclidische vectornorm.
9) Toon aan dat, als A en A + E inverteerbaar zijn, geldt:
k ( A + E )−1 − A−1 k ≤ k E k k A−1 k k ( A + E )−1 k . 10) Toon aan dat, als A inverteerbaar is en als kEk < kA−1k−1, dan geldt
k ( A + E )−1 k ≤ k A−1 k
1 − k A−1 k.k E k.
11) Zij k · k een multiplicatieve matrixnorm en zij A een matrix met kAk < 1. Bewijs dat I + A inverteerbaar is en dat geldt :
k I k
1 + k A k ≤ k ( I + A )−1 k ≤ k I k 1 − k A k .
12) Bewijs de volgende stelling : Veronderstel, dat A ∈ IRn×n een reguliere matrix is en dat b ∈ IRn een vector is en dat x ∈ IRn de oplossing is van Ax = b. Als ∆A ∈ IRn×n voldoet aan
k ∆A k k A−1 k = r < 1 , dan is het gestoorde stelsel
( A + ∆A ) y = b + ∆b ∆b ∈ IRn niet-singulier. Als er verder ook nog een δ ≥ 0 bestaat zodat
k ∆A k ≤ δ k A k en k ∆b k ≤ δ k b k , dan geldt
k x − y k
k x k ≤ 2 δ κ(A)
1 − r (mits x 6= 0) .
8.g Matrixalgoritmen
1. Zij voor gegeven vector u ∈ IRn, u 6= 0 (6= nulvector) de matrix Hu gedefinieerd door:
Hu := I − 2 uuT uTu met I de eenheidsmatrix.
a. Bepaal Hu (v) met v ⊥ u.
b. Bepaal Hu (ku) met k ∈ IR .
c. Toon aan dat deze transformatie Hu elke vector spiegelt t.o.v. u⊥, het orthogonale comple- ment van vect{u}.
d. Toon aan dat Hu een symmetrische orthogonale operator is.
e. Hoe moeten we u kiezen als we voor een gegeven vector x wensen dat Hu (x) ∈ vect{e1} .
f. Hoe moeten we u kiezen als we wensen dat
Hu (x) = k y ( k ∈ IR0 ; x , y 6= 0 ) Geef een algoritme om de matrix Hu te bepalen.
2. Zij y ∈ IRn en k ∈ IR, dan noemen we de matrices in IRn×n van de vorm N ( y , k ) := I + y · eTk
Gauss-Jordan transformaties.
a. In de veronderstelling dat N (y, k) inverteerbaar is, geef een formule voor zijn inverse.
b. Gegeven een vaste x ∈ IRn, onder welke voorwaarden kan men een y ∈ IRn vinden zodat N (y, k) (x) = ek.
c. Schrijf een algoritme om de matrix A te overschrijven met A−1 door gebruik te maken van Gauss-Jordan transformaties.
d. Welke voorwaarden op A verzekeren je dat de algoritme succesvol werkt?
3. Bepaal de LU-decompositie van
A =
µ7.00 9.00
6.00 8.00
¶
Reken in een 10-tallig getalstelsel met een mantisse van slechts drie beduidende cijfers. Er wordt niet afgerond maar afgekapt.
Bepaal H := A − ˜L ˜U met ˜L ˜U de berekende LU-decompositie.
4. Zij A een symmetrische matrix en zij J(k , k + 1 , ϑ) de Givens-rotatie rond de rijen k en k+1
J(k , k + 1 , ϑ) :=
1 0
· · · 0 0
· · · 0
0 1
· · · 0 0
· · · 0
· · ·
· · ·
· · ·
· · ·
· · ·
· · ·
· · · 0 0
· · · cos ϑ
− sin ϑ
· · · 0
0 0
· · · sin ϑ cos ϑ
· · · 0
· · ·
· · ·
· · ·
· · ·
· · ·
· · ·
· · · 0 0
· · · 0 0
· · · 1
a. Toon aan dat E = J(k , k + 1 , ϑ) A J(k , k + 1 , ϑ)T nog symmetrisch is.
b. Ga na welke matrixelementen van A onveranderd blijven (dus Aij = Eij) en bereken de andere elementen Eij die wel veranderen in functie van Aij , cos ϑ en sin ϑ .
c. Is het theoretisch en/of praktisch mogelijk Givens-rotaties Jl(kl, kl+ 1, ϑl) te vinden zodat ( Jp Jp−1· · · J3 J2 J1 ) A ( J1T J2T J3T · · · Jp−1T JpT )
een diagonaalmatrix wordt ?
d. Is het theoretisch en/of praktisch mogelijk Householder transformaties Hi te vinden zodat ( Hp Hp−1· · · H3 H2 H1 ) A ( H1T H2T H3T· · · Hp−1T HpT )
een diagonaalmatrix wordt ?
e. Wat gebeurt er als A niet symmetrisch maar antisymmetrisch is ?
5. Toon aan dat bij de QR-ontbinding van een matrix A met behulp van Givens-rotaties, het onmogelijk is dit in de volgende volgorde te doen:
for k:=1 to n-1 do for j:=0 to k-1 do ”maak A[n-j,k-j] nul ” 6. Bepaal de co¨effici¨enten a, b, c en d en bewijs per inductie:
Xn k=1
k = an2 + bn + c Xn
k=1
k2 = an3 + bn2 + cn + d Xn
k=1
Xk j=1
j = an3 + bn2 + cn + d Xn
k=0
( n − k )2 = an3 + bn2 + cn + d
7. Schrijf een effici¨ente routine om een lijn te trekken langs een aantal meetpunten ( xi , yi ), i = 1 · · · n in het vlak m.a.w. bepaal a en b zo datPni=1 ( axi + b − yi )2 minimaal is. Zoek hiervoor een matrix A en vectoren x en y zo dat het probleem zich herleidt tot
min
x k Ax − y k2
waarbij de kolommen van A liefst onderling orthogonaal zijn. Los dit op d.m.v. de nor- maalvergelijkingen.
8. Zij A een reguliere matrix en u, v vectoren. Geef een algoritme dat de ˜Q ˜R-ontbinding van A + u · vT berekent, uitgaande van de QR-ontbinding van A. Hierbij mogen i.p.v. (n2 − n)/2 Givens-rotaties slechts 2 ( n − 1 ) rotaties uitgevoerd worden (naast ook nog n2 flops ).
9. Bepaal de orthogonale matrices U en V en de diagonaalmatrix Σ van de singuliere waarde ontbinding U Σ V van de volgende matrices:
A =
a11 a21 a31
· · · an1
en B = ( b11 b12 b13· · · b1n )
10. Beschouw de matrices A =
1 10−4 10−4
10−4 1 10−4
10−4 10−4 2
B =
1 0.1 0.1
0.01 1 0
0 0.1
2
Geef een schatting van de eigenwaarden van deze matrices. Met welke nauwkeurigheid geldt deze schatting? Gebruik Gerschorin-cirkels.
8.h De Cholesky-ontbinding
Bij de algoritme van Gauss voor het oplossen van het stelsel vergelijkingen Ax = b wordt A ontbonden in een linksondermatrix L en een rechtsbovenmatrix U zo dat A = L U . Hiervoor zijn
2
3n3 + O(n2 ) flops nodig (een floating point operation of flop is een vermenigvuldiging, deling, optelling of aftrekking).
Als A symmetrisch en positief definiet is, kunnen we beter de CHOLESKY - ontbinding van A, A = LLT of de LDLT-ontbinding bepalen,
A = L D LT
waarin L een benedendriehoeksmatrix is met Lii = 1 (i = 1 · · · n) en D een diagonaalmatrix.
Dit kost maar 13 n3 + O(n2 ) flops. We zullen nu het bestaan van deze laatste ontbinding bewijzen en de algoritme voor de berekening ervan afleiden.
Definitie : Een matrix heet symmetrisch als AT = A Definitie : Een matrix heet positief definiet als
∀ x ∈ IRn : xT A x ≥ γ xT x voor zekere γ > 0 .
Neem voor het vervolg aan dat A een symmetrische en positief definiete n × n-matrix is.
opgave 1. Bewijs dat A inverteerbaar is en dat A−1 weer symmetrisch en positief definiet is.
We verdelen A nu in stukken als volgt:
A = (An) =
µ An−1 aT
a ann
¶
waarin An−1 de leidende (n − 1)×(n − 1) deelmatrix is en a een (n-1)-kolomvector is, bestaande uit de laatste kolom van A (zonder het laatste element).
opgave 2. Bewijs dat An−1 symmetrisch en positief definiet is.
opgave 3. Bewijs: ann − aT A−1n−1a > 0 Aanwijzing: beschouw xT A x met x =
µ A−1n−1a
−1
¶ .
Veronderstel dat L(n−1) een re¨ele (n − 1) × (n − 1) linksonderdriehoeksmatrix is met enen op de hoofddiagonaal,
L(n−1)ij = 0 als j > i en L(n−1)ii = 1 (i, j = 1 · · · n − 1 ) ,
veronderstel dat Dn−1 een (n − 1) × (n − 1) diagonaalmatrix is met positieve diagonaalelementen di (i = 1 · · · n − 1) en veronderstel dat dit de factoren zijn in de Cholesky-ontbinding van An−1 ,
An−1 = L(n−1)Dn−1L(n−1)T
opgave 4. Laat zien dat er een vector k ∈ IRn−1 en een positief getal dn bestaan zo dat A =
µ L(n−1) kT
0 1
¶ µ Dn−1 0T
0 dn
¶ µ L(n−1)T 0T
k 1
¶
en laat hieruit met induktie zien dat de Cholesky-ontbinding van A bestaat.
opgave 5. Schrijf de algoritme voor de bepaling van de Cholesky-ontbinding van A in (pseudo-) Pascal op basis van de hierboven afgeleide formules.
opgave 6. Nu het bestaan van de Cholesky-ontbinding bewezen is, is het mogelijk om ook andere volgorden voor de berekeningen te onderzoeken. Bepaal een volgorde, die beter ”vectoriseert” dan de boven gevonden algoritme. Neem hierbij als datastruktuur voor de matrix een lineair array van lengte 12n(n + 1), waarin de elementen van de matrix, die op of onder de hoofddiagonaal staan kolomsgewijs zijn opgeslagen. Schrijf de algoritme zo, dat de elementen van de matrix overschreven worden door de overeenkomstige niettriviale elementen van L en D. Probeer ook voor het oplossen van L D LT x = b een volgorde te vinden, die goed vectoriseert gegeven deze datastruktuur voor de matrix.
opgave 7. Laat A ∈ IRn×n een rijsgewijs diagonaaldominante matrix zijn, d.w.z.
als A = (aij )ni,j=1 dan | ajj| >
Xn i=1, i6=j
| aji |, ∀ j,
en laten alle diagonaalelementen 1 zijn, (aii = 1), bewijs dan, dat A een LU -ontbinding zonder rijverwisseling heeft, met k L k∞ ≤ n. Kunt U ook een (goede) bovengrens geven voor de norm van U ?
Aanwijzing: Laat zien dat de deelmatrix A(2 : n, 2 : n) na de eerste slag van de Gausseliminatie weer diagonaaldominant is en dat de niettriviale elementen in de bijbehorende Gausstransformatie in absolute waarde kleiner dan of gelijk aan 1 zijn.
8.i Givensrotaties etc
1. Een Givensrotatie J(k, l, ϑ) ∈ IRm×m met l > k is een lineaire transformatie, die de vectoren in het vlak opgespannen door ek en el roteert over een hoek ϑ:
J(k, l, ϑ)(αek + βel) = (αc + βs) ek + (−αs + βc) el met c := cos ϑ en s := sin ϑ , J(k, l, ϑ) x = x , ∀x ∈ vect{ek, el}⊥.
a. schrijf de matrixrepresentatie op van J(k, l, ϑ) voor l > k.
Laten c en s zo zijn, dat voor de l-de component van een zekere vector a ∈ IRm geldt:
(J(k, l, ϑ) a)l = 0 .
b. Laat zien, dat deze s en c (met c ≥ 0 en dus met −12π ≤ ϑ ≤ 12π) alsvolgt berekend kunnen worden:
if| ak| > | al| then t := al/ak; c := 1/√
1 + t2; s := t ∗ c else
t := ak/al; s := sign(t)/√
1 + t2; c := | t ∗ c | end
(8.2)
Laat ook zien, dat de afrondfout in de berekende waarde van c en s op deze manier kleiner is, dan bij gebruik van de formules:
˜c := | ak|/qa2k + a2l , ˜s := sign(ak) al/qa2k + a2l .
We kunnen het element alk(l 6= k) van een matrix A ∈ IRm×n op nul transformeren, door de matrix van links te vermenigvuldigen met een geschikte Givensrotatie J(k, l, ϑ). Als we daarna de resulterende matrix van links vermenigvuldigen met een Givensrotatie om een ander element nul te maken, kan het zijn dat het element alk weer een waarde ongelijk aan nul krijgt (ga na).
c. Geef aan, hoe de m − 1 elementen a2,1, a3,1, · · · , am,1 van de eerste kolom van de matrix A d.m.v. een serie van m − 1 Givensrotaties nul gemaakt kunnen worden. Laat ook zien, dat een QR-ontbinding van A ∈ IRm×n gemaakt kan worden met n(n − 1)/2 + n(m − n) Givensrotaties en bereken het aantal benodigde flops.
d. Laat zien, dat het i.h.a. niet mogelijk is om de n elementen am,1, am,2, · · · , am,n (m > n) van de onderste rij nul te maken met niet meer dan n Givensrotaties (toe te passen van links!).
Evenals bij een QR-ontbinding met Householdertransformaties is het onvoordelig om de matrix Q expliciet te berekenen, zowel gezien de extra hoeveelheid geheugen als het grotere aantal flops.
We willen Q dus in gefactoriseerde vorm (dus als een serie rotaties) bewaren, liefst zonder extra geheugenruimte. Volgens Stewart is het voldoende voor het reconstrueren van J(k, l, ϑ) om in het nulgemaakte matrixelement alk het getal ρ op te slaan:
if (c = 0) or (c > | s |) then ρ := s else ρ := sign(s)/c end (8.3) e. Laat zien, dat J(k, l, ϑ) hieruit gereconstrueerd kan worden op numeriek stabiele wijze (dat de
relatieve afrondfout in de gereconstrueerde c en s kleiner is dan 3η in absolute waarde).
f. In plaats van de normalisatie c ≥ 0 kunnen we ook de normalisatie c ≥ | s | or s > | c | nemen.
Bepaal de formules voor c, s en ρ, analoog aan (8.2) en (8.3) in dit geval.
2. De eerste stap in de berekening van de Singuliere-Waarden-Ontbinding van een matrix A is de transformatie van A naar een bidiagonaalmatrix B met behulp van Householder transformaties van links en van rechts. Een mogelijke algorithme voor een m × n-matrix met m ≥ n is:
for k := 1 to n − 2 do
Maak de elementen A[k + 1 : m , k] van de k-de kolom nul d.m.v. een geschikte Householder transformatie van links.
Maak de elementen A[k , k + 2 : n] van de k-de rij nul d.m.v. een geschikte Householdertransformatie van rechts.
end ;
for k := n − 1 to n do
Maak de elementen A[k + 1 : m , k] van de k-de kolom nul d.m.v. een geschikte Householdertransformatie van links.
end
a. Laat zien, dat na afloop van deze algorithme alleen de matrixelementen A[j, j] en A[j, j + 1]
(1 ≤ n) nog waarden ongelijk nul kunnen bevatten.
b. Bepaal het aantal flops, benodigd voor bovenstaande algorithme.
c. Als m veel groter is dan n is het voordeliger om eerst een QR-ontbinding A = QR te maken en dan de bidiagonalisatie-algorithme toe te passen op het bovenste n × n-blok van R (de rest bevat toch alleen maar nullen). Beschrijf deze alternatieve algorithme en bepaal het aantal flops, dat ervoor nodig is. Wanneer kost deze variant minder werk dan (1)?
d. In het alternatief van vraag 3 kunnen we de bidiagonalistie van het niet-triviale n × n-blok van R doen met Householdertransformaties als in (1). Dit blok heeft een speciale structuur, het is zelf een bovendriehoeksmatrix. Laat zien, dat we de bidiagonalisatie hiervoor ook kunnen doen met behulp van een geschikte rij van 12(n − 1)(n − 2) paren Givensrotaties. Zo’n paar bestaat uit een Givensrotatie van rechts, die ´e´en element boven de diagonaal nul maakt en ´e´en element ongelijk onder de diagonaal cre¨eert, gevolgd door een Givensrotatie van links, die dit element onder de diagonaal weer nul maakt. Beschrijf deze algorithme en bepaal het aantal f lops, dat ervoor nodig is.
3. De numerieke stabiliteit van een QR-ontbinding met Householdertransformaties of met MGS kan iets verbeterd worden door (analoog aan rijverwisselingen bij Gausseliminatie) in iedere slag de kolom met grootste lengte vooraan te plaatsen. De algorithme voor de QR-ontbinding met Householdertransformaties wordt dan:
for k := 1 to n − 2 do
Zoek onder de kolommen A[k · · · m , j] met k ≤ j ≤ n de kolom met grootste lengte en verwissel deze met de kolom A[k · · · m , k].
Maak de elementen A[k + 1 · · · m , k] nul d.m.v. een geschikte Householder transformatie van links.
end
Als we dit botweg zo implementeren en in iedere slag de normen van alle resterende kolommen moeten bepalen, betekent dit een toename van het aantal flops met 50% (ga na!). Als we echter vooraf (de kwadraten van) alle kolomnormen in een array opslaan en dit array in iedere slag, tegelijk met de matrix bijwerken, kost het niet m´e´er dan n − k flops extra in de k-de slag.
a. Beschrijf deze algorithme met kolomverwisselingen.
b. Laat zien, dat een analoge strategie mogelijk is bij MGS en beschrijf ook deze algorithme.
4. Laat J een re¨ele m × n-matrix zijn met m ≥ n en r ∈ IRm een (kolom)vector. Defini¨eer de uitgebreide (m + n) × n-matrix A en (m + n)-vector b door
A :=
µ J µ In
¶
en b :=
µr 0
¶ , waar µ een re¨eel getal is (6= 0) en Inde identiteit in IRn en defini¨eer
Dµ := ( JT J + µ2 In )−1 JT .
a. Laat zien, dat x := Dµr de oplossing levert voor het het kleinste-kwadratenprobleem min
x ∈IRn k Ax − b k2 . b. Toon aan dat
k Dµ − J† k2 = µ2 σn (σ2n + µ2 )
waarbij J† de pseudo-inverse matrix van J is, rang J = n en σn de kleinste singuliere waarde van J is.
aanwijzing: Bepaal de singuliere-waardenontbinding van Dµ uit die van J.
c. Stel dat, de QR-ontbinding van J gegeven (reeds berekend) is, J = ˜Q ˜R met ˜Q ∈ IRm×m orthogonaal en ˜R ∈ IRm×n bovendriehoeks. Beschrijf hoe de QR-ontbinding van A berekend kan worden met n Householdertransformaties die samen slechts 23n3 + O(n2) flops vergen.
aanwijzing: laat eerst zien,
A = µQ˜
0 0 In
¶ µ R˜ µ In
¶
5. Laten u1 , u2 , v1 en v2 ∈ IRn kolomvectoren zijn met u1 ⊥ u2 en v1 ⊥ v2 . Bepaal de singuliere-waardeontbinding van
a. A = u1 vT1
b. B = u1 vT1 + u2 vT2
6. Laat A een n×n matrix zijn, waarvan de QR-ontbinding reeds berekend is, A = QR. Construeer een effici¨ent algorithme, die voor gegeven vectoren u en v de QR-ontbinding berekent van A + uvT (een z.g. rank-1 update van A). De algorithme mag niet meer dan O(n2) flops vragen.
8.j Niet-lineaire problemen
1. Laat f een voldoend gladde functie zijn op IR met een nulpunt α en met f′(α) 6= 0.
a. Laten a en b twee punten (in de buurt van α) zijn en laat β het nulpunt zijn van het lineaire interpolatiepolynoom van f op de steunpunten a en b. Laat zien, dat geldt:
β = af (b) − bf(a)
f (b) − f(a) (1)
en dat er ξ en η ∈ int(a, b) bestaan zo, dat β − α = 12(α − a)(α − b) b − a
f (b) − f(a)f′′(ξ) = 12(α − a)(α − b)f′′(ξ)
f′(η ) . (2) b. Op grond van (1) kunnen we het volgende iteratieve proces (de koordenmethode of koordennewton
of secant-methode) defini¨eren:
kies startpunten x0 en x1 ;
for k = 1, 2, · · · do xk+1 := xk−1f (xk) − xkf (xk−1)
f (xk) − f(xk−1) . (3) Bewijs dat dit proces lokaal convergent is door te bewijzen met (2):
|xk+1− α| ≤ C|xk− α||xk−1− α| (4)
voor zekere C > 0 en bepaal hieruit de orde van dit proces. Aanwijzing: laat εk := xk− α de afwijking zijn van de k-de iterand. Deze voldoet dan aan de relatie εk ≈ Kεpk−1 geldend voor een proces van orde p en aan de relatie εk ≈ Cεk−1εk−2 geldend voor grote k op grond van (2).
Leid uit deze beide relaties een vergelijking voor p af.
2. Gegeven is de functie f (x) := x3− x2+ x − 2.
a. Laat zien, dat f precies ´e´en re¨eel nulpunt heeft. Noem dit punt α.
b. Voor de bepaling van α kunnen we uit de vergelijking f (x) = 0 de volgende successieve-substitu- tieprocessen afleiden:
(a) xn+1:= 2 + x2n− x3n en (b) xn+1:= 1 + 1 x2n+ 1. Bewijs, dat (a) niet en (b) wel lokaal convergent is (naar α).
c. Bewijs, dat (b) voor iedere startwaarde x0 ∈ IR convergeert naar α (dus globaal convergent is).
3. Laat ϕ tweemaal continu differentieerbaar zijn en laat α een punt zijn met α = ϕ(α) en ϕ′(α) = 0. Bewijs dan dat het proces xn+1 := ϕ(xn) kwadratisch (van orde 2) convergeert in de buurt van α.
