Aerospace Engineering - Calcolo Numerico ed Elementi di Analisi

Second partial exam

Seconda Prova in Itinere 21/06/2022 — versione 1 — }~|}|~| 32 pt – durata 1h 30’ – MS Forms Gli studenti aventi diritto a svolgere la prova ridotta del 30% secondo la L.170/2010 (indicazioni Multichance team) NON svolgono i quesiti contrassegnati con (***) TEST – 15 pt 1 — 2 pt Siano date le 3 coppie di dati {(0,2),(0.5,1),(1,1.5)}. Si costruisca il polinomio interpolante⇧ 2(x) di tali dati e si riportino i valori di⇧ 2(0.25) e d⇧2 dx (0.25). 1.3125, 2 2 — 1 pt Siano date le n+1 = 5 coppie di dati {(0,3),(0.25,0.5),(0.5,1.5),(0.75,0.5),(1,1)}. Qual `e valore minimo dell’interpolante polinomiale lineare a tratti⇧ H1(x)deidati precedenti per x2[0,1]? 0.5 3 — 1 pt Si consideri la funzione f(x)=sin( ⇡x ) e il suo interpolante polinomiale lineare a tratti⇧ H1f(x) su 4 sottointervalli equispaziati di [0 ,2]. Si riporti il valore⇧ H1f(1.6). 0.8 4 — 2 pt (***) No Multichance Si considerino le coppie di dati nella forma {(xj,y j)}nj=0,talipercui n+1 = 5 especificamente {(0,3),(0.25,0.5),(0.5,1.5),(0.75,0.5),(1,1)}. Si costruisca il polinomio p2(x) di grado 2 approssimante tali dati nel senso dei minimi quadrati e si riporti il valore dello scarto quadratico nX j=0 (p2(xj)yj)2. 2.4143 1 5 — 1 pt Si consideri l’approssimazione dell’integrale I(f)= Z2 1f(x)dx,con f(x)= p2+ |x|, tramite la formula dei trapezi composita per sottointervalli equispaziati di ampiezza H = 1. Si riporti il valore approssimato dalla formula IHt(f)cos`ıottenuto. 5.0123 6 — 1 pt (***) No Multichance Si consideri l’approssimazione dell’integrale I(f)= Zb af(x)dx,dove f2C1([a, b ]), tramite una formula di quadratura composita accurata di ordine p= 2 rispetto all’ampiezza H dei sottointervalli. Sapendo che per M1= 20 sottointervalli equis- paziati di [ a, b ] si ha un errore pari a e1(f)=10 1, per quanti sottointervalli M2 si ha un errore stimato e2(f)=10 3? 200 7 — 2 pt (***) No Multichance Si consideri il seguente problema di Cauchy: ( y0(t)= 1 81 (1 + sin( t)) ( y(t))2 t2(0,5), y(0) = 9 . Utilizzando il metodo di Eulero all’indietro basato sul metodo di Newton imple- mentato nella funzione eulero indietro newton.m epasso h=0 .2, si riporti il valore calcolato di uNt, ovvero l’approssimazione di y(5), essendo tn= nh per n=0 ,...,N t. 5.5989 8 — 2 pt Si consideri il seguente problema ai limiti: ⇢ u00(x)+2 u(x)=3+ xx 2(0,1), u(0) = 1 ,u 0(1) = 0 . Si approssimi il problema utilizzando il metodo delle di↵erenze finite centrate con passo di discretizzazione h=1 /10 e il metodo delle di↵erenze finite all’indietro (di ordine 1) per l’approssimazione di u0(1), ottenendo la soluzione numerica {uj}N+1j=0 nei corrispondenti nodi {xj}N+1j=0 per ( N + 1) = 10. Si risolva il problema e si riporti il valore della soluzione numerica u10, ovvero l’approssimazione di u(1). 1.4207 2 9 — 2 pt Si consideri il seguente problema ai limiti: ⇢ u00(x)50 u0(x)=0 x2(0,1), u(0) = 0 ,u (1) = 3 . Si approssimi il problema utilizzando il metodo delle di↵erenze finite centrate con tecnica Upwind e passo di discretizzazione h=1 /10, ottenendo la soluzione nu- merica {uj}N+1j=0 nei corrispondenti nodi {xj}N+1j=0 per ( N + 1) = 10. Si risolva il problema e si riporti il valore della soluzione numerica u1, ovvero l’approssimazione di u(x1). 2.5000 10 — 1 pt Si consideri il seguente problema di di↵usione: 8>>>< >>>: @u @t (x, t )3@2u @x 2(x, t )=0 x2(0,1),t> 0, u(0,t)= u(1,t)=0 t> 0, u(x,0) = 6 sin( ⇡x ) x2(0,1). Si consideri l’approssimazione di tale problema tramite il metodo delle di↵erenze finite centrate con passo di discretizzazione spaziale h=0 .2 e il metodo di Eulero in avanti con passo di discretizzazione temporale t> 0. Per quali valori di t`e garantita l’assoluta stabilit`a? t< 0.0067 3 ESERCIZIO – 17 pt Si consideri il sistema di Equazioni Di↵erenziali Ordinarie del primo ordine: 8< : dy dt (t)= A(y(t)) y(t)+ g(t) t2(0,tf), y(0) = y0, (1) dove y(t)=( y1(t),y 2(t))T,y0= ✓ 3 4,3 ◆T ,A(y)= 2 4 1 ✓147 16 +Ky 2 ◆ 10 3 5 per un parametro K 2R,g= ✓ 9et/2sin 2(3t)1 9 2et/4sin(3 t),0 ◆T e tf=5. Punto 1) — 2 pt (***) No Multichance Con riferimento a un generico sistema di Equazioni Di↵erenziali Ordinarie nella forma (1), si riporti la definizione di zero-stabilit`a in relazione al metodo di Eulero in avanti. Si definisca tutta la notazione utilizzata. Spazio per risposta lunga Punto 2) — 3 pt Dopo aver posto K = 1, si approssimi il problema (1) tramite il metodo di Eulero in avanti usando opportunamente la funzione Matlab reulero avanti sistemi.m con passo h=10 2. Si riportino, oltre ai comandi Matlab r: •i valori delle approssimazioni u1euNhrispettivamente di y(t1)e y(tf), dove tn=nh per n=0 ,1,...,N h,h= tf Nh; •il valore minimo u2,min =min n=0,...,N h(un)2e il tempo discreto tmin acuisi ottiene tale valore u2,min . u1=( 1.0181 ,2.9925) T, uNh=( 1.6865 ,0.7539) T, u2,min = 2.4321, tmin =1 .0300 Spazio per risposta lunga Punto 3) — 2 pt Dopo aver risposto al Punto 2) e sapendo che la soluzione esatta del problema (1) `e y(t)=3 et/4✓ 1 4cos(3 t)3sin(3 t),cos(3 t) ◆T si calcolino gli errori Eh=kuNhy(tf)kottenuti con il metodo di Eulero in avanti e corrispondenti ai passi h1=10 3,h2=5 ·104,h3=2 .5·104eh4=1 .25 ·104. Si riportino i valori Ehiper i=1 ,..., 4eicomandiMatlab rusati. 0.0182 , 0.0090 , 0.0045 , 0.0023 Spazio per risposta lunga 4 Punto 4) — 2 pt Si utilizzino gli errori Ehiottenuti al Punto 3) per stimare algebricamente l’ordine di convergenza pdel metodo di Eulero in avanti. Si giustifichi la risposta data e la si motivi alla luce della teoria. Si riportino i comandi Matlab rconsiderati. p=1 .0018 Spazio per risposta lunga Punto 5) — 3 pt Si vuole ora applicare a un sistema di Equazioni Di↵erenziali Ordinarie nella forma (1) il metodo di Runge-Kutta associato alla seguente tabella di Butcher 0 00 2/3 2/30 1/43 /4 Si scriva un’opportuna funzione Matlab rche implementi il metodo precedente per un sistema di EDO con il generico passo h> 0. Si utilizzi la funzione precedentemente implementata per risolvere il sistema di Equazioni Di↵erenziali Ordinarie (1) con K =1eusandoilpasso h=10 2.Si riportino i valori delle approssimazioni u1euNhrispettivamente di y(t1)e y(tf). Si riportino i comandi Matlab rconsiderati. u1=( 1.0171 ,2.9912) T,uNh=( 1.5135 ,0.6530) T Spazio per risposta lunga Punto 6) — 2 pt Con riferimento al sistema di Equazioni Di↵erenziali Ordinarie (1), ma col parametro K =0 e g=0, per quali valori di h> 0ilmetododi Heun risulta assolutamente stabile? Si giustifichi dettagliatamente la risposta data riportando gli eventuali comandi Matlab rusati. 0>>< >>>: Un+1i Uni t + µ h2(2 Uni Uni1 Uni+1 )=0 ,i =1 ,..., 4 Un0 = Un5 =0 ,n 0 U0i=6sin( ⇡x i),i =1 ,..., 4. In forma vettoriale, questo corrisponde a 8 >< >: U n+1 U n t + µ h2A U n= 0 n 0 U 0= U 0=6[sin( ⇡x 1),sin( ⇡x 2),sin( ⇡x 3),sin( ⇡x 4)]T, dove si `e indicato con U n=[ Un1,U n2,U n3,U n4]T il vettore delle incognite al tempo tn econ A la matrice tridiag( 1,2,1). Risolvendo per U n+1 ,siottiene U n+1 = ⇣ I µt h2 A ⌘ U n,n 0,con U 0= U 0, 5 in cui I indica la matrice identit`a di ordine 4. Questa `e una relazione ricorsiva che ha per soluzione U n= ⇣ I µt h2 A ⌘n U 0. Anch´e il metodo risulti stabile, `e necessario che la matrice B = I µth2 A abbia raggio spettrale minore di uno, i.e., ⇢(B)< 1. Dato che A `e s d p , l a c o n d i z i o n e diventa ⇢(B)< 1 () 1 µt h2 < 1 () 1< 1 µt h2 < 1, dove si `e indicato con un generico autovalore (positivo) di A,i.e, 2 sp( A). Risulta quindi 0< t< 2h2 µ , 8 2 sp( A). (2) La condizione pi`u restrittiva, corrispondente all’autovalore maggiore, pu`o essere determinata con lo script Matlab h =0.2; mu =3; n =4; A =2 ⇤eye (n)diag (ones (n 1,1), 1) diag (ones (n 1,1) ,1) ; e = eig (A); 2⇤hˆ2/( mu ⇤max (e)) efornisce t< 0.0074. Tale condizione `e meno restrittiva della stima (1), che di fatto corrisponde al caso limite della (2) per n !1 con il vincolo hn =1. Esercizio - Punto 1 Il problema di↵erenziale del testo, che qui riportiamo per averlo sotto mano: 8 < : dy dt (t)= A(y(t))y(t)+ g(t),t 2 (0 ,tf) y(0) = y0, (3) con A = A(y(t)) 2 R2⇥2,g(t):(0 ,tf)! R2,y02 R2il vettore iniziale e tf> 0il tempo finale. Introduciamo la quantit`a un 2 R2 per indicare l’approssimazione del vettore y(tn), e indichiamo con Nhil numero di intervalli in cui viene suddiviso l’intervallo [0,tf], cui corrisponde il passo temporale, h = tf/N h. Il metodo di EA applicato al problema (3) risulta pertanto: trova un+1 tale che ⇢ un+1 = un+ h(A(un)un+ g(tn)), per n =0 ,...,N h 1, u0 = y0. (4) 6 Per definire la zero-stabilit`a, o ccorre intro durre il meto do di EA p erturbato: trova zn+1 tale che ⇢ zn+1 = zn+ h(A(zn)zn+ g(tn)) + n+1 , per n =0 ,...,N h 1, z0 = y0+ 0, (5) dipendente dalle perturbazioni (vettoriali), n,per n =0 ,...,N h.Ladefinizionedi zero-stabilit`a, che ricalca quella nel caso scalare per una f(t, y )generica,confronta la soluzione di (4) con quella di (5) e sar`a pertanto la seguente: Esistono h0> 0e "0> 0taliche maxn=0 ,...,N h||zn un|| C", sotto la condizione che ||n|| ",per n =0 ,...,N h,perogni "2 (0 ," 0], per ogni h 2 (0 ,h 0], con la costante C indipendente da h,edove || · || indica una opportuna norma vettoriale. Esercizio - Punto 2 Lo script seguente h =1 e2; K =1; A = @(y)[ 1 (147/16 + K⇤y(2) ) ; 1 0 ] ; g = @(t)[ 9⇤exp ( t/2) ⇤(sin (3 ⇤t)ˆ2 1) . . . 9/2 ⇤exp ( t/4) ⇤ sin (3 ⇤t);0]; f = @(t,y) A(y)⇤y + g(t); y0 =[ 3/4;3]; tf =5; Nh = tf /h; tv =[0, tf ]; [ t, u ]= eulero avanti sistemi ( f, tv , y0 , Nh ); u(: ,[2 , end ]) [m,mm ]= min (u(2 ,:) ) ; m, t(mm ) implementa il metodo di Eulero in avanti e permette di calcolare le quantit`a: u1= (1.0181 ,2.9925) T euNh=( 1.6865 ,0.7539) T,u2,min = 2.4321 e tmin =1 .0300. Esercizio - Punto 3 Il calcolo degli errori pu`o essere e↵ettuato attraverso il seguente script 7 yex = @(t)3 ⇤exp ( t/4) ⇤[ 1/4 ⇤cos (3 ⇤t)3⇤sin (3 ⇤t); cos (3 ⇤t)]; K =1; A = @(y)[ 1 (147/16 + K⇤y(2) ) ; 1 0 ] ; g = @(t)[ 9⇤exp ( t/2) ⇤(sin (3 ⇤t)ˆ2 1) . . . 9/2 ⇤exp ( t/4) ⇤ sin (3 ⇤t);0]; f = @(t,y) A(y)⇤y + g(t); y0 =[ 3/4;3]; tf =5; tv =[0, tf ]; H =[1 e3, 5 e4, 2.5 e4, 1.25 e 4]; E = zeros (1 ,4) ; for i =1:4 h = H(i); Nh = tf /h; [ t, u ]= eulero avanti sistemi ( f, tv , y0 , Nh ); E(i)= norm (u(: , end ) yex (tf )) ; end format long E Ivaloriottenutisono: 0 .0181789, 0 .0090431, 0 .0045101 e 0 .0022522. Esercizio - Punto 4 Per stimare algebricamente l’ordine di convergenza del metodo numerico, basta os- servare che l’errore pu`o essere scritto come Eh' Ch p,perunacertacostantepositiva C econ pl’ordine di convergenza. Dal confronto fra gli errori associati a due passi diversi `e possibile eliminare C estimarel’ordine. Infatti Eh1 Eh2 ' ⇣h1 h2 ⌘p =) p' log Eh1 Eh2 log h1h2 . Applicando questa procedura al caso in esame, confrontando gli errori a due a due, otteniamo che Eh1 Eh2 =2 .0102 = ) p' log 2 .0102log 2 =1 .0074 , Eh2 Eh3 =2 .0051 = ) p' log 2 .0051log 2 =1 .0037 , Eh3 Eh4 =2 .0025 = ) p' log 2 .0025log 2 =1 .0018 . 8 Consideriamo come riferimento l’ultimo valore, perch´e associato ai valori pi`u piccoli del passo. Possiamo quindi assumere che p' 1.0018. Questo valore `e coerente con quanto ci si aspetta dalla teoria del metodo di EA, che prevede infatti ordine 1. Esercizio - Punto 5 Il metodo RK proposto `e un metodo esplicito a due stadi con coecienti: b1=1 /4, b2=3 /4, c1=0, c2=2 /3, a21 =2 /3, e tutti gli altri aij nulli. In forma estesa, si scrive: assegnato u0= y0,per n =0 ,...,N h 1, calcola 8 >< >: K 1= f(tn,un) K 2= f(tn+2 h/ 3,un+(2 h/ 3) K 1) un+1 = un+ h 4(K 1+3 K 2). Si osservi che il metodo risulta correttamente definito anche per quantit`a vettoriali. Lo script seguente clearvars K =1; A = @(y)[ 1 (147/16 + K⇤y(2) ) ; 1 0 ] ; g = @(t)[ 9⇤exp ( t/2) ⇤(sin (3 ⇤t)ˆ2 1) . . . 9/2 ⇤exp ( t/4) ⇤ sin (3 ⇤t);0]; f = @(t,y) A(y)⇤y + g(t); y0 =[ 3/4;3]; tf =5; h =1 e2; Nh = tf /h; th = linspace (0 , tf ,Nh +1) ; K1 = @(t,y) f(t,y); K2 = @(t,y) f(t +2/3 ⇤h, y +2/3 ⇤h⇤K1 (t,y)) ; u = zeros (2 , Nh +1) ; u(: ,1) = y0 ; for i =1: Nh u(: , i+1) = u(: , i)+ h⇤(K1 (th (i), u(: , i)) ... +3 ⇤K2 (th (i), u(: , i)))/4; end u(: ,[2 , end ]) implementa il metodo proposto e fornisce u1=( 1.0171 ,2.9912) T e uNh=( 1.5104 ,0.6534) T. 9 Esercizio - Punto 6 Con K =0e g = 0,lafunzionematriciale A si riduce ad una matrice algebrica costante, A =[ 1,147 /16; 1 ,0] e il problema di↵erenziale (3) diventa 8 < : dy dt (t)= A y(t),t 2 (0 ,tf) y(0) = y0, (6) che rappresenta un sistema di equazioni di↵erenziali di ordine uno a co ecienti costanti. Il metodo di Heun applicato a tale problema si scrive 8 >>< >>: u⇤n+1 = un+ hA un un+1 = un+ h 2 Aun+ Au⇤n+1 ,n 0 u0 = y0. Sostituendo u⇤n+1 dalla prima nella seconda relazione, si ottiene un+1 = un+ h 2 ⇣ Aun+ Aun+ hA un ⌘ = ⇣ I+ hA + (hA )2 2 ⌘ un, dove I `e l a m a t r i c e i d e n t i t `a d i o r d i n e d u e . Q u e s t a r e l a z i o n e r i c o r s i v a h a c o m e soluzione un= ⇣ I+ hA + (hA )2 2 ⌘n y0,n 0. La condizione di assoluta stabilit`a per il metodo di Heun equivale dunque a richiedere che ⇢(B)< 1con B = I+ hA + (hA )2 2 ,essendo ⇢(·)ilraggiospettrale. Sihaquindi ⇢(B)< 1 () 1+ h + (h )2 2 < 1 8 2 sp( A). Dato che gli autovalori della matrice A sono complessi coniugati, conviene passare aMatlab. Loscript A =[ 1 147/16; 1 0 ] ; e = eig (A); h = linspace (0 ,1 ,1000) ; R = abs (1 + h⇤e(1) + ( h⇤e(1) ) .ˆ2/2) ; plot (h,R) grid costruisce la Figura 1 in cui `e rappresentata la quantit`a R(h )= 1+ h + (h)2 2 al va r i a r e d i h,perilprimoautovaloredi A (il secondo autovalore produrrebbe lo stesso grafico). Zoomando nell’intorno dell’ascissa 0 .4sipu`ostimarechel’attraversamento della quota unitaria si ha per h ' 0.4247. In definitiva, la condizione di assoluta stabilit`a per il metodo di Heun `e 0