Mathematical Engineering - Matematica Numerica

Exercise - 20 - Solution

Divided by topic

MATEMATICA NUMERICA A.A. 2018 - 2019 Ingegneria Matematica Prof. A. Quarteroni Prof. A. Manzoni, Dr. I. Fumagalli Esercitazione 20 - Soluzione 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. Soluzione 1a)Nel nostro casop= 1; il polinomio di grado 1 che interpola i dati(t n1; f n1) e(t n; f n) si scrive 1( x) =f n1t n xh + f nx t n1h :1 Essendo di grado uno, possiamo applicare la formula del trapezio per integrarlo esattamente su [t n; t n+1] . Si trova allora Ztn+1 tn 1( x)dx=h2 f n+h2 ( f n1+ 2 f n) = h 32 f n12 f n1 ; e quindi, imponendo che risultiu n+1= u n+R tn+1 tn 1( x)dx, otteniamo il metodo AB2 un+1= u n+h2 [3 f n f n1] : b)Il metodo AB2 appartiene alla famiglia generale dei metodi ap+ 1passi lineari della forma : un+1=p X j=0a ju nj+ hp X j=1b jf nj; n =p; p+ 1; : : : dovep= 1e i coecientia je b jsono dati da : a0= 1 ; a 1= 0 ; b1= 0 ; b 0=32 ; b 1= 12 Cominciando ad analizzare la consistenza, otteniamop X j=0a j= 1 + 0 = 1 ; e inoltrep X j=0j a j+p X j=1b j=3 12 = 1 :(2) Il metodo AB2 è quindi consistente. Studiamone ora l'ordine, vericando in Matlab le condizioni aggiuntive sui coecienti, come segue :a = [1; 0; ]; b = [0; 3/2; -1/2]; p=1; for i=1:3%finoai=3val = ((-[0:p]).^i) *a + i *( ((-[-1:p]).^(i-1)) *b );display(sprintf(i=%i,valore=%f, i, val)); end i = 1, valore = 1.000000 i = 2, valore = 1.000000 i = 3, valore = -1.500000 Si vede dunque che le relazioni algebriche di interesse sono rispettate per i= 1(in tal caso esse coincidono con la relazione di consistenza (2) dato chej0 = 1,j=1; : : : ; p), sono rispettate anche peri= 2, mentre non lo sono sei= 3. Dunque il metodo ha ordine 2 (come noto per i metodi AB a p+ 1passi, l'ordine èp+ 1). Analizziamone ora la zero-stabilità. A tale scopo introduciamo il primo polinomio caratteristico, che nel caso in esame si scrive : (r) =r2 r=r(r1) e osserviamo cher 0= 1 con molteplicità 1 (la radice di consistenza) mentrer 1= 0 e dunquejr 1j > < > > > :x 0 (t) =z(t); y0 (t) =v(t); z0 (t) = 2!sin( )v(t)k2 x(t); v0 (t) =2!cos( )z(t)k2 y(t):(4) Se supponiamo che il pendolo all'istante inizialet 0= 0 sia fermo nella posizione(1;0), il sistema (4) viene completato dalle seguenti condizioni iniziali : x(0) = 1; y(0) = 0; z(0) = 0; v(0) = 0: Deniamo dunque in Matlab il vettore soluzioneu'[x; y; z; v]T e applichiamo il metodo di Eulero esplicito al problema in esame, con le seguenti istruzioni (si prenda ad esempiol= 20m)l = 20; k2 = 9.8/l; psi = pi/4; omega = 7.29e-5; f = @(x) [ u(3); u(4); 2 *omega *sin(psi) *u(4)-k2 *u(1); ...-2 *omega *cos ( psi ) *u (3) - k2 *u (2) ];u0 = [1 0 0 0]; [th,uh]=eulero_avanti(f,0,300,u0,h); figure(1); plot(uh(1,:),uh(2,:)) figure(2); plot(th,sqrt(uh(1,:).^2 + uh(2,:).^2)) Con pochi esperimenti si giunge alla conclusione che il metodo di Eulero esplicito non fornisce per questo problema soluzioni sicamente plausibili, neppure perhmolto piccolo. Ad esempio, in Figura 3 a sinistra viene riportato il graco, nel piano delle fasixy, dei movimenti del pendolo calcolati prendendoh= 1=1000. Come ci si aspetta il piano di rotazione del pendolo cambia al passare del tempo, ma, nonostante il passo di discretizzazione piccolo, l'ampiezza delle oscillazioni aumenta inaspettatamente, come si può vedere anche da Figura 4, a sinistra. Risultati analoghi si trovano anche per valori molto più piccoli dih. Ciò accade perché problemi come questo, che presentano soluzioni limitate pert! 1, ma non smorzate, hanno un comportamento analogo a quello del problema modello(y0 (t) =y(t); t >0; y(0) = 1;5 Figure 3  Traiettorie nel piano delle fasi per il pendolo di Foucault, ottenute con il metodo di Eulero esplicito (a sinistra) e conode45(a destra). Si noti che nel caso della soluzione ottenuta conode45la componente in direzioneyassume valori di due ordini di grandezza più piccoli rispetto a quella in direzionex.Figure 4  Andamento temporale del moduloj(x; y)jdel vettore posizione, ottenuto con il metodo di Eulero esplicito (a sinistra) e conode45(a destra). con valori dipuramente immaginari. In tal caso infatti la soluzione esatta è una funzione sinusoidale int. D'altra parte il metodo di Eulero esplicito ha una regione di assoluta stabilità tangente all'asse immaginario. Di conseguenza, solo il valoreh= 0garantirebbe assoluta stabilità. Osserviamo infatti che, se riscrivessimo il problema del pendolo nella forma y0 (t) =Ay(t); t >0 y0 (0) =y 0 con A=0 B B @0 0 1 0 0 0 0 1 k2 0 0 2!sin( ) 0k2 2!cos( ) 01 C C A; è immediato vericare, mediante il comandoeig, che questa matrice presenta quattro autovalori immaginari, con parte reale nulla. Per confronto, nelle Figure 3 e 4, a destra, abbiamo rappresentato la soluzione ottenuta con la funzione ode45di Matlab. Essa corrisponde ad un metodo Runge-Kutta adattivo che presenta una regione di assoluta stabilità che interseca l'asse immaginario. In eetti, sono sucienti 1089 passi (non uniformi) per ottenere la soluzione.6 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. Soluzione 3a)Per la scelta particolare dif(t; y) =y, gli stadi del metodo di Runge-Kutta proposto si riscrivono come segue : k1= u n k2=  un+h2 k 1 = +12 h 2 un k3=  un+34 hk 2 = +34 h 2 +38 h 2 3 un un+1= 1 +h+12 h 2 2 +16 h 3 3 un: Pertanto,g(z) = 1 +z+12 z2 +16 z3 . b)La regione di assoluta stabilità può essere tracciata usando i seguenti comandi :close all clear all clc %% %grigliaperplot rea = linspace(-3,3,301); img = linspace(-3,3,301); [a,b] = meshgrid(rea,img); %z z = a + i *b;%fattoredicrescita 7 Figure 5  Regioni di assoluta stabilità per Eulero Esplicito e diversi metodi di Runge-Kutta a due, tre e quattro stadi. A destra uno zoom in prossimità dell'asse immaginario.g = 1 + z + 0.5 *z.^2 + 1/6 *z.^3;gmag = abs(g); figure(1) contour(a,b,gmag,[1 1],r-); axis([rea(1),rea( end),img(1),img(end)]);axis(square); xlabel(Re(h\lambda)); ylabel(Im(h\lambda)); grid on; %% g = 1 + z + 0.5 *z.^2 ; %+1/6 *z.^3+1/24 *z.^4;gmag = abs(g); figure(2) hold on contour(a,b,gmag,[1 1],g-); hold off axis([rea(1),rea( end),img(1),img(end)]);axis(square); xlabel(Re(h\lambda)); ylabel(Im(h\lambda)); grid on; figure(1) axis([-0.2 0.2 -1.5 1.5]) legend(EE(RK1),RK2,RK3,RK4,Location,East) Tale regione è riportata in Figura 5, insieme alle regioni di assoluta stabilità di altri metodi di Runge- Kutta, ottenute in maniera analoga a quanto fatto al punto a. Si noti che nei casi a 1 e 2 stadi la regione di stabilità è tangente all'asse immaginario, mentre all'aumentare del numero di stadi, essa racchiude un intervallo più ampio di valori puramente immaginari. Come visto nell'esercizio precedente, ciò signica che diversi metodi di Runge-Kutta a più di due passi permettono di risolvere problemi dierenziali in cui la soluzione presenta oscillazioni non smorzate.8 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 ;(5) 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 ?9 Soluzione 4 Discretizzando in spazio la formulazione debole del problema con il metodo di Galerkin-elementi niti, otte- niamo un sistema (lineare, a coecienti costanti) di equazioni dierenziali ordinarie : M@ u@ t + A( t)u=f(t): Quest'ultimo viene poi approssimato in tempo con il-metodo. Otteniamo, in conclusione, il seguente pro- blema algebrico da risolvere a ogni passo di tempo, Mu n +1 un
t+ An +1 un +1 + (1)An un =fn +1 + (1)fn ; vale a dire M
t+ An +1 |{z} Cu n +1 =fn +1 + (1)fn +M
tu n (1)An un |{z} b: La funzioneParabolic1D_Solver.mrisolve il seguente problema ai valori ai limiti e iniziali (problema evo- lutivo) :@ u@ t @@ x  @ u@ x  +@ u@ x + u=f ;in = (a; b)(t 0; T ]; dotato di opportune condizioni al contorno (di tipo Dirichlet, Neumann o Robin) :u=g;su D; (Dirichlet) @ u@ x = h;su N; (Neumann) u(x) + @ u@ x = r;su R; (Robin) e condizioni iniziali :u(x;0) =u 0: Per risolvere un problema parabolico conPoliFEMoccorre specicare, oltre ai coecienti dell'operatore dierenziale e alle condizioni al contorno, anche la condizione iniziale e i dati relativi alla discretizzazione in tempo (passo di tempo, coeciente, intervallo di tempo) richiesti dal-metodo. Un template generale per il datale è fornito nel seguente script :%BCsflags(admissiblevaluesare1and2) data.flag_dirichlet = []; %$\Gamma_D$data.flag_neumann = []; %$\Gamma_N$data.flag_robin = []; %$\Gamma_R$%BCsfunctionsandparameters data.bcDir = @(x,t,param)(0. *x); %$g$data.bcNeu = @(x,t,param)(0. *x); %$h$data.bcRob_c1 = @(x,t,param)(0. *x); %$c_1$data.bcRob_c2 = @(x,t,param)(0. *x); %$c_2$data.bcRob_fun = @(x,t,param)(0. *x); %$r$%Forcingterm data.force = @(x,t,param)(0. *x); %$f(x)$%Differentialoperatorparameters data.diffusion = @(x,t,param)(0. *x); %$\alpha(x)$or$\mu(x)$data.transport = @(x,t,param)(0. *x); %$\beta(x)$data.reaction = @(x,t,param)(0. *x); %$\gamma(x)$%Exactsolutionandderivative(optionals) data.uexact = @(x,t,param)(0. *x); %$u_{ex}(x)$data.uxexact = @(x,t,param)(0. *x); %$u_{ex}(x)$10 %Stabilization(emptyorGLS) data.stabilization = []; %Forparabolicproblems data.u0 = @(x,t,param)(0. *x); %initialconditiondata.time.t0 = 0; %initialtimedata.time.tf = 1; %finaltimedata.time.dt = 0.1; %timestepdata.time.theta = 0.5; %parameterthetaintheTheta-Methoda)Iniziamo a generare il datale per il primo caso, dopo aver calcolato la derivata prima (rispetto a x) della soluzione esatta, richiesta per il calcolo degli errori.data.flag_dirichlet = [1 2]; data.flag_neumann = []; data.flag_robin = []; data.bcDir = @(x,t,param)(1./sqrt(4 *pi *0.1 *t) *exp(-x.^2./(4 *0.1 *t)) + 0. *x);data.bcNeu = @(x,t,param)(0. *x);data.bcRob_c1 = @(x,t,param)(0. *x);data.bcRob_c2 = @(x,t,param)(0. *x);data.bcRob_fun = @(x,t,param)(0. *x);data.force = @(x,t,param)(0 + 0. *x);data.diffusion = @(x,t,param)( 0.1 + 0. *x);data.transport = @(x,t,param)(0 + 0. *x);data.reaction = @(x,t,param)(0 + 0. *x);data.uexact = @(x,t,param)(1./sqrt(4 *pi *0.1 *t) *exp(-x.^2./(4 *0.1 *t)));data.uxexact = @(x,t,param)((-x./(2 *0.1 *t)). *(1./sqrt(4 *pi *0.1 *t) ...* exp(-x.^2./(4 *0.1 *t))));data.u0 = @(x,t,param)(1./sqrt(4 *pi *0.1 *t) *exp(-x.^2./(4 *0.1 *t)));data.time.t0 = 0.5; data.time.tf = 1; data.time.dt = param(1); data.time.theta = param(2); Il mainle per la chiamata della funzioneParabolic1D_Solver.mrisulta il seguente, nel caso del metodo di Eulero all'indietro :close all clear all clc addpath(Code) a = 3; N = [6/0.05 6/0.025 6/0.0125]; dt = [0.05 0.025 0.0125 ]; %dt=[0.5 *0.00250.5 *6.25e-40.5 *1.5625e-4];Nt = 0.5./dt; errorLinfL2 = zeros(1,length(dt)) ; errorLinfH1 = zeros(1,length(dt)) ; errorL2H1 = zeros(1,length(dt)) ; %%BackwardEuler fem =P1; Theta = 1; for i = 1 : length(dt)nx = N(i); 11 [vertices,elements,boundaries] = mesh1D(-a, a , nx); [solution, ~, ~, DATA, errorLinfL2(i), errorLinfH1(i), errorL2H1(i)] = ...Parabolic1D_Solver(elements,vertices,boundaries,fem,Eser20_time_dependent_heat_data,[dt(i) Theta]); %--Toplotthesolution,uncommentthefollowinglines figure(1) plot(vertices, solution(:, end));title(BackwardEuler,finaltime) for j = 1 : DATA.time.Ntfigure(3) feplot(solution(:,j+1),vertices,elements,fem); axis([-3 3 0 2]) time_plot = DATA.time.t0+j *dt(i);title([time=,num2str(time_plot)]) pause(0.25) end end figure loglog(dt, errorLinfH1,-or,LineWidth,2) hold on loglog(dt,dt/dt(1) *errorLinfH1(1),--k,LineWidth,1)grid on legend(ErrorLinfH1,dt) title(BackwardEuler) Si ottiene, come si può vedere dalla Figura 6, che l'errore converge con ordine 1, e non si hanno limiti nella scelta di
trispetto adh.Figure 6  Andamento degli errori di discretizzazioneku ex u hk L1 (t 0;T ;H1 (a;a))rispetto al passo di tempo
t, metodo di Eulero all'indietro (o Eulero implicito) b)Nel caso del metodo di Eulero in avanti, dimezzando
tquando si dimezzah, si ottiene un metodo che non risulta asintoticamente stabile. Infatti, il-metodo per 0. Anche se il dato iniziale è discontinuo in qualche punto, immediatamente dopo la soluzione è diventata continua e anzi dotata di derivate di ogni ordine. La diusione è quindi un processo regolarizzante, che tende cioè a smussare le irregolarità. Il codice usato per produrre questo risultato è il seguente :close all clear all clc a = 3; N = 6/0.05; Nt = 1/0.1; Theta = [1]; fem =P1; for i = 1 : length(Theta)[vertices,elements,boundaries] = mesh1D(-a, a , N); [solution,~,~,DATA] = Parabolic1D_Solver(elements,vertices,boundaries, ...fem,EX8_parabolic_data2,Theta(i)); %Timeplot deltat = DATA.time.dt; for j = [0:DATA.time.Nt]figure(33+i) feplot(solution(:,j+1),vertices,elements,fem); axis([-3 3 0 2]) time_plot = DATA.time.t0+j *deltat;grid on title([Theta=, num2str(Theta(i)),.Time=,num2str(time_plot)]) pause(0.5) end end 140.015 0.02 0.025 0.03 0.035 0.04 0.045 0.05 10-4 10 -3 10 -2 Crank-Nicolson Error Linf H1 dt 2