1ª Practica: Sistemas Lineales (Métodos directos)

Ejercicio 1º.
Sea el siguiente sistema:

({{2.1754` x1+4.0231` x2-2.1732` x3+5.1967` x4==17.108`}, {-4.0231` x1+6 x2+1.1973` x4==-6.1593`}, {-x1-5.2107` x2+1.1111` x3==3.0004`}, {6.0235` x1+7 x2-4.1561` x4==0}})

Resolverlo utilizando los siguientes procedimientos de Mathematica :
	- Solve
	- LinearSolve
	- Mediante las funciones preparadas de eliminación gaussiana (GaussElim []) y de pivotaje parcial escalado (PivParEsc[])
	- LUDecomposition

1º.- Comando Solve

2º.- El comando LinearSolve sirve para sistemas en que la matriz de coeficientes es tanto numérica como simbólica. La matriz puede ser cuadrada o rectangular. en el caso de sistemas indeterminados dara como resultado una de las posibles soluciones mientras que solve da la solución general.
La salida que se obtiene también es distinta.
Cuando es necesario resolver el mismo sistema con distintos vectores de términos independientes es mucho más eficiente utilizar el comandos LUDecomposition.

3º.- Se utilizan dos funciones creadas por nosotros para ilustrar el procedimiento de resolución de istemas de ecuaciones lineales mediante la eliminación gaussiana y mediante el pivotaje parcial escalado.
Se ejecutan las funciones para que mathematica las identifique y sepa como operar y a continuación se aplican a las matrices de coeficientes y termino independiente que se desee.

GaussElim [AA_,bb_]:=Module[{x,c,m,n,p,i,j,k,A,NoUniqueSoln},A=Transpose[Append[Transpose[AA],b]];n=Length[A[[All,1]]];m=A;x=A[[All,1]];
p=1  ;NoUniqueSoln=0; Print["Sistema Inicial",MatrixForm[A]];
  For
[i=1,i<n,i++,
Print["Etapa",i];
	If[A[[i,i]]≠0,p=i,
	p=i;While[(p≤n)&&(p≥i)&&(Abs[A[[p,i]]]<(10^(-15))),p=p+1]];
	If[p>n,NoUniqueSoln=1;Abort,If[p≠i,For[j=1,j≤n+1,j++,c=A[[i,j]];A[[i,j]]=A[[p,j]];A[[p,j]]=c];Print["F",p,"<->F",i,MatrixForm[A]],Continue]For[j=i+1,j≤n,j++,m[[j,i]]=A[[j,i]]/A[[i,i]]; 
A=ReplacePart[A,A[[j]]-m[[j,i]]*A[[i]] ,j];
Print["F",j,"-(",m[[j,i]],")F",i,MatrixForm[Chop[A]]]
]  ] ]  ;  
If
[(A[[n,n]]==0) ||( NoUniqueSoln==1),Abort,
		x[[n]]=A[[n,n+1]]/A[[n,n]] ;
For[i=n-1,i≥1,i--,x[[i]]=(A[[i,n+1]]-Sum[A[[i,j]]*x[[j]],{j,i+1,n}])/A[[i,i]]  ]
];
Print["Sistema Triangularizado",MatrixForm[Chop[A]]];
If[NoUniqueSoln==1 ,Print["No existe solución única"],Print["Solución=",x]]]

PivParEsc[AA_,B_]:=Module[{s,per,i,j,k,calpiv,pospiv,el,m,a,b,sol},a=AA;b=B;
s=Table[Abs[Max[Abs[a[[i]]]]],{i,1,Length[a]}];
per=Table[i,{i,1,Length[a]}];
For[i=1,i<Length[a],
calpiv=Table[Abs[a[[per[[j]],i]]/s[[per[[j]]]]],{j,i,Length[per]}];
pospiv=Flatten[Position[calpiv,Max[calpiv]]];
pospiv=pospiv[[1]]+(i-1);
el=per[[pospiv]];
per=Delete[per,pospiv];
per=Insert[per,el,i];
Print["Etapa ",i,"   Pivote ",a[[el,i]]," situado en la posición (",el,",",i,")"];
Print["Vector de permutación:",per];

For[j=i+1, j≤Length[a],m=((a[[per[[j]],i]])/(a[[el,i]]));

For[k=i+1,k≤Length[a],
a[[per[[j]],k]]=a[[per[[j]],k]]-m*a[[el,k]];k++];


a[[per[[j]],i]]=0;b[[per[[j]]]]=b[[per[[j]]]]-m*b[[el]];j++];
Print[MatrixForm[a]];
Print["Término independiente=",b];
i++];

Print["Por sustitución regresiva se obtiene: "];

sol=Table[0,{i,1,Length[a]}];
sol[[Length[a]]]=((b[[per[[Length[a]]]]])/(a[[per[[Length[a]]],Length[a]]]));
For[i=Length[a]-1,i≥1,
sol[[i]]=((b[[per[[i]]]]-(∑_(j=i+1)^(Length[a]))a[[per[[i]],j]]*sol[[j]])/(a[[per[[i]],i]]));
i--;];
Print["Solucion= ",sol];
Return[{a,b,sol}]]

4º.-  El comando LUDecomposition realiza una descomposicion de A en producto de una triangular inferior por una triangular superior con intercambio de filas. LUDecomposition, devuelve una lista tres elementos. El primero es una combinacion de las matrices triangular inferior y superior, el segundo es un vector que especifica las filas que se han usado en el pivoteo, y el tercero es una estimación del nº de condición de la matriz. El comando LUMatrices[ ] transforma el resultado de la matriz conjunta que devuelve el comando LUDecomposition, en dos matrices L y U, de modo que ya se puede realizar el proceso de resolucion del sistema en dos etapas. Se debe tener en cuenta que si el comando LUDecomposition ha realizado intercambio de filas, tambien se debera aplicar dicho intercambio a el vector b , termino independiente del sistema que se pretende resolver.

Ejercicio 2º
Resolver el siguiente sistema usando la factorización de Choleski:
6(Subscript[x, 1])+2(Subscript[x, 2])+(Subscript[x, 3])-(Subscript[x, 4])=0
2(Subscript[x, 1])+4(Subscript[x, 2])+(Subscript[x, 3])     =7
(Subscript[x, 1])+(Subscript[x, 2])+4(Subscript[x, 3])-(Subscript[x, 4])  =-1
-(Subscript[x, 1])-(Subscript[x, 3])+3(Subscript[x, 4])      =-2

La descomposición de  Cholesky se aplica a  una matriz simétrica definida positiva y da como resultado el producto de una matriz triangular superior por una triangular inferior que es la traspuesta de la anterior.