Mathematical Engineering - Matematica Numerica

Exercise - 09 - Solution

Divided by topic

MATEMATICA NUMERICA A.A. 2018 - 2019 Ingegneria Matematica Prof. A. Quarteroni Prof. A. Manzoni, Dr. I. Fumagalli Esercitazione 9 - Soluzione Metodi iterativi per sistemi lineari (2) Esercizio 1 SiaA2Rn n una matrice simmetrica e denita positiva e si consideri il metodo del gradiente coniugato per risolvere il sistema lineareAx=b: datox(0) , postor(0) =bAx(0) ep(0) =r(0) , si calcola x( k+1) =x( k) + kp( k) ; r( k+1) =r( k)  kA p( k) ; p( k+1) =r( k+1)  kp( k) ;con  k=( p( k) ;r( k) )( p( k) ;Ap( k) ); k=( p( k) ;Ar( k+1) )( p( k) ;Ap( k) ): a)Procedendo per induzione sull'interok0, si dimostrino le relazioni seguenti : (r( k+1) ;p( i) ) = 0i= 0; : : : ; k; (p( k+1) ;Ap( i) ) = 0i= 0; : : : ; k: Occorrerà osservare che per ognim0, si ha r( m) 2V m= spanfr(0) ;Ar(0) ; : : : ;Am r(0) g;p( m) 2V m: b)Sfruttando quanto trovato al punto precedente, dimostrare che (in aritmetica esatta) il metodo delgradiente coniugato termina al più doponiterazioni. c)SianoP 1 ;P 1=2 due matrici simmetriche e denite positive tali cheP 1=2 P 1=2 = P 1 , e si consideri il seguentesistema precondizionato P 1=2 AP 1=2 |{z} b AP 1 =2 x |{z} b x= P 1=2 b |{z} b b Applicando il metodo del gradiente coniugato al sistema precondizionatob Ab x=b b, si derivi il metodo del gradiente coniugato precondizionato per il sistemaAx=b. In particolare, si evidenzino il ruolo del residuo precondizionatoz( k) = P 1 r( k) e della direzione di discesap( k) = P 1=2 b p( k) , doveb p( k) denota la direzione di discesa del gradiente coniugato per il sistemab Ab x=b b. Soluzione 1a)Per dimostrare la proprietà richiesta, procediamo per induzione. Supponiamo che, per un arbitrariointerok1, sia (r( k) ;p( i) ) = 0i= 0; : : : ; k1; (p( k) ;Ap( i) ) = 0i= 0; : : : ; k1: Vogliamo dimostrare che allora vale (r( k+1) ;p( i) ) = 0i= 0; : : : ; k; (p( k+1) ;Ap( i) ) = 0i= 0; : : : ; k:(1)1 Si osservi che, per denizione dei coecienti  ke  k, vale in ogni caso (r( k+1) ;p( k) ) = 0;(p( k+1) ;Ap( k) ) = 0; il che incidentalmente prova che l'ipotesi di induzione è sicuramente vericata perk= 1e che le (1) sono vere peri=k. Pertanto, rimane da provare che (r( k+1) ;p( i) ) = 0i= 0; : : : ; k1; (p( k+1) ;Ap( i) ) = 0i= 0; : : : ; k1:(2) Cominciamo dal residuo : peri < ksi ha (r( k+1) ;p( i) ) = (r( k) ;p( i) ) |{z} 0  k(A p( k) ;p( i) ) |{z} 0= 0 ;(3) ove abbiamo usato l'ipotesi di induzione. Rimane da considerare(p( k+1) ;Ap( i) ) = (r( k+1) ;Ap( i) ) + k( p( k) ;Ap( i) ) = (r( k+1) ;Ap( i) ); sempre grazie all'ipotesi di induzione. Il passo nale è dunque provare che(r( k+1) ;Ap( i) ) = 0sei < k. Per fare ciò, si osservi che per ognim >0, si ha r( m) 2V m= spanfr(0) ;Ar(0) ; : : : ;Am r(0) g;p( m) 2V m: Lo spazioV mè anche noto come spazio di KrylovK m+1(A ;r(0) )di ordinem+ 1associato alla matrice Ae al residuo iniziale. Ora, peri < k, dap( i) 2V k1otteniamo Ap( i) 2V k: (4) Osserviamo che essendoAsimmetrica e denita positiva,(x;y) A= ( x;Ay)è un prodotto scalare inRn . Rispetto a tale prodotto scalare, tutte le direzioni di discesafp(0) ; : : : ;p( k) gsono ortogonali, dunque linearmente indipendenti. In particolare, esse formano una base perV k. Dalla (3) si desume cher( k+1) ?V k. Ma allora, dalla (4) abbiamo (r( k+1) ;Ap( i) ) = 0; i= 0; : : : ; k1; il che conclude la dimostrazione. b)Abbiamo già osservato come tutte le direzioni di discesap( i) ,i= 0; : : : ; n1sono linearmente indipendenti perché A-ortogonali. Pertanto, se allan-esima iterazione nessuna direzione di discesa si è annullata, lep( i) formano una base (A-ortogonale) diRn . Ma allora, dal punto a) abbiamo che il residuon-esimo è necessariamente nullo, in quanto esso è A-ortogonale a tutti i vettori di tale base. Pertanto il metodo termina allan-esima iterazione restituendo la soluzione esatta (in aritmetica a precisione innita). Nel caso invece in cui risultip( k) =0per un opportunok < n1, avremmo p( k) =0=)r( k) = k1p( k1) : Tuttavia sappiamo cher( k) è ortogonale ap( k1) , dunque 0 = (r( k) ;p( k1) ) = k1( p( k1) ;p( k1) ) =) k1p( k1) = 0 =)r( k) =0; ovverox( k) =xe anche in tal caso l'algoritmo termina restituendo la soluzione esatta. Questo conclude la prova. Vediamo ora di approfondire il confronto fra il metodo del gradiente ed il metodo del gradiente coniu- gato. Con riferimento al funzionale energiatale cher( k+1) =r(x( k+1) ), abbiamo che la prima delle (1) ci dice che il gradiente diinx( k+1) è ortogonale a tutte le direzionip(0) ; : : : ;p( k) . Quindi,2 la soluzione approssimata x( k+1) minimizzasul sottospazio generato da queste direzioni (che altro non è se nonK k+1(A ;r(0) )) : (x( k+1) ) = min 2R;i=0;:::;k( x( k) +p( i) ): Si noti che in questa relazione emerge la dierenza fra il metodo del gradiente e quello del gradiente coniugato : per il primo, la relazione in questione vale restringendo il minimo su2R; i=k, ovvero l'ottimalità rispetto alle precedenti direzioni di discesa non è mantenuta. Nel metodo del gradiente coniugato, invece, la soluzione approssimata al passok+ 1minimizza il funzionale energia sul sottos- pazio di tutte le direzioni di discesa precedentemente calcolate. Come si è già discusso, questo implica che il metodo del gradiente coniugato termina in al piùnpassi (in aritmetica a precisione innita). c)Osserviamo innanzitutto che, essendoA;P 1=2 simmetriche e denite positive, ancheb Alo è, dunque possiamo eettivamente applicare il metodo del gradiente coniugato al sistema precondizionato. In maniera analoga ab A;b b;b x, adottiamo la notazioneb
per identicare le quantità relative al sistema precondizionato. Inoltre, osserviamo come br( k) =b bb Ab x( k) = P 1=2 bP 1=2 AP 1=2 P1 =2 x( k) = P 1=2 (bAx( k) ) = P 1=2 r( k) L'inizializzazione del metodo si può dunque scrivere come segue :b x(0) = P1 =2 x(0) ;b r(0) =b b(0) b Ab x(0) = P 1=2 r(0) ;b p(0) =b r(0) = P 1=2 r(0) : Le direzioni di discesap( k) del metodo del gradiente coniugato precondizionato sono legate a quelle del metodo del gradiente coniugato applicato al sistemab Ab x=b b, indicate conb p( k) , dalla relazione seguente, P1 =2 p( k) =b p( k) risulta infatti logico aspettarsi, in base alle formule di aggiornamento delle iteratex( k) , che lo stesso cambio di variabili che adottiamo per queste ultime, P1 =2 x( k) =b x( k) valga anche per le direzioni di discesa. Risulta quindi, dal momento cheb p(0) =b r(0) P1 =2 p(0) =b p(0) = P 1=2 r(0) )p(0) = P 1 r(0) : Per ogni passo successivo del metodo, calcoliamo i coecienti k;  k, sostituendo i valori delle variabili b
e usando le denizioniz( k) = P 1 r( k) ;p( k) = P 1=2 b p( k) : k= b p( k) ;b r( k)
 b p( k) ;b Ab p( k) = P1 =2 p( k)
T P 1=2 r( k) P1 =2 p( k)
T P 1=2 AP 1=2 P1 =2 p( k)= p( k) ;r( k)
p( k) ;Ap( k)
; k= b p( k) ;b Ab r( k+1) b p( k) ;b Ab p( k) = P1 =2 p( k)
T P 1=2 AP 1=2 P 1=2 r( k+1) P1 =2 p( k)
T P 1=2 AP 1=2 P1 =2 p( k)= p( k) ;Az( k+1)
p( k) ;Ap( k)
: Risulta inoltrebr( k+1) =b r( k)  kb Ab p( k) =P 1=2 r( k)  kP 1=2 AP 1=2 P1 =2 p( k) =P 1=2 (r( k)  kA p( k) ) da cui ricaviamo cheP 1=2 r( k+1) = P 1=2 (r( k)  kA p( k) ))r( k+1) =r( k)  kA p( k) ; analogamente,p( k+1) = P 1=2 b p( k+1) = P 1=2 (b r( k+1)  kb p( k) ) = P 1=2 (P 1=2 r( k+1)  kP1 =2 p( k) ) = P 1 r( k+1)  kp( k) =z( k+1)  kp( k)3 da cui p( k+1) =z( k+1)  kp( k) : L'algoritmo del metodo del gradiente coniugato precondizionato risulta dunque come segue : datox(0) , sianor(0) =bAx(0) ;z(0) = P 1 r(0) ;p(0) =z(0) ; si calcolax( k+1) =x( k) + kp( k) ; r( k+1) =r( k)  kA p( k) ; Pz( k+1) =r( k+1) ; p( k+1) =z( k+1)  kp( k) ;con  k=( p( k) ;r( k) )( p( k) ;Ap( k) ); k=( p( k) ;Az( k+1) )( p( k) ;Ap( k) ):4 Esercizio 2 Costruire la matrice A1= 6:82:4 2:4 8:2 e vericare che ha come autovalori 1= 5 , 2= 10 e come autovettori v1=15  4 3 ;v 2=15  3 4 : Costruire il vettoreb= 4 8 e visualizzare la forma quadratica' 1( x) =12 xT A1x
xT bassociata e le corrispondenti linee di livello. Cosa succede alla forma quadratica se si cambia il primo autovalore in 1= 1 ? In particolare, si deniscaA 2= V DVT ;doveV= [v 1; v 2] eD=diag(1;  2) e si consideri la forma quadratica 2( x) =12 xT A2x
xT b. a)Utilizzare il metodo di Richardson stazionario scegliendo= 0:05oppure= 0:24e il metodo del gradiente (o metodo di Richardson dinamico con parametro ottimale) per risolvere il problema di minimizzazione della forma quadratica 2. In entrambi i casi si usi la funzione richardsonfornita, considerando la versione non precondizionata del metodo (P=I). Visualizzare la successione dei punti che si ottiene in ciascuno dei tre casi. (Suggerimento : modicare la funzionerichardson.m creando una nuova funzionemy_richardsonche implementi lo stesso metodo, ma restituisca in uscita la matrice che ha per colonne la successione delle iteratex( k) , e ne riporti il loro graco.) b)Considerare ora il metodo del gradiente precondizionato, indicando conPla matrice di precondizio- namento ; essa deve essere simmetrica e denita positiva, e tale cheP 1 Asia simmetrica e denita positiva. Risolvere con il metodo del gradiente precondizionato il sistema lineare associato alla forma quadratica 2, utilizzando come precondizionatore la matrice P= 1:09120:8587 0:8587 1:5921 ; visualizzare le iterate e confrontarle con quelle calcolate dal metodo del gradiente : cosa si osserva ? c)(*) Risolvere il sistemaP 1 Ax=P 1 bcon il metodo del gradiente è equivalente a risolvere il sistema Ax=bcon il metodo del gradiente precondizionato ? Si confrontino le iterate ottenute al punto precedente con quelle calcolate dal metodo del gradiente nel caso del sistemaP 1 A2x =P 1 b, a cui è possibile associare una nuova forma quadratica prec( x) =12 xT P 1 A2x
xT P 1 b. d)Usando il template fornito, implementare una funzioneconjgrad.mche realizzi il metodo del gradiente coniugato precondizionato, usando un criterio d'arresto basato sul residuo normalizzato. L'intestazione di tale funzione è[x, niter, error, flag] = conjgrad(A, b, x0, maxit, tol, P). e)Analogamente al punto 1, scrivere una nuova funzionemy_conjgradche restituisca restituisca la ma- trice che ha per colonne la successione delle iteratex( k) , e ne riporti il loro graco. Risolvere il sistema lineare associato alla forma quadratica 2con il metodo del gradiente coniugato non precondizionato e visualizzare il percorso seguito nella discesa verso il punto di minimo. Soluzione 2 Per visualizzare le forme quadratiche inR3 e le linee di livello nel piano(x; y)si usano rispettivamente i comandisurfecontour; i graci ottenuti sono riportati in Figura 1. Le linee di livello di entrambe le forme quadratiche sono delle ellissi, e quelle della forma quadratica associata alla matriceA 2sono più eccentriche, perché gli autovalori diA 2sono molto diversi fra loro.v1=([4 3]/5); v2=([ 3 4]/5);l1=5; l2=10; V = [v1,v2]; 5 D1 = diag([l1,l2]); A1 = V *D1 *V;l1=1; l2=10; D2 = diag([l1,l2]); A2 = V *D2 *V;b=[4 8]; x=linspace( 10,10,80); y=linspace(10,10,80);[X,Y]=meshgrid(x,y); Phi1 = 0.5 *( A1(1,1) *X.^2 + A1(2,2) *Y.^2 + 2 *A1(1,2) *X. *Y ) b(1) *X b(2) *Y;Phi2 = 0.5 *( A2(1,1) *X.^2 + A2(2,2) *Y.^2 + 2 *A2(1,2) *X. *Y ) b(1) *X b(2) *Y;surf(X,Y,Phi1,Lines,no); title(phi1); hold on contour(X,Y,Phi1,20) figure; surf(X,Y,Phi2,Lines,no); title(phi2); hold on contour(X,Y,Phi2,20) figure contour(X,Y,Phi1,20); title(phi1 lineedilivello)figure contour(X,Y,Phi2,20); title(phi2 lineedilivello)Figure 1  Forme quadratiche (sopra) e linee di livello (sotto) : a sinistra 1, a destra  2 a)La funzionemy_richardsonè la seguente :function [xk,k] = my_richardson(A,b,P,x0,toll,nmax,alpha)%MetododiRichardsonstazionarioodinamicoconparametroottimale(=gradiente) % %Parametridiingresso: 6−10 −5 0 5 10 −10 −5 0 5 10 −200 0 200 400 600 800 1000 1200 −10 −5 0 5 10 −10 −5 0 5 10 −200 0 200 400 600 800 1000 1200 −10 −8 −6 −4 −2 0 2 4 6 8 10 −10 −8 −6 −4 −2 0 2 4 6 8 10 −10 −8 −6 −4 −2 0 2 4 6 8 10 −10 −8 −6 −4 −2 0 2 4 6 8 10 % %A:matricedisistema %b:vettoreterminenoto %P:precondizionatoredelmetodo.SiusaP=Inelcasononprecondizionato. %x0:guessiniziale %toll:tolleranzacriteriodarresto %nmax:numeromassimodiiterazioniammesse %alpha:parametrodiaccelerazione;senonassegnatosiconsiderailmetododinamico % %Parametridiuscita: % %xk:matricelecuicolonnesonoivettoridellasoluzioneadogniiterazione %k:numerodiiterazioni n = length(b); if ((size(A,1)6 =n) || (size(A,2)6 =n) || (length(x0)6 =n))error(Dimensioniincompatibili) end x = x0; k = 0; r = b A *x;res = norm(r) / norm(b); xk = [x0]; while ((res > toll) && (k < nmax))z = P \ r; if (nargin == 6)alpha = (z *r) / (z *A *z); %alphadinamicoend x = x + alpha *z;r = b A *x;res = norm(r) / norm(b); k = k + 1; xk=[xk x]; end if (res < toll)fprintf(Richardsonconvergein%diterazioni\n, k); else fprintf(Richardsonnonconvergein%diterazioni.\n, k) end end Utilizzando come tolleranza 10 7 , il metodo di Richardson stazionario con passo= 0:05converge in331iterazioni, mentre con= 0:24diverge (si vedano i comandi riportati sotto). Dalla Figura 2 è evidente che la scelta di utilizzare un passo costante per risolvere questo tipo di sistemi può essere estremamente penalizzante anche in caso di convergenza.x0=[ 99]; nmax=1000; toll=1e7;[Sg_rich,n_it_rich] = my_richardson(A2,b,eye(2),x0,toll,nmax,0.05); figure; contour(X,Y,Phi2,80) title(\alpha=0.05); axis equal; hold on; plot(Sg_rich(1,:),Sg_rich(2,:), or,LineWidth,2,MarkerSize,4,MarkerFaceColor,r)x0=[ 99]; nmax=1000; toll=1e7;[Sg_rich2,n_it_rich2] = my_richardson(A2,b,eye(2),x0,toll,nmax,0.24); figure; contour(X,Y,Phi2,80) title(\alpha=0.24); axis equal; hold on; plot(Sg_rich2(1,1:6),Sg_rich2(2,1:6), or,LineWidth,2)Il metodo del gradiente, che adatta il passo ad ogni iterazione, è molto più eciente e converge in 877 Figure 2  Iterazioni di Richardson : a sinistra il caso convergente, a destra il caso non convergente iterazioni.x0=[ 99]; nmax=100; toll=1e7;[Sg2,n_it2] = my_richardson(A2,b,eye(2),x0,toll,nmax); figure; contour(X,Y,Phi2,80) title(Gradiente); axis equal; hold on; plot(Sg2(1,:),Sg2(2,:), or,LineWidth,2,MarkerSize,4,MarkerFaceColor,r)Tuttavia anche la discesa eettuata da questo metodo non è ottimale, come si può vedere in Figura 3, dove è evidente la proprietà del metodo del gradiente di muoversi lungo direzioni perpendicolari.Figure 3  Iterazioni del metodo del gradiente b)Nel caso del metodo del gradiente precondizionato, usiamo i seguenti comandi :x0=[ 99]; nmax=100; toll=1e7;P=V *diag([l1,5]) *V;[Sg_prec,n_it_prec] = my_richardson(A2,b,P,x0,toll,nmax); figure; contour(X,Y,Phi2,60) axis equal hold on; title(GradientePrecondiz.vsGradiente) %disegnoingrigioledirezionidelgradientenormale 8−10 −5 0 5 10 −10 −8 −6 −4 −2 0 2 4 6 8 10 −10 −5 0 5 10 −10 −8 −6 −4 −2 0 2 4 6 8 10 −10 −5 0 5 10 −10 −8 −6 −4 −2 0 2 4 6 8 10 plot(Sg2(1,:),Sg2(2,:), ,LineWidth,2,Color,[0.48 0.48 0.48])%poisovrappongoinrossoquelledelgradienteprecondizionato plot(Sg_prec(1,:),Sg_prec(2,:), or,LineWidth,2,MarkerSize,4,MarkerFaceColor,r)c)La forma quadratica associata al sistema precondizionato è  prec( x) =12 xT P 1 Ax
xT P 1 b. Con la particolare scelta diPindicata, gli autovettori della matriceP 1 A2sono gli stessi della matrice A2, ma gli autovalori sono  1= 1 ;  2= 5 , e quindi le linee di livello sono simili a delle circonferenze (vedi Figura 4).Aprec=P\A2; eig(Aprec) bprec=P\b; Phiprec = 0.5 *(Aprec(1,1) *X.^2+Aprec(2,2) *Y.^2+2 *Aprec(1,2) *X. *Y) bprec(1) *X bprec(2) *Y;figure contour(X,Y,Phi2,80); grid on title(phi2 lineedilivello)figure contour(X,Y,Phiprec,80); grid on title(phiprec lineedilivello)Figure 4  Linee di livello : a sinistra quelle di 2, a destra quelle di  prec Applicare il metodo del gradiente precondizionato al sistema lineare originalenonequivale ad applicare il metodo del gradiente al sistema lineare precondizionatoP 1 A2x =P 1 b, come si vede in Figura 5. Tuttavia le iterazioni richieste sono simili :14nel primo caso e9nel secondo.[Sg_phiprec,n_it_phiprec] = my_richardson(Aprec,bprec,eye(2),x0,toll,nmax); figure; contour(X,Y,Phiprec,60); axis equal title(GradientesuP^{ 1}A); hold on;plot(Sg_phiprec(1,:),Sg_phiprec(2,:), or,LineWidth,2,MarkerSize,4,MarkerFaceColor,r)Si può dimostrare che le direzioni utilizzate dal metodo del gradiente precondizionato sono tali che zT k+1P z k= 0 (direzioniP-ortogonali). Questa proprietà in generale non è valida per due direzioni non consecutive (cioè, ad esempio,zT 3P z 16 = 0), e quindi il metodo del gradiente precondizionato non è ancora ottimale. Numericamente, possiamo vericare tale proprietà nel seguente modo (se si ottiene 0 le direzioni sono ortogonali, altrimenti non lo sono).%direzionidelgradienteprecondizionato dir_prec=Sg_prec(:,2:5) Sg_prec(:,1:4)9−10 −8 −6 −4 −2 0 2 4 6 8 10 −10 −8 −6 −4 −2 0 2 4 6 8 10 −10 −8 −6 −4 −2 0 2 4 6 8 10 −10 −8 −6 −4 −2 0 2 4 6 8 10 Figure 5  La gura di sinistra riporta le iterazioni del gradiente precondizionato in rosso e quelle del gradiente in grigio ; quella di destra le iterazioni del gradiente applicato al sistema precondizionato.%direzionidelgradienteapplicatoalsistemaprecondizionato dir_phiprec=Sg_phiprec(:,2:5) Sg_phiprec(:,1:4)%ledirezionidelgradienteapplicatoalsistemaprecondizionatodevonoessereortogonali disp(gradienteapplicatoalsistemaprecondizionato:testortogonalita) dir_phiprec(:,1) *dir_phiprec(:,2)dir_phiprec(:,2) *dir_phiprec(:,3)dir_phiprec(:,3) *dir_phiprec(:,4)%mentrequelledelgradienteprecondizionatoingeneralenonlosaranno... disp(gradienteprecondizionato:testortogonalita) dir_prec(:,1) *dir_prec(:,2)dir_prec(:,2) *dir_prec(:,3)dir_prec(:,3) *dir_prec(:,4)%...masonoadueaduePortogonali: disp(gradienteprecondizionato:testP ortogonalita)dir_prec(:,1) *P *dir_prec(:,2)dir_prec(:,2) *P *dir_prec(:,3)dir_prec(:,3) *P *dir_prec(:,4)%direzioninonconsecutiveinveceingeneralenonsarannoPortogonali disp(gradienteprecondizionato:testP ortogonalita)dir_prec(:,1) *P *dir_prec(:,3)dir_prec(:,1) *P *dir_prec(:,4)dir_prec(:,2) *P *dir_prec(:,4)Poiché z( k+1) =P 1 r( k+1) =P 1 r( k)  kP 1 Az( k) =z( k)  kP 1 Az( k) risulta infatti(Pz( k+1) ;z( k) ) = (Pz( k)  kA z( k) ;z( k) ) = (Pz( k) ;z( k) )(z( k) ;r( k) ) = (r( k) ;z( k) )(z( k) ;r( k) ) = 0:10−10 −5 0 5 10 −10 −8 −6 −4 −2 0 2 4 6 8 10 −10 −5 0 5 10 −10 −8 −6 −4 −2 0 2 4 6 8 10 d)La funzione conjgrad.mrisulta la seguente :function [x, niter, error, flag] = conjgrad(A, b, x0, maxit, tol, P)% %[x,iter,error,flag]=conjgrad(A,b,x0,maxit,tol,P) % %ConjugategradientmethodforthesolutionofP^{ 1}Ax=P^{1}b.% %>>INPUT>OUTPUTtol flag = 0; niter = 0; bnrm2 = norm( b ); x=x0; r = b A *x; error = norm( r ) / bnrm2;if ( error < tol )return,endfor niter = 1:maxitz = P \ r; rho = (r *z);if niter > 1beta = rho / rho1; p = z + beta *p;else p = z; end q = A *p; alpha = rho / (p *q );x = x + alpha *p; r = r alpha *q;error = norm( r ) / bnrm2; if ( errortol ),break,endrho1 = rho; end if ( error > tol ) flag = 1;endreturn e)Applicare il metodo del gradiente coniugato consente di mantenere l'ottimalità fra tutte le direzioni, e di convergere alla soluzione esatta in2iterazioni. Le direzioni selezionate dal metodo del gradiente coniugato sono riportate in Figura 6, e sono ottenute con i seguenti comandi :Figure 6  La gura mostra in grigio le iterazioni del gradiente precondizionato e in rosso quelle del gradiente coniugato11−10 −5 0 5 10 −10 −8 −6 −4 −2 0 2 4 6 8 10 [Sg_conj,n_it_conj, error, flag] = my_conjgrad(A2, b, x0, nmax, toll, eye(2)); figure; contour(X,Y,Phi2,60); axis equal; hold on; title(GradienteConiugatovsPrecondizionato) plot(Sg_prec(1,:),Sg_prec(2,:), ,LineWidth,2,Color,[0.48 0.48 0.48])plot(Sg_conj(1,:),Sg_conj(2,:), or,LineWidth,2)La funzione my_conjgrad.mche memorizza la successione delle iterate nel caso del metodo del gradiente coniugato risulta la seguente :function [xk, niter, error, flag] = my_conjgrad(A, b, x0, maxit, tol, P)flag = 0; niter = 0; bnrm2 = norm( b ); x=x0; xk = [x0]; r = b A *x; error = norm( r ) / bnrm2;if ( error < tol )return,endfor niter = 1:maxitz = P \ r; rho = (r *z);if niter > 1beta = rho / rho1; p = z + beta *p;else p = z; end q = A *p; alpha = rho / (p *q );x = x + alpha *p; r = r alpha *q;error = norm( r ) / bnrm2; xk = [xk x]; if ( errortol ),break,endrho1 = rho; end if ( error > tol ) flag = 1;endreturn 12