Mathematical Engineering - Matematica Numerica

Exercise - 17 - Solution

Divided by topic

MATEMATICA NUMERICA A.A. 2018 - 2019 Ingegneria Matematica Prof. A. Quarteroni Prof. A. Manzoni, Dr. I. Fumagalli Esercitazione 17 - Soluzione Equazioni dierenziali ordinarie (2) Esercizio 1 Si consideri il problema di Cauchy (y0 (t) =f(t; y) =ty(t)2 t0< t t max y(t 0) = 2(1) La soluzione esatta di tale problema, nell'intervallo limitatot2[0;3]è y(t) =2t 2 + 1: a)Scrivere la functioneulero_avantiche implementa il metodo di Eulero in avanti per la risoluzione di un generico problema di Cauchy. La funzioneeulero_avantiriceve in ingresso : la funzionef(t; y)denita come anonymous function (tramite il comando@). Si noti chefè una funzione di due variabili ; gli estremit 0e t maxdell'intervallo di denizione ; il dato inizialey 0; il passo di avanzamento temporaleh. La funzioneeulero_avantirestituisce il vettorethdegli istanti temporali discreti e il vettoreuhdei valori della corrispondente soluzione approssimata. L'intestazione dieulero_avantisarà quindi :function [th,uh]=eulero_avanti(f,t_0,t_max,y0,h)b)Utilizzando la function appena scritta, risolvere il problema (1) con il metodo di Eulero in avanti, scegliendo per i valori del passo temporaleh 1= 3 =6eh 2= 3 =20. Rappresentare sullo stesso graco le soluzioni numeriche ottenute e confrontarle con la soluzione esatta. c)Scrivere la functioneulero_indietroche implementa il metodo di Eulero all'indietro per la risoluzione di un generico problema di Cauchy, con intestazionefunction [th,uh,iter_pf]=eulero_indietro(f,t_0,t_max,y0,h);La funzione eulero_indietrorestituisce il vettorethdegli istanti temporali discreti e il vettoreuh dei valori della corrispondente soluzione approssimata. Ad ogni passo temporale, il calcolo di un +1 =un +hf(tn +1 ; un +1 ) richiede la soluzione di un'equazione non lineare tramite un metodo di punto sso. La funzione eulero_indietrorestituirà quindi anche il vettoreiter_pfche contiene il numero delle iterazioni eettuate ad ogni passo temporale, dal metodo di punto sso, per risolvere l'equazione non lineare. Si scelgano, per le iterazioni di punto sso, una tolleranza sul criterio d'arresto basato sull'incremento pari a10 5 e un numero massimo di iterazioni pari a1000.1 d)Utilizzando la function appena scritta, risolvere il problema (1) con il metodo di Eulero all'indietro, scegliendo per i valori del passo temporaleh 1= 3 =6eh 2= 3 =20. Rappresentare sullo stesso graco le soluzioni numeriche ottenute e confrontarle con la soluzione esatta. e)Riportare su un graco l'andamento del numero di iterazioni impiegate dal metodo di punto sso perrisolvere l'equazione non lineare in funzione degli istanti temporali perh 1= 3 =6eh= 3=20. f )Denito eh= max n=0:N hj y(t n) u nj ; il massimo modulo dell'errore compiuto approssimando la soluzione esattay(t n) con la soluzione numericau n, riportare, su un graco in scala logaritmica, l'andamento di e hal variare di h, utilizzando i passi temporalih= [1=4;1=8;1=16;1=32;1=64]. Soluzione 1a)Una possibile implementazione è la seguente :function [th,uh]=eulero_avanti(f,t_0,t_max,y0,h)%[th,uh]=eulero_avanti(f,t_max,y0,h) % %RisolveilproblemadiCauchy % %y=f(t,y) %y(t_0)=y_0 % %utilizzandoilmetododiEuleroinAvanti %(EuleroEsplicito):u^(n+1)=u^n+h *f^n% %Input: %->f:functionchedescriveilproblemadiCauchy %(dichiarataadesempiotramiteinlineo@). %Devericevereiningressodueargomenti:f=f(t,y) %->t_0:listanteinizialedellintervallotemporaledisoluzione %->t_max:listantefinaledellintervallotemporaledisoluzione %->y0:ildatoinizialedelproblemadiCauchy %->h:lampiezzadelpassodidiscretizzazionetemporale. % %Output: %->th:vettoredegliistantiincuisicalcolalasoluzionediscreta %->uh:lasoluzionediscretacalcolataneinoditemporalit_h %vettoredegliistantiincuirisolvolaedo th = t_0:h:t_max; N = length(th); %uh=zeros(1,N); %uh=[y0zeros(size(y0,1),N)]; uh(:,1) = y0; for it=2:Nuh(:,it) = uh(:,it-1) + h *f(th(it-1),uh(:,it-1));end b)La soluzione si ottiene con i seguenti comandi : y0=-2; t0 = 0; tmax=3; tdis=t0:0.01:tmax; y = @(t) -2./(t.^2+1); %risolvoconeuleroinavanti f = @(t,y) t *y.^2;h1 = 3/6; 2 [EA_th1,EA_uh1] = eulero_avanti(f,t0,tmax,y0,h1); h2 = 3/20; [EA_th2,EA_uh2] = eulero_avanti(f,t0,tmax,y0,h2); figure; plot(tdis,y(tdis),k); xlabel(t); ylabel(y(t)); hold on; plot(EA_th1, EA_uh1,bo-,EA_th2, EA_uh2,ro-,Linewidth,2); legend(Soluzioneanalitica,EA,h=3/6,EA,h=3/20,Location,best); Il graco è mostrato in Figura 1. Figure 1  Confronto tra la soluzione esatta e le soluzioni numeriche ottenute con il metodo di Eulero in avanti, con passi di discretizzazioneh 1= 3 =6eh 2= 3 =20. c)Una possibile implementazione è la seguentefunction [th,uh,iter_pf]=eulero_indietro(f,t_0,t_max,y0,h)%[th,uh,iter_pf]=eulero_indietro(f,t_max,y0,h) %RisolveilproblemadiCauchy % %y=f(y,t) %y(t_0)=y_0 % %utilizzandoilmetododiEuleroIndietro %(EuleroImplicito):u^{n+1}=u^n+h *f^{n+1}.%Perciascunistantetemporalelequazioneperu^(n+1) %(apriorinonlineare)vienerisoltaconunmetododipuntofisso: %u=u^n+delta_t *f(t^{n+1},u),utilizzandolafunzionediiterazione% %phi(x)=u^n+h *f(t^{n+1},x)% %lacondizioneperlaconvergenzadelmetododipuntofisso(|phi(x)|t_max:listantefinaledellintervallotemporaledisoluzione %(listanteinizialeet_0=0) 30 0.5 1 1.5 2 2.5 3 -2 -1.8 -1.6 -1.4 -1.2 -1 -0.8 -0.6 -0.4 -0.2 0 t y(t) Soluzione analitica EA, h=3/6 EA, h=3/20 %->y0:ildatoinizialedelproblemadiCauchy %->h:lampiezzadelpassodidiscretizzazionetemporale. % %Output: %->th:vettoredegliistantiincuisicalcolalasoluzionediscreta %->uh:lasoluzionediscretacalcolataneinoditemporalit_h %->iter_pf:vettorechecontieneilnumerodiiterazioni %dipuntofissoutilizzateperrisolverelequazione %nonlineareadogniistantetemporale. th = t_0:h:t_max; N = length(th); %uh=zeros(1,N); %uh(1)=y0; uh(:,1) = y0; %parametriperleiterazionidipuntofisso Nmax = 1000; toll = 1e-5; iter_pf = zeros(1,N); for it = 2:N%ciclointempophi = @(u) uh(:,it-1) + h *f( th(it), u ); %mappaphi_nj = 0; err = toll+1; w0 = uh(:,it-1); while err>toll && j0e quindi la soluzioneysarebbe crescente in un intorno di t: non è dunque possibile cheysi annulli dal momento chey(0) = 1.80 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 -0.5 0 0.5 1 1.5 2 soluzione analitica soluzione di EA con h 2 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 soluzione analitica soluzione di EA con h 3 Dal calcolo della derivata, fy= @ f =@ y= 3=(1 + 9y2 )3 si ha3< f y< 0 per ogniy >0, da cui possiamo scegliere max= max jf yj = 3: Possiamo quindi aspettarci che le perturbazioni nel metodo di Eulero in avanti restino controllate pur- chéh