TECHNISCHE UNIVERSITEIT EINDHOVEN Faculteit Wiskunde en Informatica Tentamen Numerieke Methoden voor Werktuigbouwkunde (2N460) op donderdag 24 juni 2010, 14.00-17.00 uur. De uitwerkingen van de opgaven dienen duidelijk geformuleerd en beargumenteerd te worden.
1. Beschouw de functie f gedefinieerd door p f (x) = 2(1 − cos x), 0 ≤ x < π. We veronderstellen dat we y = f (x) exact kunnen berekenen. (a) Zij y = f π2 en yˆ = f π2 + ∆x met |∆x| ≤ ε. Toon aan dat in eerste orde √ benadering |ˆ y − y| ≤ ε/ 2. (b) Bepaal de eerste orde benadering c(x) van het conditiegetal voor het probleem y = f (x). Toon aan dat dit probleem goed geconditioneerd is voor x → 0. Voor de berekening van y = f (x) (10−8 ≤ x ≤ 10−6 ) hebben we de beschikking over het volgende M ATLAB-script: clear all; format long e; k = 6:0.001:8; x = 10.ˆ(-k); y = sqrt( 2*(1-cos(x)) ); loglog( x,abs(y),’k’ ); xlabel( ’x’,’FontSize’,16 ); ylabel( ’y’,’FontSize’,16 ); grid; set( gca,’FontSize’,16 ); Het resultaat staat in de figuur onderaan de bladzijde. (c) Verklaar het verloop van de grafiek. (d) Geef een betere formule voor de berekening van y = f (x) voor 0 < x 1. Motiveer uw antwoord. −5
10
−6
y
10
−7
10
−8
10 −8 10
−7
10 x
1
−6
10
Tentamen Numerieke Methoden voor Werktuigbouwkunde (2N460) op donderdag 24 juni 2010, 14.00-17.00 uur.
2. De variabele u = u(x, y) voldoet aan het randwaardeprobleem −∇2 u = f (u), (x, y) ∈ Ω = (0, 1) × (0, 1), u(x, y) = 0,
(x, y) ∈ ∂Ω,
(1a) (1b)
met f (u) een nader te specificeren functie. We willen van dit probleem een numerieke oplossing bepalen m.b.v. centrale differenties. Daartoe introduceren we het rooster (xi , yj ) = ((i − 1)h, (j − 1)h),
i, j = 1, 2, . . . N,
met h = 1/(N − 1) de maaswijdte in de x- en y-richting. De numerieke benadering van u(xi , yj ) geven we aan met ui,j . (a) Geef het centrale differentieschema voor de parti¨ele differentiaalvergelijking (1a) in een inwendig roosterpunt. We kunnen de differentieschema’s voor alle inwendige roosterpunten samenvoegen tot het algebra¨ısche stelsel Au = f (u),
(2)
waarbij de vector u alle onbekenden ui,j bevat. Voor de functie f (u) in (1a) kiezen we eerst f (u) = 1. (b) Wat is een geschikte oplosmethode voor (2)? Motiveer uw antwoord. In het vervolg van deze opgave beschouwen we de Gauss-Seidel methode. We noteren de (k) componenten van de kde iterand als ui,j . (c) Formuleer de Gauss-Seidel methode voor het differentieschema afgeleid in onderdeel (a). Geef hierbij ook de ordening van de roosterpunten aan. Geef een goed stopcriterium. Motiveer uw antwoord. Tenslotte nemen we f (u) = cos u. (d) Formuleer de Gauss-Seidel methode voor het differentieschema afgeleid in onder(k+1) deel (a). Beschrijf een strategie voor het berekenen van ui,j
2
Tentamen Numerieke Methoden voor Werktuigbouwkunde (2N460) op donderdag 24 juni 2010, 14.00-17.00 uur.
3. In deze opgave willen we de afgeleide f 0 (x) van een functie f benaderen m.b.v. de achterwaartse differentieformule 1 1 3f (x) − 4f (x − h) + f (x − 2h) + h2 f 000 (x) + O(h3 ) 2h 3 1 2 000 3 = Dh [f ](x) + h f (x) + O(h ). 3
f 0 (x) =
(3)
(a) Geef een afleiding van bovenstaande differentieformule m.b.v. Taylorreeksen. In het vervolg van deze opgave nemen we f (x) = ex . Voor de berekening van f 0 (1) hebben we de beschikking over het onderstaande M ATLAB-script. clear all; format long e; k = linspace( 1,40,40 ); h = 2.ˆ(-k); fd = ( 3*exp(1)-4*exp(1-h)+exp(1-2*h) )./(2*h); err = abs( fd-exp(1) ); loglog( h,err,’k’ ); grid; xlabel(’h’); ylabel(’error’); (b) Voer dit script uit, schets de grafiek en verklaar de resultaten. We kunnen m.b.v. Richardson extrapolatie een schatting voor de fout in Dh [f ](x) bepalen, en daarmee ook een betere benadering Dh∗ [f ](x) voor f 0 (x). We nemen nu h = 2−k (k = 1, 2, . . . , 20). (c) Bepaal een schatting voor de fout in Dh [f ](1) als functie van h. Bepaal ook een betere benadering Dh∗ [f ](1). Wat is de orde van nauwkeurigheid van Dh∗ [f ](1)? (d) Schrijf een M ATLAB-script welke de fout in zowel Dh [f ](1) als Dh∗ [f ](1) tekent als functie van h. Voer dit script uit, schets de grafiek en verklaar de resultaten.
3
Tentamen Numerieke Methoden voor Werktuigbouwkunde (2N460) op donderdag 24 juni 2010, 14.00-17.00 uur.
4. De midpuntregel voor de differentiaalvergelijking y 0 = f (t, y) luidt yk+1 = yk + hf tk + 21 h, 12 (yk + yk+1 ) ,
(4)
met yk de numerieke benadering van y(tk ) met tk = kh en h > 0 de stapgrootte. (a) Bepaal het stabiliteitsgebied van de midpuntregel aan de hand van het modelprobleem y 0 = λy met λ ∈ C. We willen de midpuntregel toepassen op het volgende beginwaardeprobleem x00 + ω 2 x + βx3 = cos t, x(0) = 1, x0 (0) = 0,
t > 0,
(5a) (5b)
met ω 1 en 0 < β 1, welke de gedwongen trilling van een niet-lineaire veer beschrijft. De variabele x = x(t) is de uitwijking van de veer uit evenwicht. Daartoe herschrijven we het beginwaardeprobleem (5) in de volgende vorm y0 = f (t, y), y(0) = c.
t > 0,
(6a) (6b)
(b) Geef de vectorfunctie f (t, y) en de vectoren y en c uit het beginwaardeprobleem (6). (c) Zijn de oplossingen van het stelsel differentiaalvergelijkingen (6a) stabiel? (d) Formuleer de midpuntregel voor het stelsel differentiaalvergelijkingen (6a). Laat zien dat de methode onvoorwaarlijk stabiel is. (e) Wat is een geschikte keuze voor de stapgrootte h? Motiveer uw antwoord.
Voor de vraagstukken kunnen de volgende aantallen punten worden behaald: 1.(a) 2 2.(a) 2 3.(a) 2 4.(a) (b) 3 (b) 3 (b) 2 (b) (c) 2 (c) 3 (c) 3 (c) (d) 3 (d) 2 (d) 3 (d) (e)
4
2 2 2 2 2
UITWERKINGEN 1.
. (a) In eerste orde benadering geldt yˆ − y = ∆xf 0 π2 . Differenti¨eren geeft f 0 (x) = sin x/f (x). Substitutie van x = π2 en |∆x| ≤ ε geeft de gevraagde ongelijkheid. (b) c(x) wordt gegeven door xf 0 (x) x sin x . c(x) = = f (x) 2(1 − cos x) Gebruik makend van de formules sin x = 2 sin 2 sin2 21 x kunnen we dit herschrijven als x cos 12 x , c(x) = 2 sin 12 x
1 x 2
cos
1 x 2
en 1 − cos x =
en het is duidelijk dat c(x) → 1 voor x → 0, m.a.w., het probleem y = f (x) is goed geconditioneerd voor x → 0. 1 2 (c) Uit een Taylorbenadering volgt f (x) = x 1 − 24 x + O x4 , m.a.w., f (x) = O(x) voor x → 0 en dit gedrag is af te lezen uit de grafiek voor 10−7 ≤ x ≤ 10−6 . Voor 10−8 ≤ x ≤ 10−7 treedt cijferverlies op in de term 1 − cos x, vandaar het ’staircase’ patroon. 1 2 (d) Een betere formule is f (x) = x 1 − 24 x (zie onderdeel (c)) of f (x) = 2 sin 21 x (zie onderdeel (b)) omdat hierbij geen cijferverlies optreedt. 2.
(a) Het differentieschema luidt −
1 u + u − 4u + u + u = f u . i,j−1 i−1,j i,j i+1,j i,j+1 i,j h2
(b) Het stelsel is nu lineair. Directe methoden zijn gebaseerd op de LU-decompositie A = LU. Bij het berekenen van deze LU-decompositie lopen de banden van L en U vol, zodat nzL nzA en nzU nzA, dus het oplossen van stelsels Ly = b en Ux = y is duur. In dit geval is A symmetrisch, positief definiet, zodat de geconjugeerde gradi¨entenmethode de beste methode is. De convergentie kan worden verbeterd door preconditionering toe te passen, bijvoorbeeld de onvolledige Cholesky decompositie. (c) We kiezen de volgende nummering van de (inwendige) roosterpunten (x2 , y2 ), . . . , (xN −1 , y2 ), . . . , (x2 , yN −1 ), . . . , (xN −1 , yN −1 ). De Gauss-Seidel methode wordt dan 1 (k+1) (k+1) (k+1) (k) (k) − 2 ui,j−1 + ui−1,j − 4ui,j + ui+1,j + ui,j+1 = 1 ⇔ h (k+1) (k+1) (k) (k) (k+1) ui,j = 41 ui,j−1 + ui−1,j + ui+1,j + ui,j+1 + h2 . We stoppen de iteratie wanneer het residue r := f − Au∗ met u∗ de numerieke oplossing voldoet aan ||r|| < tol. De toleratie tol kiezen we gelijk aan de geschatte discretisatiefout, dus O h2 . 5
(d) De Gauss-Seidel methode wordt nu −
1 (k+1) (k+1) (k+1) (k+1) (k) (k) . ui,j−1 + ui−1,j − 4ui,j + ui+1,j + ui,j+1 = cos ui,j 2 h (k+1)
De onbekende v = ui,j F (v) :=
4 v h2
voldoet aan de niet-lineaire vergelijking
− cos v − a = 0,
a :=
1 h2
(k+1) (k+1) (k) (k) ui,j−1 + ui−1,j + ui+1,j + ui,j+1 ,
welke we oplossen m.b.v. Newton iteratie, dus (k)
v (0) = ui,j voor
` = 0, 1, 2, . . .
v (`+1) = v (`) − F v (`) /F 0 v (`) , waarbij F 0 (v) = 4/h2 − sin v. Het stopcriterium wordt F (v (`) ≤ tol met dezelfde tolerantie als in onderdeel (c). 3.
(a) Ontwikkel f (x − h) en f (x − 2h) in Taylorreeksen. (b) We merken het volgende op 1. De afbreekfout 31 h2 f 000 (x) is dominant voor 10−6 ≤ h ≤ 1 2. In dit bereik is de helling 2, dus de afbreekfout is evenredig met h2 3. Afrondfouten t.g.v. cijferverlies, in combinatie met delen door 2h, zijn dominant voor 10−12 ≤ h ≤ 10−6 (c) Er geldt f 0 (1) = Dh + εh = Dh/2 +εh/2 met Dh = Dh [f ](1) en εh de bijbehorende fout etc. Omdat εh = Ch2 + O h3 geldt dat εh/2 ≈ 41 εh voor h voldoende klein. Combineren we dit met bovenstaande relatie, dan vinden als schatting voor de fout ε∗h = 43 Dh/2 − Dh en als betere benadering Dh∗ = Dh + ε∗h = 31 4Dh/2 − Dh . Deze benadering is derde orde nauwkeurig. (d) Een mogelijk MATLAB-script is clear all; format long e; k = linspace( 1,20,20 ); h = 2.ˆ(-k); fd = ( 3*exp(1)-4*exp(1-h)+exp(1-2*h) )./(2*h); err = abs( fd-exp(1) ); fd1 = ( 4*fd(2:20)-fd(1:19) )/3; err1 = abs( fd1-exp(1) ); loglog( h,err,’k’,h(1:19),err1,’k--’ ); grid; xlabel(’h’); ylabel(’error’); Het resultaat staat in de grafiek boven aan de volgende pagina. We merken het volgende op 1. Voor 10−3 ≤ h ≤ 1 is de afbreekfout dominant, de helling is 3, hetgeen een derde orde nauwkeurige benadering impliceert 6
0
10
−2
10
−4
error
10
−6
10
−8
10
−10
10
−12
10
−8
−6
10
10
−4
10 h
−2
0
10
10
2. Voor 10−6 ≤ h ≤ 10−3 zijn afrondfouten dominant, t.g.v. cijferverlies in combinatie met delen door 2h. 3. De minimale fout wordt bereikt voor grotere h. 4.
(a) Voor het modelprobleem geldt yk+1 = yk + hλ yk + yk+1 /2 ⇔ (1 − 12 hλ)yk+1 = (1 + 21 hλ)yk ⇔ 1 + 12 hλ yk . yk+1 = ψ(hλ)yk = 1 − 21 hλ Voor stabiliteit moet gelden ψ(hλ) ≤ 1, waaruit volgt dat |z + 2| ≤ |z − 2| voor z = hλ. Het stabiliteitsgebied is dus S = C− = {z ∈ C | Re(z) ≤ 0}. (b) Het beginwaardeprobleem wordt 0
y = f (y) =
! y2 , −ω 2 y1 − βy13 + cos t
1 y(0) = . 0
(c) De Jacobimatrix J(y) =
! 0 1 −ω 2 − 3βy12 0
p heeft de (zuiver imaginaire) eigenwaarden λ± = ±i ω 2 + 3βy12 , en derhalve zijn de oplossingen van (6a) stabiel. (d) De midpuntregel luidt y1,k+1 = y1,k + 12 h y2,k + y2,k+1 , y2,k+1 = y2,k + h − 21 ω 2 y1,k + y1,k+1 − β
1 2
y1,k + y1,k+1
3
+ cos tk + 12 h .
hλ± ∈ S voor alle h > 0, en dus is de midpuntregel onvoorwaarlijk stabiel.
7
(e) Stabiliteit speelt geen rol, immers de methode is onvoorwaarlijk stabiel. De keuze van h wordt dus bepaald door de gewenste nauwkeurigheid. We hebben te maken met twee zeer verschillende tijdschalen, t.w. τ1 = 1/ω, corresponderend met de eigentrilling van het systeem als β = 0, en τ2 = 1, corresponderend met de opgelegde kracht cos t. Het is evident dat τ1 τ2 . Een geschikte keuze voor h is bijvoorbeeld ωh = 5 × 10−2 .
8
