Mathematical Engineering - Matematica Numerica

Exercise - 11 - Solution

Divided by topic

MATEMATICA NUMERICA A.A. 2018 - 2019 Ingegneria Matematica Prof. A. Quarteroni Prof. A. Manzoni, Dr. I. Fumagalli Esercitazione 11 - SoluzioneEquazioni non lineari Esercizio 1 Si consideri il problema del calcolo del valore2I= [0; ]tale che = 1 +12 sin : a)Utilizzando opportunamente la visualizzazione di funzioni in Matlab, si determini se il metodo dibisezione è applicabile per calcolare. In seguito, si derivi una stima teorica del numero di iterazioni del metodo di bisezione necessarie per approssimarecon una tolleranza inferiore a10 3 . Inne, si verichi medianteqssbisezche il numero eettivo di iterazioni rispetta la stima. b)Si consideri ora il metodo di punto sso seguente x( k+1) =(x( k) ) = 1 +12 sin x( k) ; x(0) 2I : Si provi che il metodo èglobalmente convergentesuI, ovvero chelim k!1x( k) =8x(0) 2I. Inoltre, si provi che l'errore di approssimazione soddisfa la stima jx( k) j Ck jx(0) j: Inne, si dia una stima della costanteCe la si usi per stimare il numero di iterazioni necessarie per approssimarecon una tolleranza inferiore a10 3 . c)Per il metodo di punto sso in b), si vuole adottare un criterio di arresto basato sull'incremento, ovverol'algoritmo si arresta non appena si abbia jx( k) x( k1) j ": Si mostri come questo criterio garantisca chex( k) soddis jx( k) j "1 C: Si fornisca il valore di"che occorre scegliere per calcolare una approssimazione dicon un errore minore ditoll= 10 3 . Implementare in Matlab il metodo di punto sso scrivendo una function che riceva in ingresso la funzione,x(0) ,"e un numero massimo di iterazioni, e restituisca in uscita il valorex( k) con il numero dikdi iterazioni eettuate,function [alpha,niter]=fixed_point(phi,x0,epsilon,kmax)La funzione deve utilizzare il criterio d'arresto descritto. Si testi il proprio codice per x(0) = 0, ve- ricando che l'errorejx( k) jsia minore della tolleranzatoll= 10 3 , e che il numero di iterazioni sia approssimativamente quello stimato (usarefsolve(f, 0)per trovare lo zero esatto di una funzionef).1 Soluzione 1 a)Consideriamo la funzionef(x) =x112 sin( x); abbiamo cheè uno zero dif. Tracciamone il graco, ad esempio come segue.%Tracciamoilgraficodellafunzione xx = linspace(0,pi,100); f = @(x) x - 1 - 0.5 *sin(x);plot(xx,f(xx)); grid on Visualizzando il graco di fosserviamo la funzione assume segni opposti agli estremi dell'intervallo ; è dunque possibile utilizzare il metodo di bisezione per calcolare uno zero dif. Osserviamo anche lo zero è unico (in quantofè monotona,f0 (x) = 112 cos( x)12 ). Ad ogni iterazione del metodo di bisezione, l'errore sulla soluzione viene dimezzato. In generale, avremo quindi che, partendo dal punto medio dell'intervallo considerato,x(0) =2 , si avrà jx( k) j j Ij2 k +1=2 k +1: Occorrerà dunqueklog 2(103 )1 =)k11: In eetti, si trova :%Stimanumeroiterazioni kminBisez = ceil(log2(pi *1000) - 1)kminBisez = 11 %Metodobisezione [xvect,xdif,fx,it] = qssbisez(0,pi,100,1e-3,f); Numero di Iterazioni : 11 Radice calcolata : 1.49869923 b)Si noti che è ovviamente punto sso di, ovvero=(). Inoltre, abbiamo che i)([0; ]) = [12 ;32 ] [0; ]; ii)j0 (x)j= 12 cos( x) 12 8 x2[0; ]; da cui deduciamo, per il teorema di convergenza globale dei metodi di punto sso, chelim k!1x( k) = 8x(0) 2I= [0; ]. In eetti, si ha che per ognik2Nesiste kfra ex( k) tale che jx( k+1) j=j(x( k) )()j=j0 ( k) jjx( k) j Cjx( k) j; conC=12 . Pertanto, si ha la stima jx( k) j 12 kj x(0) j 2 k: Si noti che la stima dell'errore è analoga a quella ottenuta per il metodo di bisezione. Dunque, ci aspettiamo che con 12 iterazioni l'approssimazione calcolata soddis sicuramente la tolleranza di10 3 . Bisogna dire, tuttavia, che il coeciente di riduzione dell'erroreasintoticonon èC=12 , che è una stima per eccesso valida su tutto l'intervalloI, ma piuttostoj0 ( k) j ' j0 ()j '0:09, pertanto molto più piccolo ; ovvero, il numero di iterazioni richieste potrebbe essere signicativamente minore (se prendessimo0:09come coeciente di riduzione dell'errore, troveremmo che bastano circa 3 iterazioni). c)Si ha che Cjx( k) j  jx( k+1) j  jx( k+1) x( k) +x( k) j  jx( k) j jx( k+1) x( k) j da cui, sfruttandoC epsilon && k < kmaxk=k+1; xnew=phi(x); err=abs(xnew-x); x=xnew; end alpha = x; niter=k; return e si ha : %Metodopuntofisso phi = @(x) 1 + 0.5 *sin(x);epsilon = 0.5e-3; x0 = 0; kmax = 100; [alpha,niter] = fixed_point(phi,x0,epsilon,kmax) alpha = 1.4987 iterazioni = 5 %Errore abs(alpha - fsolve(f,0)) ans = 6.3651e-06 Si vede come il numero eettivo di iterazioni sia compreso fra 11 (stima per eccesso) e 3 (stima ottenuta considerando il fattore di abbattimento dell'errore asintotico, perksucientemente grande). Inoltre, l'errore è entro la tolleranza stabilita. Esercizio 2a)Datoa >0, scrivere il metodo di Newton per l'approssimazione della soluzione dell'equazione ln(ax) = 0:3 b)Si dimostri che nel caso considerato vale jx( k+1) 1=aj< Mjx( k) 1=aj2 , pur di prenderex(0) su- cientemente vicino a1=a, per una costante opportunaM >0. c)A partire da uno sviluppo di Taylor dig(x) =xln(ax)centrato in1=a, costruire uno schema di Newton modicato che sfrutti solamente le operazioni di somma-sottrazione e moltiplicazione per approssimare il reciproco di un numero positivoa. Cosa si può dire sulla convergenza di tale metodo ? Soluzione 2a)Cerchiamo lo zero= 1=adella funzionef(x) = ln(ax). Essendof0 (x) = 1=x, il metodo di Newton si scrive come segue : x( k+1) =x( k) f (x( k) )f 0 (x( k) )= x( k) x( k) ln(ax( k) ): b)Consideriamo lo sviluppo in serie dig(x) =xln(ax)centrato in1=ae troncato al secondo ordine. Essendo g(1=a) = (1=a) ln(1) = 0; g0 (x) = ln(ax) + 1; g0 (1=a) = 1; g00 (x) = 1=x; otteniamo : xln(ax) = x1a  +12  x1a  2 ;con2[x;1=a]punto di Lagrange: Sostituendo l'espressione trovata nel metodo di Newton si ha,x( k+1) 1a = ( x( k) 1a ) (x( k) 1a ) 12  k( x( k) 1a ) 2 ; ovvero, prendendo i moduli,jx( k+1) 1=aj=j1=(2 k) jjx( k) 1=aj2 ; k2 [x k; 1=a]:(1) Si noti che avremmo potuto ottenere tale espressione direttamente, sapendo che la funzione di itera- zione N( x) =xg(x)associata al metodo di Newton è tale che0 N( ) = 0; da uno sviluppo di Taylor con resto in forma di Lagrange abbiamo, allora, x( k+1) = N( x( k) )() =12!  00 N(  k)( x( k) )2 ; da cui, prendendo i moduli e sostituendo l'espressione di00 N, si ottiene la (1). Sex(0) è sucientemente vicino a1=a, il metodo di Newton convergerà ; in particolare, k! 1=ae dalla (1) si avràjx( k+1) 1=aj< Mjx( k) 1=aj2 conM=12 j kj'1a . c)Sostituiamo il terminelogaxcon il suo sviluppo al second'ordine centrato in1=a, ovvero logax'a(x1=a) +12 a 2 (x1=a)2 perx!1=a: Otteniamo il metodo seguente :x( k+1) =(x( k) ); con(x) =x 1(ax1) +12 ( ax1)2 : Osserviamo che, per essere valutata, la funzione di iterazionenon richiede altre operazioni oltre a somma, sottrazione e moltiplicazione. Inoltre1=aè un punto sso die si ha 0 (1=a) = 0: Grazie al Teorema di Convergenza Locale e ai risultati sull'ordine di convergenza, abbiamo allora che il nostro metodo converge a1=a(purchéx(0) sia sucientemente vicino a1=a) e risulta, come il metodo di Newton, di ordine 2.4 Esercizio 3 Siauna radice doppia di una funzionef:R!R, ovvero tale chef() =f0 () = 0,f00 ()6 = 0. Si supponga fdi classeC2 (I)in un intornoIdi. a)A partire dalla relazionef(x) = (x)2 h(x)oveh()6 = 0, provare che il metodo di Newton per l'approssimazione dinon può avere ordine 2. b)Provare che il metodo modicato x( k+1) =x( k) 2f (x( k) )f 0 (x( k) ) converge localmente adcon ordine 2. c)Indicare un possibile test d'arresto per il metodo di Newton (non modicato), ricavando una stima dell'errore in funzione della tolleranza. d)Come si estendono i punti precedenti al caso in cuisia una radice di molteplicitàm2N? Soluzione 3a)Il metodo di Newton si denisce come segue x( k+1) = N( x( k) ) =x( k) f (x( k) )f 0 (x( k) ): Nel caso in cuiè una radice di molteplicità2, la funzionef(x)si può riscrivere come segue, f(x) =h(x)(x)2 ; h()6 = 0: Sia N( x)la funzione di iterazione associata al metodo di Newton : allora 0 N( x) = 1f 0 (x)2 f(x)f00 (x)f 0 (x)2=f (x)f00 (x)( f0 (x))2: Abbiamo poif(x) =h(x)(x)2 ; f0 (x) =h0 (x)(x)2 + 2h(x)(x); f00 (x) =h00 (x)(x)2 + 4h0 (x)(x) + 2h(x): Sostituendo nell'espressione di0 N( x)otteniamo 0 N( x) =[ h(x)(x)2 ][h00 (x)(x)2 + 4h0 (x)(x) + 2h(x)][ h0 (x)(x)2 + 2h(x)(x)]2; da cui si ricavalim x! 0 N( x) =2 h2 ()4 h2 ()= 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( )j=12 : b)Il metodo modicato ha una funzione di iterazione denita da M( x) =x2f (x)f 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.5 c)Dal teorema di Lagrange, x( k+1) =(x( k) )() =0 ( k)( x( k) ); con kcompreso fra x( k) e. Dunque si ha 0 ( k)( x( k) ) =x( k+1) x( k) +x( k) , (x( k) ) =11 0 ( k)( x( k) x( k+1) ): Pertanto, assumendo0 ( k) '0 (), si ha la seguente stima dell'erroree( k) =jx( k) jin funzione dell'incrementod( k) =jx( k+1) x( k) j: e( k) '11 0 ()d ( k) ; che per il metodo di Newton è ottimale seè radice singola, in quanto risultae( k) 'd( k) . Per una radice doppia, da quanto trovato avremoe( k) '2d( k) . In ogni caso, il test d'arresto sarà del tipo d( k) < toll, essendotolluna tolleranza sull'errore ssata. d)Nel caso generale in cuiè uno zero di molteplicitàm2dell'equazione non linearef(x) = 0, procederemo come sopra. In particolare, per una opportuna funzionehtale cheh()6 = 0, avremo : f(x) =h(x)(x)m ; f0 (x) =h0 (x)(x)m +mh(x)(x)m 1 ; f00 (x) =h00 (x)(x)m + 2mh0 (x)(x)m 1 +m(m1)h(x)(x)m 2 : Sostituendo nell'espressione di0 N( x)otteniamo 0 N( x) =[ h(x)(x)m ][h00 (x)(x)m + 2mh0 (x)(x)m 1 +m(m1)h(x)(x)m 2 ][ h0 (x)(x)m +mh(x)(x)m 1 ]2; da cui si ricavalim x! 0 N( x) =m (m1)h2 ()m 2 h2 ()= m 1m = 1 1m : In particolare, si evince che il metodo di Newton (non modicato) è localmente convergente, ma solo al primo ordine. Il metodo modicato del secondo ordine corrisponde ad un metodo iterativo di punto sso per la seguente funzione di iterazione : M( x) =xtf (x)f 0 (x); ove la costantet2Rè scelta in maniera da avere0 M( ) = 0. Dato che M( x) =xtf (x)f 0 (x)= t N( x)(t1)x; per avere0 M( ) = 0imporremo 0 = limx! 0 M( x) =t
m 1m (t1); da cuit=m.6 Esercizio 4 Per il calcolo dello zero della funzionef(x) =x3 x2 + 8x8si usano 4 diversi metodi di punto sso, descritti rispettivamente dalle seguenti funzioni di iterazione : 1( x) =x3 +x2 7x+ 8,  2( x) =8 x38 x,  3( x) =110 x3 +110 x2 +15 x +45 ,  4( x) =2 x3 x2 +83 x2 2x+8. Nella tabella seguente sono riportate, in ordine sparso, le iterate successive ottenute dai 4 metodi. Giusti- cando adeguatamente la risposta, abbinare le quattro successioni a ciascuno dei quattro metodi.Metodo A Metodo B Metodo C Metodo D 5.000000000000000e-01 5.000000000000000e-01 5.000000000000000e-01 5.000000000000000e-01 9.125000000000001e-01 1.032258064516129e+00 4.625000000000000e+00 1.050000000000000e+00 9.897857421875000e-01 1.000235245684712e+00 -1.019160156250000e+02 9.845143884892086e-01 9.989578145726552e-01 1.000000012299503e+00 1.069697123778202e+06 1.004312677086027e+00 9.998955643403695e-01 1.000000000000000e+00 -1.224001861234915e+18 9.987590594698483e-01 9.999895542527895e-01 1.000000000000000e+00 1.833775789385161e+54 1.000353832012369e+00 9.999989554034564e-01 1.000000000000000e+00 -6.166499545700052e+162 9.998988463640411e-01Soluzione 4 Osserviamo prima di tutto chef(x) = (x1)(x2 + 8), quindi l'unico zero reale dif(x)è= 1. Si noti che, grazie al teorema del valor medio, si ha : x( k+1) =0 ( k)( x( k) ) ove, almeno perksucientemente grande,0 ( k) '0 (). In particolare, sej0 ()j0le iteratex( k) si avvicineranno adin maniera monotona, mentre per0 ()