Metodo del gradiente coniugato

Da Wikipedia, l'enciclopedia libera.

Questa voce di matematica è solo un abbozzo: contribuisci a migliorarla secondo le convenzioni di Wikipedia.

In analisi numerica, il metodo del gradiente coniugato è un algoritmo per la risoluzione di un sistema lineare della forma $A x = b$ .

Per poter applicare questo metodo è necessario che A sia una matrice simmetrica definita positiva cioè $\forall x, \, x^t A x \geq 0$ .

[modifica] Spiegazione del metodo

Quando $A$ è simmetrica si ha che la minimizzazione della funzione $f(x) = {1 \over 2} x^t A x -b^t \cdot x + c$ è equivalente alla risoluzione del sistema iniziale $A x = b$ , ed avviene quando
$\bigtriangledown f(x) = A x - b = 0$
cioè quando:
$\bigtriangledown f(x) = 0 \Leftrightarrow A x - b = 0 \Leftrightarrow A x = b$

Questo significa che risolvere il sistema lineare è equivalente a trovare lo zero di $\bigtriangledown f$ e quindi il minimo di $f$ .

Per trovare il minimo si procede iterativamente come nel metodo del gradiente cioè, partendo da un punto $x i$ si sceglie una direzione $d i$ e si trova il punto successivo minimizzando $f$ sulla retta avente direzione $d i$ e passante per $x i$ .

La differenza sostanziale rispetto al metodo del gradiente è nella scelta della direzione $d i$ che non corrisponde alla direzione di massima pendenza come nel metodo del gradiente, ma è scelta $A$ -ortogonale alle superfici di livello.

La ragione di questa scelta è che $f (x)$ è un paraboloide dilatato con coefficiente $c i$ sulla direzione $e i$ e poi traslato. I $c i$ sono gli autovalori di $A$ e gli $e i$ i corrispondenti autovettori che esistono per il teorema spettrale. La direzione scelta coincide con il gradiente nel paraboloide originale. Questo assicura la convergenza in $n - 1$ passi.

Il metodo del gradiente biconiugato fornisce una generalizzazione per matrici non simmetriche.

[modifica] Implementazione d'esempio nel linguaggio Octave

 function [x] = conjgrad(A,b,x0)
  
 r = b - A*x0;
 w = -r;
 z = A*w;
 a = (r'*w)/(w'*z);
 x = x0 + a*w;
 B = 0;
  
 for i = 1:size(A)(1);
    r = r - a*z;
    if( r < 1e-10)
         break;
    endif
    B = (r'*z)/(w'*z);
    w = -r + B*w;
    z = A*w;
    a = (r'*w)/(w'*z);
    x = x + a*w;
 end