MÉTODO ITERATIVO (JACOBI)

Cuando hay muchas ecuaciones y , sobre todo , si hay muchos 0-s y no se puede resolver con el método directo, se empleará el MÉTODO ITERATIVO.

Los métodos iterativos se basan en la siguiente idea: La ecuación f(x)=0 se modificará para obtener ! X(n+1)=g(Xn)

Teniendo el siguiente sistema de n ecuaciones y n incógnitas:

a11x1+a12x2+...+a1nxn=b1

: : ! Ax=b

: :

an1x1+an2x2+...+annxn=bn

Se quiere calcular una solución para el sistema definido, para ello , consideramos A una matriz en cuya diagonal no tenga ningún elemento nulo.

Para poder realizar el calculo es recomendable dividir la matriz A en dos diferentes matrices: A= D+E

Considerando un sistema de 3 ecuaciónes:

a11 a12 a13 a11 0 0 0 a12 a13

a21 a22 a23 = 0 a22 0 + a21 0 a23

a31 a32 a33 0 0 a33 a31 a32 0

! ! !

A D E

Sustituyendo A por (D+E) en la ecuación Ax=b obtendremos la siguiente

manera de representar la ecuación (D+E)x=b. Transformando la ecuación a la siguiente forma (Dx=b-Ex) :

a11x1= b1-a12x2-a13x3

a22x2= b2-a21x1-a23x3

a33x3= b3-a31x1-a32x2

!

Términos de la diagonal

-1

(D+E)x=b ! Dx=b-Ex ! x=D (b-Ex)

-1 1/a11 0 0

Siendo D = 0 1/a22 0

0 0 1/a33

Sustituyendo de la manera que hemos deducido anteriormente, obtendremos el siguiente sistema de ecuaciones:

X1=1/a11 (b1-a12x2-a13x3)

X2=1/a22 (b2-a21x1-a23x3)

X3=1/a33 (b3-a31x1-a32x2)

Generalizando para “n” ecuaciónes:

(m+1) (m) (m)

x1 =1/a11 (b1-a12x2 -...-a1nxn )

(m+1) (m) (m)

X2 =1/a22 (b2-a21x21-...-a2nxn )

:

(m+1) (m) (m)

Xn =1/ann (bn-an1x1-...-ann-1xn-1 )

(m+1) -1 (m)

Representándolo en forma matricial: x = D (b-Ex ) !

(m+1) n (m)

! xi = 1/aii (bi- aij xj ) (i=1,2,3...n)

j=1

i"j

A continuación el programa que hemos realizado para calcular el método de Jacobi:

program jacobi;

type

lista=array [1..10,1..10] of real;

lista1=array [1..10] of real;

var

l:lista;

l1,l2:lista1;

fila,columna,i:integer;

error:real;

procedure introducir_datos (var l:lista;fila,columna:integer);

var

i,j:integer;

begin

for i:=1 to fila do

for j:=1 to columna do

begin

writeln ('Introduce valor de la fila ',i,' y de la columna ',j,':');

readln (l[i,j]);

end;

end;

function dominante (l:lista;fila,columna:integer):boolean;

var

cont:real;

i,j:integer;

resultado:boolean;

begin

resultado:=true;

for i:=1 to fila do

begin

cont:=0.0;

for j:=1 to (columna-1) do

if (j<>i) then

cont:=cont+abs(l[i,j]);

if (abs(l[i,i])>cont) and (resultado) then

resultado:=true

else

resultado:=false;

end;

dominante:=resultado;

end;

procedure calcular (l:lista;var l1:lista1;fila,columna:integer;error:real);

var

cont,i,j:integer;

resultado1,resultado:boolean;

begin

cont:=1;

if (dominante(l,fila,columna)) then

begin

resultado1:=false;

while (not resultado1) do

begin

for i:=1 to fila do

l2[i]:=0.0;

for i:=1 to fila do

begin

for j:=1 to columna do

if (i<>j) and (j<>columna) then

l2[i]:=l2[i]+((-1)*(l[i,j]*l1[j]))

else

if (i<>j) and (j=columna) then

l2[i]:=l2[i]-(l[i,j]);

l2[i]:=l2[i]/l[i,i];

end;

resultado:=true;

for i:=1 to fila do

if (abs((l2[i]-l1[i]))< error) and (resultado) then

begin

resultado:=true;

resultado1:=true;

end

else

begin

resultado:=false;

resultado1:=false;

end;

if (not resultado) then

begin

for i:=1 to fila do

l1[i]:=l2[i];

cont:=cont+1;

end;

end;

writeln ('EL NUMERO DE INTERACCIONES: ',cont);

for i:=1 to fila do

writeln ('EL VALOR DE X',i,' ES: ',l2[i]);

end

else

writeln ('ES DIVERGENTE');

end;

begin

writeln ('Introduce el numero de filas:');

readln (fila);

writeln ('Introduce el numero de columnas:');

readln (columna);

writeln ('Introduce el error:');

readln (error);

for i:=1 to fila do

begin

writeln ('Introduce el valor inicial de X',i,' :');

readln (l1[i]);

end;

introducir_datos (l,fila,columna);

calcular (l,l1,fila,columna,error);

end.

A continuación introduciremos tres ejemplos diferentes, analizando sus respectivas soluciones:

1/ x1+6x2-3x3=4 (0) (0) (0)

5x1+2x2-x3=6 ! x1=x2=x3=0 ! Error máximo= 0.005

2x1+x2+4x3=7

Solución: ES DIVERGENTE

Este problema no se puede resolver del modo que está planteado, para obtener una solución la matriz A tiene que ser diagonalmente dominante. Para ello debe cumplir la siguiente condición:

n

aii >  aij ( i=1...n)

j=1

j"i

2/ Si cambiamos el orden de las ecuaciones para que se cumpla dicha condición, obtendremos unos resultados.

Intercambiaremos las dos primeras ecuaciones obteniendo el siguiente sistema:

5x1+2x2-x3=6

x1+6x2-3x3=4

2x1+x2+4x3=7

Solución: NÚMERO DE INTERACCIONES: 11

EL VALOR DE X1 ES: 1.0007418833

EL VALOR DE X2 ES: 0.9998272877

EL VALOR DE X3 ES: 1.0017926613

3/ A continuación introduciremos un sistema de ecuaciones al azar, teniendo en cuenta que tiene que ser la matriz diagonalmente dominante:

5x1+x2-0.5x3=2 (0) (0) (0)

3x1+7.25x2+2.5x3=4 ! x1=x2=x3= 0 ! Error máximo= 0.005

-0.75x1-3.7x2+9x3=-7.5

Solución: NÚMERO DE INTERACCIONES: 5

EL VALOR DE X1 ES: 0.21537555233

EL VALOR DE X2 ES: 0.65072165553

EL VALOR DE X3 ES: -0.54767634828

MÉTODO DIRECTO (GAUSS)

Muchos de los problemas requieren de la solución de sistemas de ecuaciones lineales de la forma:

a11x1+a12x2+a13x3+...+a1nxn=b1

a21x1+a22x2+a23x3+...+a2nxn=b2

... ... ... ... ...

an1x1+an2x2+an3x3+...+annxn=bn

Donde aij y bi son las constantes, siendo xi los términos a calcular.

Mediante el método directo debemos procurar transformar el anterior sistema de ecuaciones a la siguiente forma:

A11x1+0+0+...+0= b1

0+a22x2+0+...+0= b2 ! creando un sistema, cuyos valores son

... ... nulos excepto en la diagonal.

0+0+0+...+annxn= bn

Como es muy complicado crear esta diagonal, se opta por ponerlo de forma triangular:

a11x1+a12x2+a13x3+...+a1nxn=b1

0 +a22x2+a23x3+...+a2nxn=b2

0 + 0 +a33x3+...+a3nxn=b3

... ... ...

0 + 0 + 0 +...+annxn=bn ! xn=bn/ann

Tal y como muestra la flecha , se debe resolver hacia atrás el sistema triangular superior.

Generalizando: xn-1= (b(n-1) - a(n-1)nxn) / a(n-1)(n-1)

!

MÉTODO DE GAUSS

Se debe tener en cuenta, que las ecuaciones deben ser independientes las unas de las otras, es decir, el determinante debe ser diferente a 0.

En primera instancia los pasos a seguir para lograr la solución serán:

1/ Triangularizar superiormente el sistema.

2/ Resolver hacia atrás el sistema triangular superior.

El programa diseñado para ello es el siguiente:

program metodo_directo;

type

lista=array [1..10,1..10] of real;

lista1=array [1..10] of real;

var

l:lista;

l1:lista1;

fila,columna,i,j:integer;

procedure introducir_datos (var l:lista;fila,columna:integer);

var

i,j:integer;

valor:real;

begin

for i:=1 to fila do

for j:=1 to columna do

begin

writeln ('Introduce el valor de la fila ',i,' y de la columna ',j);

readln (valor);

l[i,j]:=valor;

end;

end;

procedure hacer_ceros (var l:lista;fila,columna:integer);

var

j,i,u:integer;

resultado:real;

begin

for j:=1 to (fila-1) do

for i:=(j+1) to fila do

begin

resultado:=(l[i,j]/l[j,j]);

for u:=1 to columna do

l[i,u]:=(l[i,u]-(resultado*l[j,u]));

end;

end;

procedure calcular (l:lista;var l1:lista1;fila,columna:integer);

var

cont,resultado:real;

begin

for i:=1 to fila do

l1[i]:=1;

for i:=fila downto 1 do

begin

cont:=0.0;

resultado:=0.0;

for j:=(columna-1) downto (i+1) do

if (l[i,j]<>0.0) then

cont:=cont+(l1[j]*l[i,j]);

resultado:=((-1)*(l[i,columna])-cont)/l[i,i];

l1[i]:=resultado;

end;

end;

begin

writeln ('Introduce el numero de filas de la matriz:');

readln (fila);

writeln ('Introduce el numero de columnas de la matriz:');

readln (columna);

introducir_datos (l,fila,columna);

for i:=1 to fila do

begin

for j:=1 to columna do

write (l[i,j],' ');

writeln;

end;

hacer_ceros (l,fila,columna);

writeln ('La matriz despues de hacer los ceros:');

for i:=1 to fila do

begin

for j:=1 to columna do

write (l[i,j],' ');

writeln;

end;

calcular (l,l1,fila,columna);

for i:=1 to fila do

writeln ('El resultado de la X',i,' es: ',l1[i]);

end.

A continuación para comprobar la validez del programa realizaremos dos ejemplos:

1/ 2x1+x2+3x3=11

4x1+3x2+10x3=28

2x1+4x2+17x3=31

Solución: EL VALOR DE X1 ES=3

EL VALOR DE X2 ES=2

EL VALOR DE X3 ES=1

2/ 0.610x1+1.23x2+1.72x3=0.792

1.02x1+2.15x2-5.51x3=12.0

-4.34x1+11.2x2-4.25x3=16.3

Solución: EL VALOR DE X1 ES=1.6074749447

EL VALOR DE X2 ES=1.6019476216

EL VALOR DE X3 ES=-1.2552065644

¿ Y si intercambiásemos la tercera ecuación con la primera ?

-4.34x1+11.2x2-4.25x3=16.3

1.02x1+2.15x2-5.51x3=12.0

0.610x1+1.23x2+1.72x3=0.792

Solución: EL VALOR DE X1 ES=1.6074749448

EL VALOR DE X2 ES=1.6019476216

EL VALOR DE X3 ES=-1.2552065644

En nuestro caso, la diferencia es insignificante, pero si utilizáramos pocos dígitos significativos, el error cometido en el caso anterior sería enorme (casi del 100%), ya que la magnitud del error se multiplica por el mismo.

Para evitarlo, debemos colocar el coeficiente de mayor módulo en lo más alto.

Por lo tanto hay que volver a modificar el programa para que en cada fila que haya que hacer 0 se compare y se ponga el de mayor módulo en lo más alto: ( MÉTODO DE GAUSS CON PIVOTACIÓN PARCIAL )

program metodo_directo;

type

lista=array [1..9,1..10] of real;

lista1=array [1..10] of real;

var

l:lista;

l1:lista1;

fila,columna,i,j:integer;

procedure introducir_datos (var l:lista;fila,columna:integer);

var

i,j:integer;

valor:real;

begin

for i:=1 to fila do

for j:=1 to columna do

begin

writeln ('Introduce el valor de la fila ',i,' y de la columna ',j);

readln (valor);

l[i,j]:=valor;

end;

end;

procedure pivotear (var l:lista;fila,columna:integer;posfila:integer);

var

aux:lista1;

i,j:integer;

begin

for i:=posfila to fila do

if (abs(l[posfila,posfila])<abs(l[i,posfila])) then

for j:=1 to columna do

begin

aux[j]:=l[i,j];

l[i,j]:=l[posfila,j];

l[posfila,j]:=aux[j];

end;

end;

procedure hacer_ceros (var l:lista;fila,columna:integer);

var

j,i,u:integer;

resultado:real;

begin

for j:=1 to (fila-1) do

begin

pivotear (l,fila,columna,j);

for i:=(j+1) to fila do

begin

resultado:=(l[i,j]/l[j,j]);

for u:=1 to columna do

l[i,u]:=(l[i,u]-(resultado*l[j,u]));

end;

end;

end;

procedure calcular (l:lista;var l1:lista1;fila,columna:integer);

var

cont,resultado:real;

begin

for i:=1 to fila do

l1[i]:=1;

for i:=fila downto 1 do

begin

cont:=0.0;

resultado:=0.0;

for j:=(columna-1) downto (i+1) do

if (l[i,j]<>0.0) then

cont:=cont+(l1[j]*l[i,j]);

resultado:=((-1)*(l[i,columna])-cont)/l[i,i];

l1[i]:=resultado;

end;

end;

begin

writeln ('Introduce el numero de filas de la matriz:');

readln (fila);

writeln ('Introduce el numero de columnas de la matriz:');

readln (columna);

introducir_datos (l,fila,columna);

for i:=1 to fila do

begin

for j:=1 to columna do

write (l[i,j],' ');

writeln;

end;

hacer_ceros (l,fila,columna);

writeln ('La matriz despues de hacer los ceros:');

for i:=1 to fila do

begin

for j:=1 to columna do

write (l[i,j],' ');

writeln;

end;

calcular (l,l1,fila,columna);

for i:=1 to fila do

writeln ('El resultado de la X',i,' es: ',l1[i]);

end.

Nos basaremos en un sencillo ejemplo para explicar en que se basa y por qué es necesario el escalamiento:

Si tenemos ! x1+2x2=5

0.1x1+10x2=5.4

!

Estaría correctamente posicionado según la pivotación de GAUSS

Sin embargo si multiplicásemos por 100 la segunda ecuación obtendríamos:

x1+2x2=5

10x1+1000x2=540

!

Sería necesaria la modificación de posiciones

Para evitar estas ambigüedades, se deben poner todas las ecuaciones en el mismo orden de magnitud, es decir , en orden unitario.

Para ello se debe buscar el módulo más grande y el menor de cada ecuación y realizar la media geométrica (c*d). Una vez realizada, se divide cada ecuación por su media geométrica, logrando un orden unitario y evitando el problema anteriormente comentado. Esto se denomina escalamiento.

NOTA: El elemento bi no se tendrá en cuenta a hora de buscar el máximo y mínimo de cada ecuación, pero sí a la hora de dividirlo por la media.

Para solucionar el problema, y poner todas las ecuaciones en el mismo orden de magnitud, al programa anteriormente citado, le añadiremos la función que se denomina media_geometrica:

program metodo_directo;

type

lista=array [1..9,1..10] of real;

lista1=array [1..10] of real;

var

l:lista;

l1:lista1;

fila,columna,i,j:integer;

procedure introducir_datos (var l:lista;fila,columna:integer);

var

i,j:integer;

valor:real;

begin

for i:=1 to fila do

for j:=1 to columna do

begin

writeln ('Introduce el valor de la fila ',i,' y de la columna ',j);

readln (valor);

l[i,j]:=valor;

end;

end;

procedure media_geometrica (var l:lista;fila,columna:integer);

var

i,j:integer;

max,min,media:real;

begin

for i:=1 to fila do

begin

max:=abs(l[i,1]);

min:=abs(l[i,1]);

for j:=2 to (columna-1) do

if (abs(l[i,j])>max) then

max:=abs(l[i,j])

else

if (abs(l[i,j])<min) then

min:=abs(l[i,j]);

media:=sqrt(max*min);

for j:=1 to columna do

l[i,j]:=l[i,j]/media;

end;

end;

procedure pivotear (var l:lista;fila,columna:integer;posfila:integer);

var

aux:lista1;

i,j:integer;

begin

for i:=posfila to fila do

if (abs(l[posfila,posfila])<abs(l[i,posfila])) then

for j:=1 to columna do

begin

aux[j]:=l[i,j];

l[i,j]:=l[posfila,j];

l[posfila,j]:=aux[j];

end;

end;

procedure hacer_ceros (var l:lista;fila,columna:integer);

var

j,i,u:integer;

resultado:real;

begin

for j:=1 to (fila-1) do

begin

pivotear (l,fila,columna,j);

for i:=(j+1) to fila do

begin

resultado:=(l[i,j]/l[j,j]);

for u:=1 to columna do

l[i,u]:=(l[i,u]-(resultado*l[j,u]));

end;

end;

end;

procedure calcular (l:lista;var l1:lista1;fila,columna:integer);

var

cont,resultado:real;

begin

for i:=1 to fila do

l1[i]:=1;

for i:=fila downto 1 do

begin

cont:=0.0;

resultado:=0.0;

for j:=(columna-1) downto (i+1) do

if (l[i,j]<>0.0) then

cont:=cont+(l1[j]*l[i,j]);

resultado:=((-1)*(l[i,columna])-cont)/l[i,i];

l1[i]:=resultado;

end;

end;

begin

writeln ('Introduce el numero de filas de la matriz:');

readln (fila);

writeln ('Introduce el numero de columnas de la matriz:');

readln (columna);

introducir_datos (l,fila,columna);

for i:=1 to fila do

begin

for j:=1 to columna do

write (l[i,j],' ');

writeln;

end;

media_geometrica (l,fila,columna);

writeln ('La matriz despues de calcular la media geometrica:');

for i:=1 to fila do

begin

for j:=1 to columna do

write (l[i,j],' ');

writeln;

end;

hacer_ceros (l,fila,columna);

writeln ('La matriz despues de hacer los ceros:');

for i:=1 to fila do

begin

for j:=1 to columna do

write (l[i,j],' ');

writeln;

end;

calcular (l,l1,fila,columna);

for i:=1 to fila do

writeln ('El resultado de la X',i,' es: ',l1[i]);

end.

Para comprobar el correcto funcionamiento del programa realizado, introduciremos varios ejemplos, analizando los resultados obtenidos:

1/ -0.04x1+0.04x2+0.12x3=3

0.56x1-1.56x2+0.32x3=1

-0.24x1+1.24x2-0.28x3=0

Solución: EL VALOR DE X1 ES=7

EL VALOR DE X2 ES=7

EL VALOR DE X3 ES=25

2/ 5x1+x2-0.5x3=2

3x1+7.25x2+2.5x3=4

-0.75x1-3.7x2+9x3=-7.5

Solución: EL VALOR DE X1 ES=0.21492612704

EL VALOR DE X2 ES=0.65159742392

EL VALOR DE X3 ES=-0.5475438818

Para finalizar con este método sería apropiado mencionar que si se detecta en la diagonal en la cual voy a trabajar que hay un cero, debemos chequear en las ecuaciones que tiene debajo y cuando se encuentre alguna en la que no es cero, intercambiar ambas ecuaciones. Si no hay ninguna que no sea diferente de cero, significa que no se puede resolver ( No es compatible determinado ).