Mathematical Engineering - Matematica Numerica

Exercise 20

Divided by topic

MATEMATICA NUMERICA A.A. 2018 - 2019 Ingegneria Matematica Prof. A. Quarteroni Prof. A. Manzoni, Dr. I. Fumagalli Esercitazione 20 Metodi multipasso (2) e multistadio per eq. di. ordinarie Esercizio 1a)Ricavare il metodo AB2 di Adams-Bashforth ap+ 1 = 2passi per l'approssimazione numerica del generico problema di Cauchy y0 (t) =f(t; y(t)); t >0; y(0) =y 0: In particolare, siah >0il passo di tempo, eu nuna approssimazione di y(t n) al tempot n= nh. Detto p2 P pil polinomio che interpola i valori f nj= f(t nj; u nj) nei nodit nj, j= 0; : : : ; p, si imponga un+1= u n+Z tn+1 tn p( x)dx: Si verichi che lo schema così costruito risulta scritto nella forma seguente :un+1= u n+h2 [3 f n f n1] : b)Studiare l'ordine e la zero-stabilità del metodo AB2. Si verichi in Matlab, mediante la funzioneqssmultistep, l'ordine di convergenza considerando il problema y0 (t) =(y(t)1)2 ; t2(0;10); y(0) = 2;(1) la cui soluzione esatta èy(t) = 1 + 1=(t+ 1), e calcolando perh= [2 6 ;2 7 ;2 8 ]l'errore di appros- simazioneE h= max n=1;:::;N hj u n y(t n) j, ovet n= nh,N h= 10 =h. Per semplicità, si utilizzi la soluzione esatta per fornire i valori iniziali al metodo AB2. c)Si visualizzi la regione di assoluta stabilità del metodo AB2 utilizzando il comandoqssstabreg. Si determini quindi il valore del passo di discretizzazioneh 0tale che il metodo AB2 risulti assoluta- mente stabile per ogni0< h < h 0considerando il problema y0 (t) =10y(t); t >0; y(0) = 1: Riportare l'andamento della soluzione numerica assieme alla soluzione esatta nei due casih= 1:1h 0, h= 0:5h 0, commentando i risultati.1 Esercizio 2 Le equazioni che descrivono il moto di un pendolo di Foucault in assenza di attrito sono x00 2!sin( )y0 +k2 x= 0; y00 + 2!cos( )x0 +k2 y= 0;(2) dove è la latitudine del luogo in cui si trova il pendolo,!= 7:29
10 5 s-1 è la velocità angolare della Terra, k=pg=l cong= 9:8m/s2 . Si consideri la latitudine di Milano ( '=4) e si calcolino numericamente x=x(t)ey=y(t)pert2[0;300]s, con il metodo di Eulero esplicito. a)Si applichi la funzione Matlabeulero_avantiper la soluzione del problema in esame per diversi valori dih(per evitare calcoli troppo onerosi, si prenda sempreh1=1000). Si tracci il graco delle traiettorie nel piano delle fasixye il graco del modulo del vettore posizione(x; y). Cosa si può osservare ? Il metodo scelto è adatto per risolvere questo problema ? b)Si ripeta ora il punto precedente utilizzando il comando Matlabode45, che ha la rma seguente :[t, u] = ode45(f, [tstart, tend], u0) dove tstartetendsono gli istanti iniziale e nale,tè un vettorecolonnacontentente gli istanti di valutazione della soluzione numerica, e ognicolonnadiucontiene i valori di una componente della soluzione nel tempo. Si noti cheode45è un metodo adattivo, con scelta automatica del passoh. Esercizio 3 Si consideri un generico problema di Cauchy della forma (y0 (t) =f(t; y(t)); t2[t 0; T ]; y(t 0) = y 0; e l'approssimazionefu ngN ndella sua soluzione su una griglia di passo h= [Tt 0] =N, ottenuta mediante il seguente metodo di Runge-Kutta a 3 stadi : k1= f(t n; u n) k2= f tn+h2 ; u n+h2 k 1 k3= f tn+34 h; u n+34 hk 2 un+1= u n+29 hk 1+13 hk 2+49 hk 3 Per tale metodo, si tracci la regione di assoluta stabilità in Matlab, secondo i seguenti passi :a)Si consideri il problema modello per il qualef(t; y) =y, doveè una costante, e si trovi l'espressione della funzioneg:C!Rtale che un+1= g(h)u n: b)Si consideri la rappresentazione inR2 dei numeri complessizmediante le parti reale e immaginaria : z=a+ib. Alla funzionegrisulta associata la funzionee g:R2 !R,e g(a; b) =g(a+ib). Utilizzando i comandi Matlabmeshgridecontour, si tracci la regione di assoluta stabilità del metodo proposto, che corrisponde all'insieme di livello 1 dijg(a+ib)j.2 Esercizio 4 (Fuori programma) Una barra metallica rappresentata dall'intervallo monodimensionaleI= [a; a]scambia calore con l'ambiente circostante nell'intervallo di tempo[t 0; T ]. Consideriamo come noti : il coeciente di diusioneD==c p, che rappresenta la risposta termica del materiale ; dove la conduttività termica >0, la densitàe il calore specicoc psono costanti ; la sorgente di caloref(x; t)alla quale la barra è soggetta ; la temperaturag(t)alla quale gli estremi della barra sono mantenuti, per ognit2[t 0; T ]; la distribuzione di temperaturau 0( x)della barra al tempo inizialet 0, per ogni x2[a; a]. Il problema ai valori ai limiti e ai valori iniziali per l'equazione del calore nel caso in esame assume la seguente forma :8 > > < > > :@ u@ t D@ 2@ x 2u =f(x; t);in = (a; a)(t 0; T ]; u(a; t) =u(a; t) =g(t); t2(t 0; T ]; u(x; t 0) = u 0( x); x2(a; a): La formulazione debole di tale problema è la seguente : per ognit2(t 0; T ]trovareu(t)2V=H1 0( a; a) tale che(@ tu (t); v) +a(u(t); v) =F(v);8v2V ;(3) conu(0) =u 0, dove a(u; v) =DZ a au 0 (x)v0 (x)dx; F(v) =Z a af (x; t)vdxa(R g( x; t); v) essendoR g( x; t)2H1 (a; a)una funzione tale cheR g( x; t) =g(t)perx=a. Il secondo termine che compare nella denizione della formaF(v)è dovuto al rilevamento della condizione al contorno di Dirichlet non omogenea, e viene trattato automaticamente all'interno della libreria. a)Siaf(x; t) = 0e g(t) =1p 4 Dte a 24 Dt etu 0( x) =1p 4 Dt 0e x 24 Dt 0: In questo caso, la soluzione esatta risulta data da uex( x; t) =1p 4 Dte x 24 Dt : Usare una discretizzazione in spazio con elementi niti lineari(P 1) e il metodo di Eulero all'indietro (= 1) per la discretizzazione in tempo. Siano inoltrea= 3,t 0= 0 :5,T= 1eD= 0:1, e scegliere (h;
t)come segue : (h;
t) = (0:05;0:05);(0:025;0:025);(0:0125;0:0125) dovehè il passo di discretizzazione in spazio e
til passo di discretizzazione in tempo. Tracciare il graco dell'errore in funzione del passo di tempo (in scala logaritmica)ku ex u hk L1 (t 0;T ;H1 (a;a)) (memorizzato nella variabileerrorLinfH1restituita daParabolic1D_Solver.m) b)Ripetere lo stesso calcolo del punto precedente considerando elementi niti lineari (P 1) e = 0(Eulero in avanti) ; successivamente, scegliere elementi niti quadratici (P 2) e = 1=2(Crank-Nicolson). Cosa si osserva ? c)Considerare orat 0= 0 ,T= 1,D= 1, eu(3;0) = 0pert2(0;1], con u0( x) = 1x2[1;1] 0altrimenti; calcolare la soluzione numerica con il metodo di Eulero all'indietro(= 1)con
t= 0:1eh= 0:05. Cosa si osserva nel comportamento della soluzione rispetto ai dati iniziali ?3