– Typeset by GMNI & FoilTEX –
METODOS SEMI-ITERATIVOS PARAGRANDES SISTEMAS DE ECUACIONES LINEALES:
DIRECCIONES CONJUGADAS. PCGF. Navarrina, I. Colominas, M. Casteleiro, H. Gomez, J. Parıs
GMNI — GRUPO DE METODOS NUMERICOS EN INGENIERIA
Departamento de Metodos Matematicos y de RepresentacionEscuela Tecnica Superior de Ingenieros de Caminos, Canales y Puertos
Universidad de A Coruna, Espana
e-mail: [email protected] web: http://caminos.udc.es/gmni
UNIVERSIDAD DE A CORUNA — GRUPO DE METODOS NUMERICOS EN INGENIERIA
INDICE
I Direcciones Conjugadas
I Metodos de Direcciones Conjugadas
I Metodo de Gradientes Conjugados [CG]
I Planteamiento de Mınima Energıa
I Precondicionamiento
I Gradiente Conjugado Precondicionado [PCG]
I Implementacion del Algoritmo PCG
UNIVERSIDAD DE A CORUNA — GRUPO DE METODOS NUMERICOS EN INGENIERIA
Direcciones Conjugadas (Ia)
Sea el sistema
A˜ x = b, con
{A˜ = A˜T (A˜ SIMETRICA),
vTA˜ v > 0 ∀v 6= 0 (A˜ DEF. +).
Sean los vectores conjugados respecto a la matriz A˜ (o A˜–conjugados),
{si}i=1,n, tales que sTi A˜ sj{
= 0 si i 6= j,
6= 0 si i = j,
=⇒ FORMAN UNA BASE.
UNIVERSIDAD DE A CORUNA — GRUPO DE METODOS NUMERICOS EN INGENIERIA
Direcciones Conjugadas (Ib)
Pues. . .
REDUCCION AL ABSURDO
Hipotesis:n∑i=1
λisi = 0 con algun λj 6= 0,
entonces sTj A˜(
n∑i=1
λisi
)=
n∑i=1
λi(sTj A˜ si
)= λj s
Tj A˜ sj︸ ︷︷ ︸>0
= 0,
luego λj = 0 ∀j (ABSURDO).
Por tanto, los vectores conjugados son linealmente independientes y(puesto que su numero n iguala al orden de la matriz) forman una basedel correspondiente espacio vectorial.
UNIVERSIDAD DE A CORUNA — GRUPO DE METODOS NUMERICOS EN INGENIERIA
Direcciones Conjugadas (IIa)
EXPRESION DE LA SOLUCIONSean
x = A˜−1b la solucion del problema, y
x0 una aproximacion a x.
El vector (x− x0) se puede escribir como combinacion lineal de loselementos de la base de vectores conjugados.Luego,
x− x0 =
n∑i=1
αisi =⇒ x = x0 +
n∑i=1
αisi.
UNIVERSIDAD DE A CORUNA — GRUPO DE METODOS NUMERICOS EN INGENIERIA
Direcciones Conjugadas (IIb)
Pero,A˜ x = b =⇒ A˜
(x0 +
n∑i=1
αisi
)= b
=⇒n∑i=1
αiA˜ si =(b− A˜ x0
)y premultiplicando por sTj obtenemos
sTj
(n∑i=1
αiA˜ si)
= sTj
(b− A˜ x0
)=⇒
n∑i=1
αi(sTj A˜ si
)= s
Tj
(b− A˜ x0
),
=⇒ αj(sTj A˜ sj
)= s
Tj
(b− A˜ x0
).
Por tanto
αj =sTj(b−A˜ x0
)sTj A˜ sj , j = 1, . . . , n.
UNIVERSIDAD DE A CORUNA — GRUPO DE METODOS NUMERICOS EN INGENIERIA
Metodos de Direcciones Conjugadas (I)
FORMULACION SEMI–ITERATIVA:
Dados x0 y s1 −→ x1 = x0 + α1s1, con α1 =sT1(b− A˜ x0
)sT1A˜ s1
,
x1 y s2 −→ x2 = x1 + α2s2, con α2 =sT2(b− A˜ x0
)sT2A˜ s2
,
...
xn−1 y sn −→ xn = xn−1 + αnsn, con αn =sTn(b− A˜ x0
)sTnA˜ sn .
Finalmente
xn =
(. . .((x0 + α1s1) + α2s2
)+ . . .
)+ αnsn
= x0 +
n∑i=1
αisi = x = A˜−1b.
UNIVERSIDAD DE A CORUNA — GRUPO DE METODOS NUMERICOS EN INGENIERIA
Metodos de Direcciones Conjugadas (IIa)
EN GENERAL:
Dados x0, {si}i=1,n
xk+1 = xk + αk+1sk+1, con αk+1 =sTk+1
(b− A˜ x0
)sTk+1A˜ sk+1
, k = 0, . . . , n− 1
y xn = x = A˜−1b (SALVO ERRORES DE REDONDEO).
Pero tambien,
αk+1 =sTk+1
(b− A˜ x0
)sTk+1A˜ sk+1
=sTk+1
(b− A˜ xk)
sTk+1A˜ sk+1
,
pues
xk = x0 + α1s1 + . . .+ αksk =⇒ sTk+1A˜ xk = s
Tk+1A˜ x0+
α1 sTk+1A˜ s1︸ ︷︷ ︸
0
+ . . .+ αk sTk+1A˜ sk︸ ︷︷ ︸
0
.
UNIVERSIDAD DE A CORUNA — GRUPO DE METODOS NUMERICOS EN INGENIERIA
Metodos de Direcciones Conjugadas (IIb)
Y FINALMENTE:
Dados x0, {si}i=1,n
rk = b−A˜ xk, αk+1 =sTk+1rk
sTk+1A˜ sk+1,
xk+1 = xk + αk+1sk+1,
k = 0, . . . , n− 1,
y xn = x = A˜−1b (SALVO ERRORES DE REDONDEO).
PROBLEMA:
¿Como se genera la base de vectores conjugados {si}i=1,n?
UNIVERSIDAD DE A CORUNA — GRUPO DE METODOS NUMERICOS EN INGENIERIA
Metodos de Direcciones Conjugadas (IIIa)
OBTENCION DE UNA BASE DE VECTORES CONJUGADOS {si}i=1,n
I METODO DE GRAM-SCHMIT
I METODO DE GRADIENTES CONJUGADOS(caso particular de Gram-Schmit)
UNIVERSIDAD DE A CORUNA — GRUPO DE METODOS NUMERICOS EN INGENIERIA
Metodos de Direcciones Conjugadas (IIIb)
Metodo de Gram-Schmit
Se eligen los vectores {vi}i=1,n linealmente independientes.
Se obtienen los vectores conjugados {vi}i=1,n en la forma:
s1 = v1
s2 = v2 + β21s1
s3 = v3 + β31s1 + β
32s2
. . .sk = vk + β
k1 s1 + β
k2 s2 + . . .+ β
kk−1sk−1
. . .sn = vn + β
n1 s1 + β
n2 s2 + . . . + β
nn−1sn−1
con las condiciones (que permiten calcular los coeficientes{βki}i=1,k−1
)
sTj A˜ sk = 0, j = 1, . . . , k − 1.
UNIVERSIDAD DE A CORUNA — GRUPO DE METODOS NUMERICOS EN INGENIERIA
Metodos de Direcciones Conjugadas (IIIc)
Luego
sk = vk +
k−1∑i=1
βki si
sTj A˜ sk = 0,
=⇒ sTj A˜(vk +
k−1∑i=1
βki si
)= 0, j = 1, . . . , k−1,
y por consiguiente
sTj A˜ vk +
βkj
(sTj A˜ sj)︷ ︸︸ ︷
k−1∑i=1
βki
(sTj A˜ si
)= 0 =⇒ β
kj = −
sTj A˜ vksTj A˜ sj j = 1, . . . , k − 1.
UNIVERSIDAD DE A CORUNA — GRUPO DE METODOS NUMERICOS EN INGENIERIA
Metodos de Direcciones Conjugadas (IIId)
Y FINALMENTE:
Dados {vi}i=1,n
s1 = v1,
βkj = −sTj A˜ vksTj A˜ sj , j = 1, . . . , k − 1,
sk = vk +
k−1∑i=1
βki si,
k = 2, . . . , n.
PROBLEMA: Gran coste computacional, T (n3).
UNIVERSIDAD DE A CORUNA — GRUPO DE METODOS NUMERICOS EN INGENIERIA
Metodos de Direcciones Conjugadas (IV)
PREGUNTA:
¿Es posible elegir (habilmente) los vectores {vi}i=1,n de forma que lamayor parte de los coeficientes
{βki}i=1,k−1; k=2,n
sean nulos?
Respuesta:
SI: Metodo de Gradientes Conjugados
UNIVERSIDAD DE A CORUNA — GRUPO DE METODOS NUMERICOS EN INGENIERIA
Metodo de Gradientes Conjugados [CG] (Ia)
Dado el sistema
A˜ x = b,
definimos la funcion cuadratica
f(x) =1
2xTA˜ x− bT x+ c,
cuyo gradiente es
∇f =
(df
dx
)T= A˜ x− b.
Luego,
∇f(xi) = −ri, siendo ri = b− A˜ xi el residuo en la aproximacion xi.
UNIVERSIDAD DE A CORUNA — GRUPO DE METODOS NUMERICOS EN INGENIERIA
Metodo de Gradientes Conjugados [CG] (Ib)
Si elegimos como vectores {vi}i=1,n los gradientes (cambiados de signo)
v1 = −∇f(x0) = r0 =(b− A˜ x0
)v2 = −∇f(x1) = r1 =
(b− A˜ x1
)v3 = −∇f(x2) = r2 =
(b− A˜ x2
). . .vk = −∇f(xk−1) = rk−1 =
(b− A˜ xk−1
). . .vn = −∇f(xn−1) = rn−1 =
(b− A˜ xn−1
)entonces sucede que (*)
βkj = −
sTj A˜ rk−1
sTj A˜ sj{
= 0, j = 1, . . . , k − 2
6= 0, j = k − 1
(*) Se comprueba que esto es ası aunque no es evidente (ver equivalencias).
UNIVERSIDAD DE A CORUNA — GRUPO DE METODOS NUMERICOS EN INGENIERIA
Metodo de Gradientes Conjugados [CG] (Ic)
FINALMENTE, LA BASE DE GRADIENTES CONJUGADOS ES (*):
{r0 = b−A˜ x0
s1 = r0,
}
rk−1 = b−A˜ xk−1
βkk−1 = −sTk−1A˜ rk−1
sTk−1A˜ sk−1,
sk = rk−1 + βkk−1sk−1,
k = 2, . . . , n.
(*) EN LO SUCESIVO PRESCINDIREMOS DEL SUPERINDICE k EN βkk−1
(ya que no es necesario).
UNIVERSIDAD DE A CORUNA — GRUPO DE METODOS NUMERICOS EN INGENIERIA
Metodo de Gradientes Conjugados [CG] (IIa)
PRIMERA ITERACION k = 0
Dado x0 −→ r0 = b−A˜ x0,
s1 = r0, −→ α1 =sT1 r0
sT1 A˜ s1,
x1 = x0 + α1s1.
UNIVERSIDAD DE A CORUNA — GRUPO DE METODOS NUMERICOS EN INGENIERIA
Metodo de Gradientes Conjugados [CG] (IIb)
ITERACIONES SIGUIENTES k = 1, . . . , n− 1
Dado xk −→ rk = b−A˜ xk, −→ βk = −sTk A˜ rksTk A˜ sk ,
sk+1 = rk + βksk, −→ αk+1 =sTk+1 rk
sTk+1 A˜ sk+1,
xk+1 = xk + αk+1sk+1,
Y xn verifica A˜ xn = b (SALVO ERRORES DE REDONDEO).
UNIVERSIDAD DE A CORUNA — GRUPO DE METODOS NUMERICOS EN INGENIERIA
Metodo de Gradientes Conjugados [CG] (IIc)
El algoritmo se detendra al finalizar el paso k. . .
I si k = n− 1, pues xk+1 = xn y por tanto verifica
rk+1 = b−A˜ xk+1 = 0 (salvo errores de redondeo *),
I o si ha convergido segun un criterio tipo{‖xk+1 − xk‖ ≤max{εx, rx ‖xk‖} y
‖rk+1‖ ≤max{εr, rr∥∥b∥∥}
(*) Se suelen realizar algunas iteraciones mas para refinar la solucion.
UNIVERSIDAD DE A CORUNA — GRUPO DE METODOS NUMERICOS EN INGENIERIA
Metodo de Gradientes Conjugados [CG] (IIIa)
Equivalencias
1) sTk+1A˜ sk+1 = sTk+1A˜ rk,pues sk+1 = rk + βksk,
y sTk+1A˜ sk = 0.
2) rk+1 = rk − αk+1A˜ sk+1,pues rk+1 = b− A˜ xk+1 = b− A˜
(xk + αk+1sk+1
)=(b− A˜ xk)− αk+1A˜ sk+1 = rk − αk+1A˜ sk+1.
3) sTk+1rk+1 = 0, LUEGO EL AVANCE ES OPTIMO (*),
pues sTk+1rk+1 = s
Tk+1
(rk − αk+1A˜ sk+1
)= s
Tk+1rk − αk+1
(sTk+1A˜ sk+1
)
y αk+1 =sTk+1rk
sTk+1
A˜ sk+1
.
(*) Vease el planteamiento de mınima energıa.
UNIVERSIDAD DE A CORUNA — GRUPO DE METODOS NUMERICOS EN INGENIERIA
Metodo de Gradientes Conjugados [CG] (IIIb)
4) A˜ sk = − 1αk
(rk − rk−1),pues rk = rk−1 − αkA˜ sk (por la equivalencia 2).
5) (rk − rk−1)Tsk = (rk − rk−1)
Trk−1,
pues sTk A˜ sk = s
Tk A˜ rk−1 (por la equivalencia 1) =⇒ s
Tk(A˜ sk) = r
Tk−1
(A˜ sk) ,
luego sTk
(−
1
αk
(rk − rk−1
))= r
Tk−1
(−
1
αk
(rk − rk−1
))(por la equivalencia 4).
6) sTk+1rk = rTk rk,
pues sk+1 = rk + βksk =⇒ sTk+1rk = r
Tk rk + βks
Tk rk,
y sTk rk = 0 (por la equivalencia 3).
UNIVERSIDAD DE A CORUNA — GRUPO DE METODOS NUMERICOS EN INGENIERIA
Metodo de Gradientes Conjugados [CG] (IIIc)
7) rTk+1rk = 0,
pues rTk+1rk =
(rk − αk+1A˜ sk+1
)Trk (por la equivalencia 2)
= rTk rk − αk+1
(sTk+1A˜ rk
)= r
Tk rk − αk+1
(sTk+1A˜ sk+1
)(por la equivalencia 1)
y αk+1 =sTk+1rk
sTk+1
A˜ sk+1
=rTk rk
sTk+1
A˜ sk+1
(por la equivalencia 6).
UNIVERSIDAD DE A CORUNA — GRUPO DE METODOS NUMERICOS EN INGENIERIA
Metodo de Gradientes Conjugados [CG] (IIId)
Luego hay muchas formulas equivalentes. . .
βk = −sTkA˜ rksTkA˜ sk = −
rTkA˜ sksTkA˜ sk = −
rTk (A˜ sk)sTk
(A˜ sk) = −rTk
(− 1αk
(rk − rk−1
))sTk
(− 1αk
(rk − rk−1
)) (por la equivalencia 4)
=
(rk − rk−1
)Trk
−(rk − rk−1
)Tsk
=
(rk − rk−1
)Trk
−(rk − rk−1
)Trk−1
(por la equivalencia 5)
=rTk rk
rTk−1
rk−1
(por la equivalencia 7)
=
(rk − rk−1
)Trk
rTk−1
rk−1
(por la equivalencia 7)
=rTk rk
−(rk − rk−1
)Tsk
(por la equivalencia 7).
UNIVERSIDAD DE A CORUNA — GRUPO DE METODOS NUMERICOS EN INGENIERIA
Metodo de Gradientes Conjugados [CG] (IV)
Formulas equivalentes para el calculo del coeficiente βk
βk = −sTk A˜ rksTk A˜ sk
βk =[rk − rk−1]T rk−[rk − rk−1]T sk
HESTENES–STIEFEL
βk =rTk rk
rTk−1 rk−1FLETCHER–REEVES
βk =[rk − rk−1]T rkrTk−1 rk−1
POLAK–RIBIERE
βk =rTk rk
−[rk − rk−1]T skMYERS
UNIVERSIDAD DE A CORUNA — GRUPO DE METODOS NUMERICOS EN INGENIERIA
Planteamiento de Mınima Energıa (I)
Sea el problema:
Hallar x,
que VERIFICA r = 0,
siendo r = b−A˜ x,con
A˜ = A˜T (A˜ SIMETRICA),
vTA˜ v > 0 ∀v 6= 0 (A˜ DEF. +).
El problema anterior tambien puede escribirse en la forma siguiente...
UNIVERSIDAD DE A CORUNA — GRUPO DE METODOS NUMERICOS EN INGENIERIA
Planteamiento de Mınima Energıa (II)
MINIMA ENERGIA
Hallar x,
que MINIMIZA f(x) =1
2xTA˜ x− xT b, (∗)
con
A˜ = A˜T (A˜ SIMETRICA),
vTA˜ v > 0 ∀v 6= 0 (A˜ DEF. +).
(*) Pues ∇f(x) =[dfdx
]T= −r, siendo r = b− A˜ x.
UNIVERSIDAD DE A CORUNA — GRUPO DE METODOS NUMERICOS EN INGENIERIA
Planteamiento de Mınima Energıa (IIIa)
AVANCE OPTIMO
Si aplicamos al problema anterior un metodo numerico tipo
xk+1 = xk + αk+1sk+1
diremos que el avance αk+1 es optimo cuando se cumpla
αk+1 = minα
Φ(α), con Φ(α) = f(x)
∣∣∣∣∣x=xk+αsk+1
,
o lo que es lo mismo . . .
UNIVERSIDAD DE A CORUNA — GRUPO DE METODOS NUMERICOS EN INGENIERIA
Planteamiento de Mınima Energıa (IIIb)
dΦ
dα
∣∣∣α=αk+1
= 0, siendodΦ
dα=df
dx
∣∣∣∣∣x=xk+αsk+1
sk+1
= −(b−A˜ (xk + αsk+1)
)Tsk+1.
Por tanto, la condicion de avance optimo es
rTk+1sk+1 = 0 ⇐⇒ sTk+1rk+1 = 0,
que tambien puede escribirse como
αk+1 =sTk+1rk
sTk+1A˜T sk+1
,
luego LOS METODOS DE DIRECCIONES CONJUGADAS SON DE AVANCE OPTIMO.
UNIVERSIDAD DE A CORUNA — GRUPO DE METODOS NUMERICOS EN INGENIERIA
Planteamiento de Mınima Energıa (IIIc)
AVANCE OPTIMO
UNIVERSIDAD DE A CORUNA — GRUPO DE METODOS NUMERICOS EN INGENIERIA
Precondicionamiento (I)
Sea el sistemaA˜ x = b.
Sea el precondicionador
P˜ = E˜E˜T , con det(E˜) 6= 0.
Reescribimos el sistema original en la forma(E˜−1
A˜E˜−T)
︸ ︷︷ ︸A˜
(E˜T x
)︸ ︷︷ ︸x
=(E˜−1
b)
︸ ︷︷ ︸b
,
y esperamos que al aplicar el Metodo de Gradientes Conjugados al sistemaprecondicionado
A˜ x = b
el proceso converja mas rapidamente a la solucion que en el caso delsistema original.
UNIVERSIDAD DE A CORUNA — GRUPO DE METODOS NUMERICOS EN INGENIERIA
Precondicionamiento (IIa)
PRIMERA ITERACION k = 0
Dado x0 = E˜T x0 −→ r0 = b− A˜ x0,
s1 = r0, −→ α1 =sT1 r0
sT1 A˜ s1
,
x1 = x0 + α1s1.
UNIVERSIDAD DE A CORUNA — GRUPO DE METODOS NUMERICOS EN INGENIERIA
Precondicionamiento (IIb)
ITERACIONES SIGUIENTES k > 0
Dado xk −→ rk = b− A˜ xk, −→ βk = −sTk A˜ rksTk A˜ sk ,
sk+1 = rk + βksk, −→ αk+1 =sTk+1 rk
sTk+1 A˜ sk+1
,
xk+1 = xk + αk+1sk+1,
Y xn = E˜−T xn verifica A˜ xn = b (SALVO ERRORES DE REDONDEO).
UNIVERSIDAD DE A CORUNA — GRUPO DE METODOS NUMERICOS EN INGENIERIA
Precondicionamiento (IIIa)
Relacion del problema precondicionado con el problema original
1)
xk = E˜T xk −→ rk = b− A˜E˜T xk= E˜−1
b−(E˜−1
A˜E˜−T)E˜T xk
= E˜−1b− E˜−1
A˜(E˜−TE˜T
)xk
= E˜−1b− E˜−1
A˜ xk= E˜−1 (
b− A˜ xk)= E˜−1
rk, con rk = b− A˜ xk.
UNIVERSIDAD DE A CORUNA — GRUPO DE METODOS NUMERICOS EN INGENIERIA
Precondicionamiento (IIIb)
2)
xk+1 = E˜T xk+1 −→ xk+1 = E˜−T (xk + αk+1sk+1)
=(E˜−TE˜T
)xk + αk+1E˜−T sk+1
= xk + αk+1sk+1, con sk+1 = E˜−T sk+1.
sk+1 = E˜−T(rk + βksk
)= E˜−T
(E˜−1
rk + βkE˜T sk)
=(E˜E˜T
)−1
rk + βksk
= P˜−1rk + βksk.
UNIVERSIDAD DE A CORUNA — GRUPO DE METODOS NUMERICOS EN INGENIERIA
Precondicionamiento (IIIc)
3)
αk+1 =sTk+1rk
sTk+1A˜ sk+1
=
(E˜T sk+1
)T (E˜−1rk
)(E˜T sk+1
)T (E˜−1A˜E˜−T) (E˜T sk+1
)=
sTk+1
(E˜E˜−1
)rk
sTk+1
(E˜E˜−1
)A˜ (E−T˜ E˜T) sk+1
=sTk+1rk
sTk+1A˜ sk+1
.
UNIVERSIDAD DE A CORUNA — GRUPO DE METODOS NUMERICOS EN INGENIERIA
Precondicionamiento (IIId)
4)
βk = −sTk A˜ rksTk A˜ sk
= −(E˜T sk)T (E˜−1A˜E˜−T) (E˜−1rk
)(E˜T sk)T (E˜−1A˜E˜−T) (E˜T sk)
= −sTk(E˜E˜−1
)A˜ (E−T˜ E˜−1
)rk
sTk(E˜E˜−1
)A˜ (E˜−TE˜T) sk
= −sTkA˜ (E˜E˜T)−1
rk
sTkA˜ sk= −
sTkA˜P˜−1rk
sTkA˜ sk .
UNIVERSIDAD DE A CORUNA — GRUPO DE METODOS NUMERICOS EN INGENIERIA
Gradiente Conjugado Precondicionado [PCG] (Ia)
PRIMERA ITERACION k = 0
Dado x0 −→ r0 = b−A˜ x0,
s1 = P˜−1r0, −→ α1 =sT1 r0
sT1 A˜ s1,
x1 = x0 + α1s1.
UNIVERSIDAD DE A CORUNA — GRUPO DE METODOS NUMERICOS EN INGENIERIA
Gradiente Conjugado Precondicionado [PCG] (Ib)
ITERACIONES SIGUIENTES k > 0
Dado xk −→ rk = b−A˜ xk, −→ βk = −sTk A˜ P˜−1 rksTk A˜ sk ,
sk+1 = P˜−1rk + βksk, −→ αk+1 =sTk+1 rk
sTk+1 A˜ sk+1,
xk+1 = xk + αk+1sk+1,
Y xn verifica A˜ xn = b (SALVO ERRORES DE REDONDEO).
UNIVERSIDAD DE A CORUNA — GRUPO DE METODOS NUMERICOS EN INGENIERIA
Gradiente Conjugado Precondicionado [PCG] (Ic)
El algoritmo se detendra al finalizar el paso k. . .
I si k = n− 1, pues xk+1 = xn y por tanto verifica
rk+1 = b−A˜ xk+1 = 0 (salvo errores de redondeo *),
I o si ha convergido segun un criterio tipo{‖xk+1 − xk‖ ≤max{εx, rx ‖xk‖} y
‖rk+1‖ ≤max{εr, rr∥∥b∥∥}
(*) Se suelen realizar algunas iteraciones mas para refinar la solucion.
UNIVERSIDAD DE A CORUNA — GRUPO DE METODOS NUMERICOS EN INGENIERIA
Gradiente Conjugado Precondicionado [PCG] (II)
Observaciones:
• Los sk son A˜–conjugados, pues
sTi A˜ sj =
(E˜−T si
)T (E˜A˜E˜T
)(E˜−T sj
)= s
Ti
(E˜−1
E˜)A˜(E˜TE˜−T
)sj
= sTi A˜ sj.
• El avance es optimo, pues
sTk+1rk+1 =
(E˜−T sk+1
)T(E˜ rk+1) = s
Tk+1
(E˜−1
E˜)rk+1
= sTk+1rk+1.
• En generalrTk+1rk = (E˜ rk+1)
T(E˜ rk) = r
Tk+1
(E˜TE˜
)rk
6= 0.
UNIVERSIDAD DE A CORUNA — GRUPO DE METODOS NUMERICOS EN INGENIERIA
Gradiente Conjugado Precondicionado [PCG] (III)
CONCLUSIONES:
♥ Suele utilizarse el precondicionador diagonal P˜ = D˜A.
♣ Para calcular xn el tiempo de computacion es T (n3).♣ NO ES MEJOR QUE LOS METODOS DIRECTOS (GAUSS, CHOLESKY), A MENOS QUE
CONVERJA ANTES.
♥ Es un metodo directo, pero funciona como los metodos iterativos ya que♥ PUEDE DETENERSE ANTICIPADAMENTE, LO QUE PROPORCIONA UN APROXIMACION A
LA SOLUCION.♦ NO MODIFICA LOS COEFICIENTES DE LA MATRIZ,♦ POR LO QUE PUEDE SER UTILIZADO CON ESQUEMAS DE ALMACENAMIENTO MINIMO
PARA MATRICES DISPERSAS.
♣ Puede ser muy ventajoso (o la unica alternativa posible) cuando lamatriz es no-estructurada y muy dispersa.
♦ EN TALES CASOS SUELE UTILIZARSE SIN ENSAMBLAR LA MATRIZ.
UNIVERSIDAD DE A CORUNA — GRUPO DE METODOS NUMERICOS EN INGENIERIA
Implementacion del Algoritmo PCG (Ia)
INICIALIZACION (k = 0)
Inicializamos:x0 → r0 = b− A˜ x0,
w0 = P˜−1r0 → s1 = w0
z1 = −A˜ s1 → s1 = −zT1 s1 → α1 = rT0 s1/s1
∆x0 = α1s1 → x1 = x0 + ∆x0 → r1 =
r0 + α1z1 (*)
b− A˜ x1 (**)
cnvgc=.false.
k = 0
(*) Pues r1 = b− A˜ x1 = b− A˜(x0 + α1s1) = (b− A˜ x0)− α1A˜ s1 = r0 + α1z1.
(**) De vez en cuando se recalcula el residuo para evitar la acumulacion de errores de redondeo.
UNIVERSIDAD DE A CORUNA — GRUPO DE METODOS NUMERICOS EN INGENIERIA
Implementacion del Algoritmo PCG (Ib)
ITERACIONES SIGUIENTES (k > 0)
do while (k.lt.niter.and.not.cnvgc)
wk = P˜−1rk → βk = z
Tk wk/sk → sk+1 = wk + βksk
zk+1 = −A˜ sk+1 → sk+1 = −zTk+1sk+1 → αk+1 = rTk sk+1/sk+1
∆xk = αk+1sk+1 → xk+1 = xk + ∆xk → rk+1 =
rk + αk+1zk+1 (*)
b− A˜ xk+1 (**)
cnvgc=(‖∆xk‖.le. max(εx, rx ‖xk‖)) .and. (‖rk+1‖.le. max(εr, rr∥∥b∥∥))
k = k + 1
enddo
(*) Pues rk+1 = b− A˜ xk+1 = b− A˜(xk + αk+1sk+1) = (b− A˜ xk)− αk+1A˜ sk+1 = rk + αk+1zk+1.
(**) De vez en cuando se recalcula el residuo para evitar la acumulacion de errores de redondeo.
UNIVERSIDAD DE A CORUNA — GRUPO DE METODOS NUMERICOS EN INGENIERIA
Implementacion del Algoritmo PCG (II)
TRATAMIENTO DE LAS CONDICIONES DE VINCULACION
Dado el sistema A˜ x = b+ R, con algunas xv = pv, se procede de lasiguiente manera:
1) Se inicializan los valores de los g.d.l. prescritos, xv = pv.
2) Cuando se calcula el residuo r = b− A˜ x se opera con toda la matriz.
3) Al actualizar los valores de x en cada iteracion se ignoran (no se modifican) lasincognitas correspondientes a los g.d.l. prescritos, xv = pv.
4) Al comprobar la condicion de convergencia en r hay que tener en cuenta que losterminos corespondientes a los g.d.l. prescritos no tienden a 0 sino a los valores delas reacciones correspondientes.
5) Finalmente, las reacciones se obtienen a partir del ultimo residuo, ya que R = −r.
UNIVERSIDAD DE A CORUNA — GRUPO DE METODOS NUMERICOS EN INGENIERIA