Mathematical Engineering - Matematica Numerica

Exercise - 12 - Solution

Divided by topic

MATEMATICA NUMERICA A.A. 2018 - 2019 Ingegneria Matematica Prof. A. Quarteroni Prof. A. Manzoni, Dr. I. Fumagalli Esercitazione 12 - SoluzioneEquazioni non lineari (2) Esercizio 1 Si considerino le due equazioni seguenti, nell'intervallo[0:5;2]: f(x) = ln(x) = 0;(1) g(x) = 1x+ ln(x) = 0:(2) a)Dal graco dif ; ge delle loro derivate, discutere le proprietà di convergenza del metodo di Newton per tutti gli zeri, valutando l'opportunità di applicare il metodo di Newton modicato. Noto il valore esatto delle radici, se ne valutino le molteplicità. b)Considerando lo sviluppo di Taylor al second'ordine difegin un intorno dix= 1, trovare gli esponentipeqtali che f(x) = (1x)p F(x);(3) g(x) = (1x)q G(x);(4) conF; Gfunzioni tali cheF(1)6 = 0; G(1)6 = 0. Basandosi su tale risultato, provare che il metodo di Newton per l'approssimazione della radice= 1può avere ordine 2 solo se applicato af. c)Per la funzioneg, provare che il metodo modicato x( k+1) =x( k) 2g (x( k) )g 0 (x( k) ) converge localmente ad= 1con ordine 2. d)Scrivere una funzione Matlabnewton.mche implementi il metodo di Newton ; la funzione riceve in input la guess inizialex0, la tolleranzatoll, il numero massimo di iterazioninmax, la funzionefun, la sua derivatadfune la molteplicità della radicemol. restituisce in output il vettorex_vectdelle iterate e il numeroitdi iterazioni eettivamente eseguite. Si utilizzi un criterio d'arresto basato sul modulo della dierenza tra due iterate successive. La function newton.mdeve comportarsi come il metodo di Newton classico o come il metodo di Newton modicato x( k+1) =x( k) molf (x( k) )f 0 (x( k) ) a seconda del valore della molteplicità passata come argomento. L'intestazione della funzione sarà ad esempio la seguente : [x_vect,it]=newton(x0,nmax,toll,fun,dfun,mol) e)Risolvere separatamente i problemi della ricerca degli zeri dife digtramite il metodo di Newton, utilizzando la funzionenewton.mcon tolleranza10 6 enmax= 100. Si utilizzi per il calcolo di ogni radice un opportuno valorex0. Si riporti su un graco in scala semilogaritmica l'andamento degli errori in funzione del numero di iterazioni per il metodo di Newton standard e per quello modicato.1 Figure 1  Graco delle funzionif(x)eg(x)in blu, delle loro derivatef0 (x)eg0 (x)in rosso e della retta y= 0in nero tratteggiato. Soluzione 1a)Supponiamo di inserire tutti i seguenti comandi, necessari per risolvere l'esercizio, in uno script diMatlab. La denizione dif(x); g(x)ef0 (x); g0 (x)come@function e i loro graci si ottengono con i seguenti comandi :f = @(x) log(x); g = @(x) 1 - x + log(x); df = @(x) 1./x; dg = @(x) -1 + 1./x; x = linspace(0.5, 2, 100); figure(1); plot(x,f(x), x,df(x), x,0. *x,k--)legend(f,f); xlabel(x); grid on figure(2) plot(x,g(x), x,dg(x), x,0. *x,k--)legend(g,g); xlabel(x); grid on I graci risultanti sono mostrati in Figura 1. Osserviamo che la radice = 1è una radice semplice per f, mentre è doppia perg, in quanto ancheg0 () = 0. Utilizzando il metodo di Newton classico per calcolareci aspettiamo di avere una convergenza quadratica perfe lineare perg. Verichiamo la molteplicità della radice= 1tramite i seguenti comandi :%discussionesulleradicidalgraficodifedf d2f = @(x) -1./x.^2; alpha = 1; if (df(alpha) == 0)if (d2f(alpha) == 0)disp(laradicealphahamolteplicitamaggioredidue) else disp(laradicealphahamolteplicitaugualeadue) end else disp(laradicealphahamolteplicitaugualeaduno) end d2g = @(x) -1./x.^2; if (dg(alpha) == 0)if (d2g(alpha) == 0)disp(laradicealphahamolteplicitamaggioredidue) else 20.511.52x-1 -0.5 0 0.5 1 1.5 2ff' 0.511.52x -0.5 0 0.5 1gg' disp(laradicealphahamolteplicitaugualeadue) end else disp(laradicealphahamolteplicitaugualeaduno) end e otteniamo la seguente stampa a schermo : la radice alpha ha molteplicita uguale a uno la radice alpha ha molteplicita uguale a due b)Poiché f(1) =g(1) = 0e f0 (x) =1x ; f 00 (x) =1x 2; le funzionifegammettono i seguenti sviluppi di Taylor, in un intorno dix= 1: f(x) =x112 ( x1)2 +o (x1)2
g(x) =12 ( x1)2 +o (x1)2
e dunque abbiamop= 1eq= 2, conF(x) =112 ( x1) +o (x1)2
eG(x) =12 + o (x1)2
. Per valutare l'ordine di convergenza del metodo di Newton, consideriamo inizialmente l'applicazione di tale metodo alla funzioneg: x( k+1) = N( x( k) ) =x( k) g (x( k) )g 0 (x( k) ): La derivata della funzione di iterazione associata al metodo di Newton sarà dunque 0 N( x) = 1g 0 (x)2 g(x)g00 (x)g 0 (x)2=g (x)g00 (x)g 0 (x)2: Poichég(x) = (x1)2 G(x); g0 (x) = 2(x1)G(x) + (x1)2 G0 (x); g00 (x) = 2G(x) + 4(x1)G0 (x) + (x1)2 G00 (x); Sostituendo nell'espressione di0 N( x)otteniamo 0 N( x) =( x1)2 G(x) 2G(x) + 4(x1)G0 (x) + (x1)2 G00 (x)[2( x1)G(x) + (x1)2 G0 (x)]2 ; da cui si ricavalim x!1 0 N( x) =2 G(x)24 G(x)2=12 : Questa relazione dimostra che, nel caso in esame, il metodo di Newton converge solo linearmente, con coeciente di riduzione lineare dell'errore dato da j0 N(1) j=12 : Se ora applichiamo il metodo di Newton alla funzionef, la funzione di iterazione risulta essere N( x) =xf (x)f 0 (x):3 Ripercorrendo i passaggi di cui sopra, possiamo mostrare che 0 N( x) =( x1)F(x) [2F0 (x) + (x1)F00 (x)][ F(x) + (x1)F0 (x)]2 : Poiché il denominatore non si annulla inx= 1, possiamo calcolare 0 N(1) : dal momento che otteniamo 0 N(1) = 0 , in questo caso il metodo di Newton risulta di ordine superiore al primo (in particolare, dalla teoria sappiamo che sarà di ordine 2). c)Il metodo modicato ha una funzione di iterazione denita da M( x) =x2g (x)g 0 (x)= 2  N( x)x e dunque0 M( x) = 20 N( x)1; da cuilim x! 0 M( x) = 2
12 1 = 0; ovvero il metodo modicato ha ordine 2. d)Una possibile implementazione della funzione richiesta è la seguente :function [x_vect,it] = newton(x0,nmax,toll,fun,dfun,mol)%[x_vect,it]=newton(x0,nmax,toll,fun,dfun,mol) % %MetododiNewtonperlaricercadeglizeridellafunzionefun. %Testdarrestobasatosullincrementotradueiterazionisuccessive. % %Parametridiingresso: % %x0Puntodipartenza %nmaxNumeromassimodiiterazioni %tollTolleranzasultestdarresto %fundfunAnonymousfunctioncontenentilafunzioneelasuaderivata %molMolteplicitadellozerochesivuoletrovare % %Parametridiuscita: % %x_vectVett.contenentetutteleiteratecalcolate %(lultimacomponenteelasoluzione) %itIterazionieffettuate x_vect = x0; err = toll+1; it = 0; x_old = x0; while it < nmax && err > tollfx = fun(x_old); dfx = dfun(x_old); if dfx == 0%errorfaterminarelesecuzionedellafunction error(Arrestoperazzeramentodidfun); end %nuovaiterazionedinewton x_new = x_old-mol *fx/dfx;err = abs(x_new-x_old); %salvolequantitadiinteresseneglioutput x_vect = [x_vect; x_new]; it = it+1; 4 x_old = x_new; end fprintf(NumerodiIterazioni:%d\n, it); fprintf(Radicecalcolata:%-12.8g\n, x_new); return e)Con i seguenti comandi si ottengono le operazioni richeste : figure(12) plot(x,g(x), x,dg(x), x,d2g(x), x,0. *x,k--)legend(g,g,g); xlabel(x); grid on %%1.e toll = 1e-6; nmax = 100; disp(Calcolodellaradiceconmolteplicita1) x0_f = 0.5; [xvect_f,it_f]=newton(x0_f,nmax,toll,f,df,1); disp(Calcolodellaradiceconmolteplicita2) x0_g = 0.5; disp(MetododiNewtonsemplice); [xvect_g,it_g]=newton(x0_g,nmax,toll,g,dg,1); disp(MetododiNewtonmodificato); [xvect_gm,it_gm]=newton(x0_g,nmax,toll,g,dg,2); %% alpha = 1; %soluzioneesattaerr_f = abs(xvect_f - alpha); err_g = abs(xvect_g - alpha); err_gm = abs(xvect_gm - alpha); figure(13); semilogy(0:it_f,err_f,ob-); grid on title(Radicesemplice) Il precedente codice produce a schermo : Calcolo della radice con molteplicita 1 Convergenza al passo k : 5 Radice calcolata : 1.00000000 Calcolo della radice con molteplicita 2 Metodo di Newton semplice Convergenza al passo k : 20 Radice calcolata : 0.99999926 Metodo di Newton modificato Convergenza al passo k : 5 Radice calcolata : 1.00000000 I graci dell'errore sono mostrati in Figura 2 : nel caso di radice doppia, il metodo di Newton modicato permette di recuperare l'andamento quadratico dell'errore.5 Figure 2  Andamento dell'errore del metodo di Newton nel caso di radice semplice (funzionef, a sinistra) e confronto con il metodo modicato nel caso di radice doppia (funzioneg, a destra). Esercizio 2 Si vogliono calcolare le radici 1= 2 e 2= 3 dell'equazionef(x) =x2 5x+ 6. a)A partire dall'espressione dif, proporre due metodi di punto sso per il calcolo di 1e  2. b)Analizzare la consistenza e la convergenza dei due metodi.c)Scrivere unafunctionche implementi il metodo di punto sso, con la seguente rma[xvect, niter] = fixed_point_iterations( phi, x0, tol, kmax ) dove x0è il punto di partenza,kmaxetolil numero massimo di iterazioni e la tolleranza per il criterio d'arresto,phila funzione di iterazione. In uscita, la funzione deve restituire il vettore delle iteratexvecte il numeroniterdi iterazioni eettuate. d)Vericare sperimentalmente le conclusioni tratte al punto b) utilizzando la funzionefixed_point_iterations con tolleranza10 10 e variando il dato inizialex(0) =1.5, 2.5, 2.98, 3.1, per approssimare 1. Ripe- tere considerandox(0) =2.5, 3.5, 2.1, 1.9 per approssimare 2. A partire dai vettori delle iterate, stimare anche l'ordine di convergenza con la funzioneqssstimap, che calcola l'ordinepe il fattore di riduzioneCper una successione convergentex( k) assumendo che valga jx( k+1) x( k) j 'C x( k) x( k1) p : Soluzione 2 Le radici della funzionefsi calcolano facilmente e sono 1= 2 , 2= 3 . a)Innanzitutto ribadiamo che non ci sono ricette per costruire delle buone funzioni di iterazione. Unprimo algoritmo può essere ottenuto isolando, nell'espressione della funzionef(x), il termine5xrispetto agli altri termini : f(x) =x2 5x+ 6 = 0,5x=x2 + 6,x=x 2 + 65 ; che conduce al metodo iterativo associato alla funzione 1( x) =x 2 +65 , ossia x( k+1) =( x( k) )2 + 65 ; k 0: Un secondo algoritmo si può costruire risolvendo il termine quadraticox2 rispetto agli altri termini f(x) =x2 5x+ 6 = 0,x2 = 5x6,x=p5 x6; ove si è presa soltanto la radice positiva dato che sia 1che  2sono positive. La funzione di iterazione che si ottiene in questo caso è 2( x) =p5 x6alla quale corrisponde l'algoritmo : x( k+1) =p5 x( k) 6k0:611.522.533.54 iterazioni 10-1010-810-610-410-2100 errore 05101520 iterazioni 10-2010-1510-1010-5100 errore NewtonNewton modificato b) Consistenza Per vericare la consistenza dei metodi dobbiamo accertarci che la radicedifsia anche punto sso della funzione di iterazione: =(): Osserviamo che ogni qualvolta la funzione di iterazione è ottenuta tramite manipolazioni algebriche della funzionefla consistenza è garantita per costruzione. In particolare i due casi in esame sono consistenti per costruzione. Convergenza Applichiamo il Teorema di Convergenza Locale, che ci assicura la convergenza in un intorno del punto sso purché valgaj0 ()j 1: Si conclude che il metodo 1 converge al punto sso 1e non converge al punto sso  2. Per il secondo metodo, si ha0 2( x) =52 p5 x6, dunque 0 2(  1) =54 > 1; 0 2(  2) =56 2 (1;1): Il metodo 2 converge ad 2mentre non converge ad  1. Inoltre entrambi i metodi, quando convergono, hanno ordine di convergenza 1 essendo0 ()6 = 0. c)Verica sperimentale Verichiamo la convergenza confixed_point_iterationse determiniamo l'ordine di convergenza conqssstimap.toll = 1e-10; nmax = 100; phi_1 = @(x) (x.^2 + 6)/5; fun = @(x) x.^2 - 5 *x + 6;x0 = 1.5; [xvect,niter] = fixed_point_iterations(phi_1,x0,tol,kmax); subplot(2,2,1); plot(xvect,.-); title([x^{(0)}=,num2str(x0)]) xvect( end),niterans = 2.0000 niter = 92 [p,c] = qssstimap(xvect); Ordine stimato : 1.00000437 Fattore di riduzione : 0.80007770 x0 = 2.5; [xvect,niter] = fixed_point_iterations(phi_1,x0,tol,kmax); subplot(2,2,2); plot(xvect,.-); title([x^{(0)}=,num2str(x0)]) xvect( end),niterans = 2.0000 niter = 98 [p,c] = qssstimap(xvect); Ordine stimato : 1.00000000 Fattore di riduzione : 0.80000000 7 Figure 3  Graco delle iteratex( k) generate dalle funzioni d'iterazione 1(a sinistra) e  2(a destra), per diversi valori dix(0) .x0 = 2.98; nmax = 1000; [xvect,niter] = fixed_point_iterations(phi_1,x0,tol,kmax); subplot(2,2,3); plot(xvect,.-); title([x^{(0)}=,num2str(x0)]) xvect( end),niterans = 2.0000 niter = 119 [p,c] = qssstimap(xvect); Ordine stimato : 1.00001623 Fattore di riduzione : 0.80029612 x0 = 3.1; [xvect,niter] = fixed_point_iterations(phi_1,x0,tol,kmax); subplot(2,2,4); plot(xvect,.-); title([x^{(0)}=,num2str(x0)]) xvect( end),niterans = Inf niter = 25 Si verica pertanto che il metodo converge solo alla radice  1, che l'ordine di convergenza è 1 e che il fattore di abbattimento dell'errore è0 1(  1) = 4 =5 = 0:8. Si osserva inoltre che il punto 2è un repulsore. Ciò è confermato dagli andamenti delle iterate generate da 1e riportate Figura 3, in cui si nota anche che la convergenza (quando vericata) è monotona e che il metodo diverge per valori di x(0) >3: in 24 iterazioni raggiunge un valore talmente elevato chex(25) = 1( x(24) )superarealmax ed è dunque assimilato a+1da Matlab.phi_2 = @(x) sqrt(5 *x - 6);x0 = 2.5; [xvect,niter] = fixed_point_iterations(phi_2,x0,tol,kmax); subplot(3,2,1); plot(xvect,.-); title([x^{(0)}=,num2str(x0)]) xvect( end),niterans = 3.0000 80501001.41.61.8 2x(0)=1.5 0501002 2.22.42.6 x(0)=2.5 0501001502 2.5 3x(0)=2.98 01020300246810158x(0)=3.1 0501001502.62.8 3x(0)=2.5 0501001503 3.5 x(0)=3.5 0501001502 2.5 3x(0)=2.1 050100150024Re(x(k)), x(0)=1.9 050100150024Im(x(k)), x(0)=1.9 niter = 117 [p,c] = qssstimap(xvect); Ordine stimato : 0.99998671 Fattore di riduzione : 0.83308253 x0 = 3.5; [xvect,niter] = fixed_point_iterations(phi_2,x0,tol,kmax); subplot(3,2,2); plot(xvect,.-); title([x^{(0)}=,num2str(x0)]) xvect( end),niterans = 3.0000 niter = 112 [p,c] = qssstimap(xvect); Ordine stimato : 0.99999579 Fattore di riduzione : 0.83325376 x0 = 2.1; [xvect,niter] = fixed_point_iterations(phi_2,x0,tol,kmax); subplot(3,2,3); plot(xvect,.-); title([x^{(0)}=,num2str(x0)]) xvect( end),niterans = 3.0000 niter = 128 [p,c] = qssstimap(xvect); Ordine stimato : 0.99998923 Fattore di riduzione : 0.83312650 x0 = 1.9; [xvect,niter] = fixed_point_iterations(phi_2,x0,tol,kmax); subplot(3,2,4); plot(real(xvect),.-); title([Re(x^{(k)}),x^{(0)}=,num2str(x0)]) subplot(3,2,6); plot(imag(xvect),.-); title([Im(x^{(k)}),x^{(0)}=,num2str(x0)]) xvect( end),niterans = 3.0000 + 0.0000i niter = 129 [p,c] = qssstimap(xvect); Ordine stimato : 0.99999686 Fattore di riduzione : 0.83327328 Il secondo metodo converge in tutti i casi alla radice  2e non converge mai ad  1. Inoltre l'ordine di convergenza è 1 e il fattore di abbattimento dell'errore è circa quello previsto teoricamente,0 2(  2) = 5=6'0:8333. Ancora una volta, i graci di Figura 3 confermano quanto detto e mostrano una convergenza monotona nei primi tre casi. Notiamo inne che, perx(0) = 1:9, il metodo di punto sso si muove sulla retta reale solo perk= 0; : : : ;8, mentre poi le iterate hanno una parte immaginaria non9 nulla no a convergenza : Im(x(129) )'2:971
10 10 . Questo perché, durante le iterazioni, si arriva ad un puntox(8) '0:9369 tollxv = x; Jxv = J(xv); if det(Jxv) == 0error(ArrestoperJacobianosingolare); end x = xv - Jxv \ fun(xv); err = norm(x - xv, Inf); it = it + 1; end fprintf(NumerodiIterazioni:%d\n, it); fprintf(Zerocalcolato:%-12.13f\n, x(1)); for i=2:length(x)fprintf(%-12.13f\n, x(i)); end return c)La funzione newtonsys.mpuò essere utilizzata nel seguente modo :fun = @(x) [x(1)^2-4 *x(2)^2+3; x(1)-x(2)^2];11−1 0 1 2 3 4 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 J = @(x) [2 *x(1) (-8 *x(2)); 1 (-2 *x(2))];max_iter = 1000; tol = 1e-15; x_0 = [0; -0.5]; [xvect,it] = newtonsys(x_0, max_iter, tol, fun, J); In questo caso, l'output visualizzato nella command window sarà : Numero di Iterazioni : 6 Zero calcolato : 1.0000000000000 -1.0000000000000 d)Come noto dalla teoria, la convergenza del metodo di Newton ad un particolare zero della funzione f si ottiene scegliendo comex 0un punto sucientemente vicino alla soluzione cercata. Per determinare quindi le altre intersezioni fraC 1e C 2è possibile scegliere dei punti x 0vicini all'intersezione cercata, sfruttando ad esempio la rappresentazione graca costruita al punto 1 dell'esercizio. Un modo possibile è :x_0 = [3.5;-1.5]; [xvect,it] = newtonsys(x_0, max_iter, tol, fun, J); x_0 = [3.5; 1.5]; [xvect,it] = newtonsys(x_0, max_iter, tol, fun, J); x_0 = [0; 0.5]; [xvect,it] = newtonsys(x_0, max_iter, tol, fun, J); L'output sarà : Numero di Iterazioni : 6 Zero calcolato : 3.0000000000000 -1.7320508075689 Numero di Iterazioni : 6 Zero calcolato : 3.0000000000000 1.7320508075689 Numero di Iterazioni : 6 Zero calcolato : 1.0000000000000 1.0000000000000 Esercizio 4 Si consideri l'equazione f(x) =xarctan (ex 1) = 0; che ammette due soluzioni >1e= 0. Si consideri il metodo di punto sso seguente : x( k+1) =(x( k) ); (x) = arctan (ex 1): a)Dimostrare che esiste una costante positivaC 10 (x)è decrescente : in particolare 00 (x)K=0 (1)'0:68778x2[1;+1):(6) Si noti che, sex1, allora(x)1; dunquemanda[1;+1)in sé. Inoltre si ha x( k+1) =(x( k) )() =0 ( k)( x( k) );(7) ove kè compreso fra x( k) e, dunque è ancora in[1;+1). Pertanto, sex(0) >1, si ha che le successioni di valorix( k) e ksono tutte contenute in [1;+1). Dalla (7) e grazie alla (6) abbiamo dunque jx( k+1) j Kjx( k) j 8k= 0;1; : : : da cui la stima di convergenzajx( n) j Kn jx(0) j:Figure 5  Graci delle curve (A) :y=0 (x), (B) :y=x, (C) :y=(x)130 1 2 3 4 5 0 1 2 3 4 5 (A) (B) (C) Dalla Figura 5 si osserva come 2(1;2). Dunque, sex(0) = 1, si ha jx( n) j Kn ; da cui l'errore sarà certamente minore di10 8 pur di prendere n > 8 log 10K' 49:2: Eettuiamo quindi 50 iterazioni di punto sso :phi = @(x) atan(exp(x)-1); x = 1; fori = 1:50; x = phi(x);endx x = 1.1195 Otteniamo il valore x( n) = 1:1195(pern= 50). Sebbene non si conosca esattamente, grazie alla stima a priori sull'errore siamo sicuri che la precisione è almeno pari a 8 cifre decimali. Alternativamente, in mancanza di una stimaa prioridell'errore, potremmo utilizzare una stimaa posteriori(ovvero basata su quantità direttamente calcolabili, come l'incremento o il residuo ad ogni iterazione) come discusso nel punto c). Inne, il metodo di punto sso non è adeguato per approssimare= 0in quanto, essendo0 () = 1, il metodo non può essere convergente. b)Metodo di Newton : x( k+1) =x( k) f (x( k) )f 0 (x( k) ); ove nel nostro casof(x) =x(x). Dato chef0 ()>0,f0 () = 0, si ha cheè radice semplice dif mentreè radice doppia. Sappiamo allora che il metodo di Newton converge localmente per entrambi gli zeri, tuttavia l'ordine di convergenza sarà pari a2per, mentre sarà pari a 1 per. c)Il criterio d'arresto per un generico metodo di punto sso è generalmente basato sull'incremento, ossia si termina il processo quandojx( k+1) x( k) j< ; oveè la tolleranza stabilita. A partire dalle stime seguenti, j(x( k) )x( k) j=j (x( k) )+x( k) j  jx( k) j j(x( k) )j=j x( k) j j(x( k) )()j=j x( k) j j0 ( k)( x( k) )j=jx( k) j j0 ( k) jjx( k) j; ove kè compreso fra x( k) eper il teorema del valor medio, assumendo che il criterio d'arresto basato sull'incremento sia soddisfatto, e ricordando chej0 ( k) j C, otteniamo jx( k) j 11 Cj x( k+1) x( k) j