TECHNISCHE UNIVERSITEIT EINDHOVEN Faculteit Wiskunde en Informatica Tentamen Numerieke Methoden voor Werktuigbouwkunde (2N460) op donderdag 23 juni 2011, 14.00-17.00 uur. Deel 1: Van 14.00 uur tot uiterlijk 15.30 uur. Het gebruik van het notebook is niet toegestaan. De uitwerkingen van de opgaven dienen duidelijk geformuleerd en beargumenteerd te worden.
1.
(a) Geef de definities van de norm ||A|| en het conditiegetal cond(A) van een matrix A. (b) Toon aan dat cond(A) = max ||A x||/ min ||A x||. ||x||=1
||x||=1
In het vervolg van deze opgave nemen we cos α − sin α A= , sin α cos α met
π 2
≤ α ≤ π. We noteren A = A(α). Beschouw het stelsel 1 A x = e1 met e1 = . 0
(c) Toon aan dat ||A(α)||2 = 1 en cond2 A(α) = 1. ˆ = A(ˆ ˆ de oplossing van Veronderstel dat α ˆ = α + ε met |ε| α en definieer A α). Zij x ˆ het verstoorde stelsel Aˆ x = e1 . (d) Toon aan dat in eerste orde benadering sin α cos α ˆ A = A − εB, B = . − cos α sin α (e) Toon aan dat in eerste orde benadering ||ˆ x − x||2 ≤ |ε|. ||x||2
1
Tentamen Numerieke Methoden voor Werktuigbouwkunde (2N460) op donderdag 23 juni 2011, 14.00-17.00 uur.
2. De variabele u = u(x) voldoet aan het volgende randwaardeprobleem − u00 = f (x), u(0) = 1,
0 < x < 1,
u0 (1) = 0,
met f (x) = 12 1 + tanh 10(x − 21 ) . We willen een numerieke benadering voor u(x) berekenen op het rooster xi = (i − 1)h,
i = 1, 2, · · · , N,
met h = 1/(N − 1) de maaswijdte. Voor de discretisatie van u0 en u00 maken we gebruik van centrale differenties. De numerieke benadering van u(xi ) geven we aan met ui . (a) Toon aan m.b.v. Taylorreeksen dat u00 (xi ) =
1 u(x ) − 2u(x ) + u(x ) − i+1 i i−1 h2
1 2 (4) h u (xi ) 12
+ O(h4 ).
(b) Geef het differentieschema voor het randwaardeprobleem in de roosterpunten xi (i = 2, 3, . . . , N ). Zij u = (u2 , u3 , . . . , uN )T de numerieke oplossing. Hiervoor kunnen we het lineaire stelsel Au = b afleiden. (c) Bepaal de matrix A en de vector b. (d) Formuleer een effici¨ent algoritme voor het oplossen van het stelsel Au = b.
2
Tentamen Numerieke Methoden voor Werktuigbouwkunde (2N460) op donderdag 23 juni 2011, 14.00-17.00 uur. Deel 2: Direct aansluitend op deel 1. Het gebruik van het notebook is wel toegestaan.
3. Een bol met straal R en soortelijke massa ρb drijft op het water. Het volume van de bol is 4 πR3 . Gegeven is dat ρb = ερw met 0 < ε < 0.5 en ρw de soortelijke massa van water. 3 1 2 Het volume van het verplaatste water is gelijk aan 3 πh 3R − h met h de hoogte van het gedeelte van de bol onder water, gemeten langs de verticaal door het middelpunt. We defini¨eren z = h/R. (a) Toon aan dat z voldoet aan de vergelijking f (z) = z 3 − 3z 2 + 4ε = 0. (b) Toon aan dat f precies een nulpunt z ∗ ∈ (0, 1) heeft. (c) Formuleer de methode van Newton voor het berekenen van z ∗ . Geef een beginschatting en een stopcriterium. (d) Implementeer de methode van Newton uit onderdeel (c) in een M ATLAB-script. Schrijf uw script volledig uit. Bereken z ∗ in vier significante cijfers voor het geval ε = 0.3. 4. Voor de berekening van de integraal Z b I(f ) = f (x) dx, (a < b),
(1)
a
maken we gebruik van de (samengestelde) midpuntregel. De enkelvoudige midpuntregel luidt 1 I(f ) = (b − a)f (m) + (b − a)3 f 00 (ξ), ξ ∈ (a, b), 24 met m = 21 (a + b). De samengestelde midpuntregel voor (1) kunnen we afleiden door het interval [a, b] in k gelijke deelintervallen van lengte h te splitsen en op elk van deze intervallen de midpuntregel toe te passen. Uiteraard geldt nu dat h = (b − a)/k. De samengestelde midpuntregel Mh (f ) met bijbehorende foutterm εh (f ) wordt gegeven door: I(f ) = Mh (f ) + εh (f ), Mh (f ) = h
k X
(2a)
f (mi ),
(2b)
i=1
εh (f ) =
1 (b − a)h2 f 00 (ξ), 24
ξ ∈ (a, b),
(2c)
waarbij mi = (xi +xi−1 )/2 (i = 1, 2, . . . , k) en xi = a+ih (i = 0, 1, 2, · · · , k). Formule (2c) geldt mits f twee keer continu differentieerbaar is op (a, b), hetgeen we aannemen voor deze opgave.
3
Tentamen Numerieke Methoden voor Werktuigbouwkunde (2N460) op donderdag 23 juni 2011, 14.00-17.00 uur.
4.
(a) Geef een afleiding van de samengestelde midpuntregel (2). (b) Toon aan m.b.v. Richardson extrapolatie dat voor h voldoende klein qh =
Mh/4 (f ) − Mh/2 (f ) 1 ≈ . Mh/2 (f ) − Mh (f ) 4
Voor de berekening van Z 1 cosh x dx
(3)
0
hebben we de beschikking over het volgende, onvolledige M ATLAB-script. clear all; format long e; f = @(x) cosh(x); a = 0; b = 1; n = 6; H = [ ]; Mf = [ ]; for j = 1:n k = 5*(2ˆj); h = (b-a)/k; H = [ H;h ]; . . . end qh = ( Mf(3:n)-Mf(2:n-1) )./( Mf(2:n-1)-Mf(1:n-2) ); (c) Completeer het script voor de berekening van de integraal in (3). Schrijf uw script volledig uit. Voer het resulterende script uit en geef de resultaten voor Mh (f ) en qh in vier significante cijfers. Verklaar de resultaten. (d) Geef een schatting voor de fout εh (f ) voor h = 1/320.
Voor de opgaven kunnen de volgende aantallen punten worden behaald: 1.(a) (b) (c) (d) (e)
2 2 2 2 2
2.(a) 3 (b) 3 (c) 2 (d) 2
3.(a) 2 (b) 2 (c) 3 (d) 3
4.(a) 2 (b) 3 (c) 3 (d) 2
Het eindcijfer wordt bepaald door 0.3p + 0.7t af te ronden, met p het practicumcijfer en t het tentamencijfer. 4
UITWERKINGEN 1.
(a) De definities van norm en conditiegetal zijn: ||A|| = max x6=0
||Ax|| ||x||
en
cond(A) = ||A|| ||A−1 ||.
(b) Zie lecture3.pdf. (c) Stel y = Ax, dan is x1 cos α − x2 sin α y= . x1 sin α + x2 cos α De 2-norm van y wordt gegeven door ||y||22 = y12 +y22 = (x1 cos α−x2 sin α)2 +(x1 sin α+x2 cos α)2 = x21 +x22 = ||x||22 , zodat ||Ax||2 = ||x||2 voor elke vector x. Uit de definitie in onderdeel (a) volgt dan dat ||A||2 = 1. Uit het voorgaande concluderen we dat ook ||A−1 ||2 = 1, zodat cond2 (A) = 1. (d) Dit volgt rechtstreeks uit de eerste orde Taylorbenaderingen cos(α + ε) = cos α − ε sin α + O ε2 , sin(α + ε) = sin α + ε cos α + O ε2 . (e) De relatieve fout in x t.g.v. een verstoring in A wordt gegeven door |ε|||B||2 ||ˆ x − x||2 ≤ cond2 (A) = |ε|, ||x||2 ||A||2 omdat ||A||2 = 1, cond(A)2 = 1 (zie onderdeel (c)) en ook ||B||2 = 1. Het bewijs van de laatste gelijkheid verloopt volkomen analoog aan het bewijs van ||A||2 = 1. 2.
(a) Ontwikkel u(xi±1 ) in Taylorreeksen. (b) Het differentieschema luidt 1 u − 2u + u = f (xi ), i+1 i i−1 h2 1 − 2 − 2uN + 2uN −1 = f (xN ), h −
i = 2, 3, . . . , N − 1,
met als randvoorwaarde u1 = 0. (c) De matrix A en de vector b worden gegeven door −2 1 f (x2 ) 1 −2 1 f (x3 ) 1 1 .. ... ... ... A=− 2 , b = + 2 e1 . . h h f (xN −1 ) 1 −2 1 2 −2 f (xN ) 5
(d) De matrix A is een tridiagonaalmatrix, waarvoor LU-decompositie de meest effici¨ente oplosmethode is (geen fill in). 3.
(a) Volgens de wet van Archimedes geldt Fz − Fopw = 0 met Fz de zwaartekracht en Fopw de opwaartse kracht werkend op de bol, welke worden gegeven door Fz = ρb 34 πR3 g,
Fopw = ρw 31 πh2 (3R − h)g,
met g de gravitatieversnelling. Substitutie van bovenstaande uitdrukkingen in de krachtenbalans geeft 4εR3 − h2 (3R − h) = 0 en delen door R3 levert tenslotte de gevraagde vergelijking. (b) De functie f is continu op (0, 1) en voldoet aan f (0) = 4ε > 0,
f (1) = 4ε − 2 < 0,
f 0 (z) = 3z(z − 2) < 0.
De functie f is monotoon dalend op (0, 1). Uit het voorgaande concluderen we dat f precies een nulpunt z ∗ ∈ (0, 1) heeft. (c) De methode van Newton luidt zk+1
zk2 zk − 3 + 4ε f (zk ) . = zk − 0 = zk − f (zk ) 3zk zk − 2
Een geschikte beginschatting is z = 0.5, corresponderend met ε = 0.5 en een mogelijk stopcriterium is |zk+1 − zk | < tol en/of |f (zk )| < tol voor een zekere tolerantie tol, bijvoorbeeld tol = 10−8 . (d) Een mogelijk script is: clear all; format long e; fid = fopen( ’sc3juni11.out’, ’w’ ); fprintf( fid, ’ k z(k) \n\n’ ); epsilon = 0.3; maxit = 50; tol = 1e-8; k = 0; z = input(’beginschatting: ’); conv = 0; while (˜conv) zold = z; k = k+1; fz = zˆ3 - 3*zˆ2 + 4*epsilon; dfz = 3*zˆ2 - 6*z; dz = fz/dfz; z = z - dz; fprintf( fid, ’%4i %20.12e\n’, k, z ); conv = ( abs((z-zold)/zold)
met als uitvoer k 1 2 3 4 5
z(k) 7.555555555556e-01 7.267429193900e-01 7.265149975282e-01 7.265149821811e-01 7.265149821811e-01
We vinden nu z ∗ = 7.265 × 10−1 . 4.
(a) De afleiding van de samengestelde midpuntregel gaat als volgt I(f ) =
=
k Z X i=1 k X
f (x) dx
xi−1 1 3 00 h f (ξi ) 24
hf (mi ) +
i=1 k X
=h
=h
xi
i=1 k X
f (mi ) +
1 3 h 24
k X
ξi ∈ (xi−1 , xi )
f 00 (ξi )
i=1
f (mi ) +
1 (b 24
− a)h2 f 00 (ξ)
ξ ∈ (a, b).
i=1
(b) Toepassen van (2) voor h en h/2 levert I(f ) = Mh (f ) + εh (f ) = Mh/2 (f ) + εh/2 (f ). Voor h voldoende klein geldt dat εh/2 (f ) ≈ 41 εh (f ) omdat de samengestelde midpuntregel tweede orde nauwkeurig is. Hieruit volgt 3 ε (f ) 4 h
≈ Mh/2 (f ) − Mh (f ).
Op analoge manier vinden we 3 ε (f ) 4 h/2
≈ Mh/4 (f ) − Mh/2 (f ).
Delen van beide relaties geeft dan de gevraagde benadering.
7
(c) Het volledige script luidt clear all; format long e; fid = fopen( ’integraal.out’,’w’ ); f = @(x) cosh(x); a = 0; b = 1; n = 6; H = [ ]; Mf = [ ]; for j = 1:n k = 5*(2ˆj); h = (b-a)/k; H = [ H;h ]; x = a:h:b; m = 0.5*( x(1:k) + x(2:k+1) ); fm = f(m); mf = h*sum( fm ); Mf = [ Mf;mf ]; fprintf( ’%4i %15.4e %15.4e \n’,k,h,mf ); end fprintf( ’\n’ ); qh = ( Mf(3:n)-Mf(2:n-1) )./( Mf(2:n-1)-Mf(1:n-2) ); for k = 1:n-2 fprintf( ’%15.4e \n’,qh(k) ); end fclose( fid ); met als uitvoer: 10 20 40 80 160 320
1.0000e-01 5.0000e-02 2.5000e-02 1.2500e-02 6.2500e-03 3.1250e-03 2.5007e-01 2.5002e-01 2.5000e-01 2.5000e-01
1.1747e+00 1.1751e+00 1.1752e+00 1.1752e+00 1.1752e+00 1.1752e+00
De tweede kolom bevat de waarden van h, de derde kolom de bijbehorende waarden Mh (f ) en de vierde kolom de quoti¨enten qh , gedefinieerd in onderdeel (b). Hieruit concluderen we dat de methode inderdaad tweede orde nauwkeurig is. (d) Toepassen van Richardson extrapolatie geeft εh (f ) ≈ 13 Mh (f ) − M2h (f ) ; zie ook onderdeel (b). Passen we dit toe voor h = 1/320 dan vinden we in vier significante cijfers de foutschatting 4.782 × 10−7 .
8
