Ecuaciones lineales

•
SIN SIGLA

Jorge Luis Cerrada
24/5/2024
¡Este material tiene más páginas!
Vista previa del material en texto
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
 ij
 
	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
 ji
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 seguirpara 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 ).