PDF de programación - Implementación en Fortran95 de los métodos iterativos básicos para la solución de sistemas de ecuaciones lineales

Implementación en Fortran95 de los métodos iterativos básicos

para la solución de sistemas de ecuaciones lineales

Sea

Pablo Santamaría
v0.1 (Mayo 2014)

A x = b,

un sistema de n ecuaciones lineales con n incógnitas, donde A es una matriz cuadrada de orden n de elementos
reales y no singular (por lo cual el sistema admite una solución y ésta es única). Un método iterativo para
resolver dicho sistema comienza con una aproximación inicial x(0) a la solución x, y genera una sucesión de
vectores {x(k)}∞
k=0 que se espera converja a x. Los dos métodos iterativos básicos son el método de Jacobi y
el método de Gauss–Seidel cuyas implementaciones algorítmicas se describen a continuación.

Método de Jacobi
Dada A de n × n con elementos diagonales no nulos.
Escoger x(0) de orden n, usualmente x(0) = 0.
Para k = 1, 2, . . .

bi −Pn

aijx

(k−1)
j

j=1
(j6=i)
aii

Tomar x

(k)
i =

, para i = 1, . . . , n.

Método de Gauss–Seidel
Dada A de n × n con elementos diagonales no nulos.
Escoger x(0) de orden n, usualmente x(0) = 0.
Para k = 1, 2, . . .

i = bi −Pi−1

(k)

j=1 aijx

j −Pn

(k)

aii

Tomar x

j=i+1 aijx

(k−1)
j

, para i = 1, . . . , n.

Como en todo método iterativo, debe establecerse un criterio de paro para las iteraciones, el cual, en este

contexto, puede ser

kx(k)k
respecto de una norma y cierta tolerancia
prefijada.

kx(k) − x(k−1)k

<
,

Nótese que en el método de Jacobi cada componente de una nueva iteración es generada a partir de las
componentes de la iteración anterior en forma independiente y simultáneamente, de aquí que el método
se conozca también como método de desplazamientos simultáneos. En el método de Gauss-Seidel, por otra
parte, las componentes de una nueva iteración son sucesivamente actualizadas tan pronto una es conocida, y
por ello al método se lo conoce también como método de desplazamientos sucesivos.

Desde el punto de vista de la implementación computacional, lo anterior implica que el método de Jacobi
requiere del almacenamiento de dos vectores, el de la iteración actual y el de la anterior, en tanto que el
método de Gauss-Seidel permite considerar un solo vector cuyas componentes se actualizan al momento de
la iteración. Por otra parte, la naturaleza simultánea de las correcciones en el método de Jacobi, hacen del

1

mismo un método inherentemente paralelizable, mientras que, en principio, el método de Gauss-Seidel es
secuencial.

La naturaleza paralelizable del método de Jacobi queda manifiesta en Fortran 95 donde podemos manipular
los arreglos como un todo. En efecto, reescribiendo las ecuaciones de la iteración k-ésima del método como
sigue

(k)
i =

x

bi −Pn

(k−1)
j

aijx

j=1
(j6=i)
aii
(k−1)
i

−Pn
+ bi −Pn

aii

= bi + aiix

= x

(k−1)
i

(k−1)
j=1 aijx
j

(k−1)
j

,

j=1 aijx

aii

notemos quePn

j=1 aijx
iteración procede como

(k−1)
j

es la componente i-ésima del producto matricial A x(k−1), y por lo tanto la

(k)
i = x

x

(k−1)
i

= x

(k−1)
i

+ bi − (A x(k−1))i
+ (b − A x(k−1))i

aii

aii

,

para i = 1, . . . , n. Si introducimos el vector d = (a11, a22, . . . , ann)t, podemos escribir, finalmente, la iteración
k-ésima del método de Jacobi como

x(k) = x(k−1) + b − A x(k−1)

d

.

donde la “división” de arreglos se interpreta componente a componente.

Por otra parte, en las ecuaciones de la iteración k-ésima del método de Gauss-Seidel

i = bi −Pi−1

(k)

x

j=1 aijx

j −Pn

(k)

aii

(k−1)
j=i+1 aijx
j

,

la primer sumatoria puede pensarse como el producto escalar del vector de i−1 elementos (ai1, ai2, . . . , aii−1)t
conformado por las primeras i−1 columnas de la fila i-ésima, por el vector de i−1 elementos (x
(k)
i−1)t
conformado por los valores obtenidos hasta el momento en la iteración. Análogamente, la segunda sumatoria
puede escribirse como el producto escalar del vector (aii+1, aii+2, . . . , ain)t por el vector conformado por
valores obtenidos en la iteración anterior (x

(k)
2 , . . . , x

(k)
1 , x

, . . . , x

)t.

(k−1)
n

, x

(k−1)
i+1

(k−1)
i+2

El siguiente módulo Fortran implementa sendos métodos según lo expuesto, considerando la norma

máxima para el criterio de paro de las iteraciones.

2

Código 1. Implementación de los métodos iterativos de Jacobi y Gauss-Seidel

MODULE iterative_methods

USE, INTRINSIC :: iso_fortran_env, ONLY: WP => REAL64
IMPLICIT NONE

CONTAINS

SUBROUTINE iterative_jacobi(a,b,x,tol,iter,clave)

! -----------------------------------------------------------------
! Argumentos de la subrutina
! -----------------------------------------------------------------
REAL(WP), INTENT(IN)
:: a(:,:) ! Matriz del sistema
REAL(WP), INTENT(IN)
:: b(:)
REAL(WP), INTENT(INOUT) :: x(:)
REAL(WP), INTENT(IN)
INTEGER,
INTEGER,

! Término independiente
! Aproximación inicial/solución
! Tolerancia para el error
! Max iteraciones/iter realizadas
! Clave de éxito/error
! clave = 0, OK.
! clave < 0, error en asig. memoria
! clave > 0, max iter. alcanzado

INTENT(INOUT) :: iter
INTENT(OUT)

:: clave

:: tol

! -----------------------------------------------------------------
! Variables locales
! -----------------------------------------------------------------
REAL(WP), ALLOCATABLE :: diag(:) ! Arreglo para guardar la

REAL(WP), ALLOCATABLE :: x0(:)

! diagonal de la matriz del sistema
! Arreglo para guardar la
! iteración anterior
! Contador de iteraciones
! Indice

:: k
:: i

INTEGER
INTEGER
! -----------------------------------------------------------------
! Procedimiento
! -----------------------------------------------------------------
ALLOCATE(diag(SIZE(a,1)),x0(SIZE(a,1)), STAT=clave)
IF (clave /=0 ) THEN

clave = -1
RETURN

ENDIF

clave = 1
diag
DO k=1,iter

= [ (a(i,i), i=1,SIZE(a,1)) ]

x0 = x
x
IF ( MAXVAL(ABS(x-x0)) <= MAXVAL(ABS(x))*tol ) THEN

= x0 + (b - MATMUL(a,x0))/diag

clave = 0
iter
= k
EXIT

ENDIF

ENDDO

DEALLOCATE(diag,x0)

END SUBROUTINE iterative_jacobi

SUBROUTINE iterative_gauss_seidel(a,b,x,tol,iter,clave)

! -----------------------------------------------------------------
! Argumentos de la subrutina
! -----------------------------------------------------------------
REAL(WP), INTENT(IN)
:: a(:,:) ! Matriz del sistema
REAL(WP), INTENT(IN)
:: b(:)
REAL(WP), INTENT(INOUT) :: x(:)
REAL(WP), INTENT(IN)
INTEGER,

! Término independiente
! Aproximación inicial/solución
! Tolerancia para el error
! Max iteraciones/iter realizadas

INTENT(INOUT) :: iter

:: tol

3

INTEGER,

INTENT(OUT)

= 0, OK.

:: clave

! Clave de éxito/error
! clave
! clave /= 0, max iter. alcanzado
! -----------------------------------------------------------------
! Variables locales
! -----------------------------------------------------------------
INTEGER
INTEGER
INTEGER
REAL(WP) :: xi
REAL(WP) :: xnorma
REAL(WP) :: difnorma

! Dimensión del problema
! Contador de iteraciones
! Indice
! Componente del vector iterado
! Norma del vector iterado
! Norma de la diferencia entre dos
! iteraciones del vector

:: n
:: k
:: i

! -----------------------------------------------------------------
! Procedimiento
! -----------------------------------------------------------------
clave = 1
n
DO k=1,iter

= SIZE(a,1)

= 0.0_WP
xnorma
difnorma = 0.0_WP
DO i=1,n

xi

&

= (b(i) - DOT_PRODUCT(a(i,1:i-1),x(1:i-1)) &

- DOT_PRODUCT(a(i,i+1:n),x(i+1:n)))/a(i,i)

= MAX(xnorma,ABS(xi))

xnorma
difnorma = MAX(difnorma,ABS(xi-x(i)))
x(i)

= xi

END DO
IF (difnorma <= xnorma*tol) THEN

iter
= k
clave = 0
EXIT

ENDIF

ENDDO

END SUBROUTINE iterative_gauss_seidel

END MODULE iterative_methods

4