RESOLUCIÓN DE SISTEMAS LINEALES. METODOS ITERATIVOS.

Condicionamiento de un sistema lineal

Resolver el siguiente sistema lineal mediante el comando LinearSolve:
(({{10, 7, 8, 7}, {7, 5, 6, 5}, {8, 6, 10, 9}, {7, 5, 9, 10}})).(({{(Subscript[x, 1])}, {(Subscript[x, 2])}, {(Subscript[x, 3])}, {(Subscript[x, 4])}})) = (({{32}, {23}, {33}, {31}}))
Perturbando ligeramente el término independiente se obtiene:
(({{10, 7, 8, 7}, {7, 5, 6, 5}, {8, 6, 10, 9}, {7, 5, 9, 10}})).(({{(Subscript[x, 1])}, {(Subscript[x, 2])}, {(Subscript[x, 3])}, {(Subscript[x, 4])}})) = (({{32.1}, {22.9}, {33.1}, {30.9}})). Resolver este sistema utilizando el mismo comando.
Perturbando ligeramente la matriz de coeficientes se obtiene:
(({{10, 7, 8.1, 7.2}, {7.08, 5.04, 6, 5}, {8, 5.98, 9.89, 9}, {6.99, 4.99, 9, 9.98}})).(({{(Subscript[x, 1])}, {(Subscript[x, 2])}, {(Subscript[x, 3])}, {(Subscript[x, 4])}})) =  (({{32}, {23}, {33}, {31}})). Resolver utilizando el mismo comando.
¿Qué se puede decir a la vista de los tres resultados obtenidos?

A=(({{10, 7, 8, 7}, {7, 5, 6, 5}, {8, 6, 10, 9}, {7, 5, 9, 10}})); b= (({{32}, {23}, {33}, {31}}));

Normas vectoriales y normas matriciales. Número de condición de una matriz.

Cuando se construye una sucesión de vectores aproximación de la solución de un sistema lineal, para estudiar la convergencia interesa medir la distancia entre las dos últimas aproximaciones, y para medir esta distancia es necesario utilizar el concepto de norma de un vector . Tanto en el estudio del condicionamiento de un sistema de ecuaciones lineales cuando se resuelve mediante métodos directos, como en el estudio de la convergencia de los métodos iterativos para la resolución de sistemas lineales, es necesario el concepto de norma de una matriz. El comando Norm de Mathematica nos da la norma tanto en el caso de un vector como de una matriz.

Por otra parte, algunos sistemas son muy sensibles a los errores de redondeo y el vector solución puede ser bastante inexacto. En este caso se dice que el sistema es inestable o que está mal condicionado. En este tipo de sistemas lo que suele ocurrir es que pequeños cambios en los coeficientes o en los términos independientes dan lugar a cambios apreciables en la solución. El número de condición de una matriz A se define como :
	K(A) = ||A|| . ||(A^(-1))||
y sirve para tener una medida del condicionamiento de una matriz.

Ejemplo 1º.-
Sea el vector b = (({{-1.63`}, {7.23`}, {-3.38`}})), y la matriz  : A = (({{3.02`, -1.05`, 2.53`}, {4.33`, 0.56`, -1.78`}, {-0.83`, -0.54`, 1.47`}})). Hallar las tres normas vectoriales de b: (Subscript[‖b‖, 1]),  (Subscript[‖b‖, 2]) y  (Subscript[‖b‖, ∞]). Idem para la matriz A y para su inversa. Calcular el número de condición de la matriz utilizando la norma del máximo.

A = {{3.02,-1.05,2.53}, {4.33,0.56,-1.78},{-0.83,-0.54,1.47}};
b={-1.63,7.23,-3.38};

Las tres normas del vector b:

euclideavect=Norm[b]

maxvect=Norm[b,∞]

sumavect=Norm[b,1]

Las tres normas de la matriz A

euclideamat=Norm[A]

Comprobamos que es igual a (Sqrt[ρ(A.(A^t))])

Max[Eigenvalues[A.Transpose[A]]]^.5

maxmat=Norm[A,∞]

Comprobamos que es igual a la mayor de las sumas de los elementos de las filas en valor absoluto

Max[{(∑_(j=1)^3)Abs[A[[1,j]]],(∑_(j=1)^3)Abs[A[[2,j]]],(∑_(j=1)^3)Abs[A[[3,j]]]}]

sumamat=Norm[A,1]

Comprobamos que es igual a la mayor de las sumas de los elementos de las columnas en valor absoluto

Max[{(∑_(i=1)^3)Abs[A[[i,1]]],(∑_(i=1)^3)Abs[A[[i,2]]],(∑_(i=1)^3)Abs[A[[i,3]]]}]

Repetimos el cálculo de la norma para la inversa de A

euclideainv=Norm[Inverse[A]]

maxinv=Norm[Inverse[A],∞]

sumainv=Norm[Inverse[A],1]

m=maxmat*maxinv

Ejemplo 2º.-
Sea el sistema A. x = b, siendo A y b la matriz y vector respectivamente del ejemplo anterior. Sea x = (({{1}, {2}, {-1}})) la solución exacta y (OverTilde[x]) = (({{0.88`}, {-2.34`}, {-2.66`}}))la solución aproximada. Calcular el error relativo y el residuo relativo. Comprobar que el error relativo está acotado superiormente por el producto del residuo relativo por el nº de condición de la matriz, e inferiormente por el cociente del residuo relativo entre el nº de condición de la matriz.

xe={1.,2.,-1.};xa={.88,-2.34,-2.66};e=xe-xa;errorel=Norm[e]/Norm[xe]

resi=b-A.xa;resirelat=Norm[resi]/Norm[b]

errorel
m*resirelat

errorel
resirelat/m

Hallar el Número de Condición de la matriz :
(({{10, 7, 8, 7}, {7, 5, 6, 5}, {8, 6, 10, 9}, {7, 5, 9, 10}}))

El método de Jacobi

Sea el sistema :
	(({{(Subscript[a, 11]), (Subscript[a, 12]), .., (Subscript[a, 1n])}, {(Subscript[a, 21]), (Subscript[a, 22]), .., (Subscript[a, 2n])}, {.., .., .., ..}, {(Subscript[a, n1]), (Subscript[a, n2]), .., (Subscript[a, nn])}})).(({{(Subscript[x, 1])}, {(Subscript[x, 2])}, {..}, {(Subscript[x, n])}}))= (({{(Subscript[b, 1])}, {(Subscript[b, 2])}, {..}, {(Subscript[b, n])}}))
Para resolverlo mediante el método de Jacobi se implementa el siguiente proceso iterativo:
	(Subscript[x, i]^((k)))=(1/(Subscript[a, ii])).((Subscript[b, i])-(∑_(j=1)^(i-1))(Subscript[a, ij]) (Subscript[x, j]^((k-1)))-(∑_(k=i+1)^n)(Subscript[a, ij]) (Subscript[x, j]^((k-1))))
El superíndice (k) correponde a la etapa de iteración.

Ejemplo 3º. Consideremos el sistema de ecuaciones : A . x = b, siendo A = (({{4, -1, 1}, {4, -8, 1}, {-2, -1, 5}})), y b = (({{7}, {-21}, {15}})). Realizar 5 iteraciones del mdo. de Jacobi partiendo de la aproximación inicial : x0 = (({{1}, {2}, {2}})).Calcular, en cada iteración, el error absoluto y relativo.

A={{4,-1,1},{4,-8,1},{-2,-1,5}};

b={7,-21,15};x0={1.,2,2};xn=x0;

Do[Do[xn[[i]]=(b[[i]]-(∑_(k=1)^(i-1))A[[i,k]]*x0[[k]]-(∑_(k=i+1)^3)A[[i,k]]*x0[[k]])/A[[i,i]],{i,3}];Print[xn];erabs=Norm[xn-x0];errel=erabs/Norm[xn];Print["El error absoluto es: ",erabs];Print["El error relativo es: ",errel];x0=xn,{iter,5}]

El método de Gauss-Seidel

Considerando el mismo sistema anterior:
	(({{(Subscript[a, 11]), (Subscript[a, 12]), .., (Subscript[a, 1n])}, {(Subscript[a, 21]), (Subscript[a, 22]), .., (Subscript[a, 2n])}, {.., .., .., ..}, {(Subscript[a, n1]), (Subscript[a, n2]), .., (Subscript[a, nn])}})).(({{(Subscript[x, 1])}, {(Subscript[x, 2])}, {..}, {(Subscript[x, n])}}))= (({{(Subscript[b, 1])}, {(Subscript[b, 2])}, {..}, {(Subscript[b, n])}}))
Para resolverlo mediante el método de Gauss-Seidel se implementa el siguiente proceso iterativo:
	(Subscript[x, i]^((k)))=(1/(Subscript[a, ii])).((Subscript[b, i])-(∑_(j=1)^(i-1))(Subscript[a, ij]) (Subscript[x, j]^((k)))-(∑_(k=i+1)^n)(Subscript[a, ij]) (Subscript[x, j]^((k-1))))
El superíndice (k) correponde a la etapa de iteración.

Ejemplo 4º. Consideremos el sistema de ecuaciones del ejemplo anterior . Realizar 5 iteraciones del mdo. de Gauss-Seidel partiendo de la misma aproximación inicial .Calcular el error absoluto y relativo.

A={{4,-1,1},{4,-8,1.},{-2,-1,5}};

b={7,-21,15};x0={1.,2,2};xn=x0;

Do[Do[xn[[i]]=(b[[i]]-(∑_(k=1)^(i-1))A[[i,k]]*xn[[k]]-(∑_(k=i+1)^3)A[[i,k]]*x0[[k]])/A[[i,i]],{i,3}];Print[xn];erabs=Norm[xn-x0];errel=erabs/Norm[xn];Print["El error absoluto es: ",erabs];Print["El error relativo es: ",errel];x0=xn,{iter,5}]

La condición necesaria y suficiente para que un método iterativo de resolución de sistemas lineales sea convergente es que el radio espectral de la matriz de paso sea menor que la unidad. (Se recuerda que el radio espectral de una matriz es igual al mayor de los valores popios de la matriz en módulo o valor absoluto) con Mathematica resulta muy sencillo calcular el radio espectral de una matriz utilizando el comando Eigenvalues.
Sea el sistema :
	(({{(Subscript[a, 11]), (Subscript[a, 12]), .., (Subscript[a, 1n])}, {(Subscript[a, 21]), (Subscript[a, 22]), .., (Subscript[a, 2n])}, {.., .., .., ..}, {(Subscript[a, n1]), (Subscript[a, n2]), .., (Subscript[a, nn])}})).(({{(Subscript[x, 1])}, {(Subscript[x, 2])}, {..}, {(Subscript[x, n])}}))= (({{(Subscript[b, 1])}, {(Subscript[b, 2])}, {..}, {(Subscript[b, n])}})), siendo A = (({{(Subscript[a, 11]), (Subscript[a, 12]), .., (Subscript[a, 1n])}, {(Subscript[a, 21]), (Subscript[a, 22]), .., (Subscript[a, 2n])}, {.., .., .., ..}, {(Subscript[a, n1]), (Subscript[a, n2]), .., (Subscript[a, nn])}})).

Definamos las matrices :
	D= (({{(Subscript[a, 11]), 0, .., 0}, {0, (Subscript[a, 22]), .., 0}, {.., .., .., ..}, {0, 0, .., (Subscript[a, nn])}})),  	L = (({{0, 0, .., 0}, {-(Subscript[a, 21]), 0, .., 0}, {.., .., .., ..}, {-(Subscript[a, n1]), -(Subscript[a, n2]), .., 0}})) 	y 
	
	U = (({{0, -(Subscript[a, 12]), .., -(Subscript[a, 1n])}, {0, 0, .., -(Subscript[a, 2n])}, {.., .., .., ..}, {0, 0, .., 0}})).
Entonces la matriz de paso del método de Jacobi es :
	(Subscript[T, j ]) = (D^(-1)).(L+U)
y la del método de Gauss-Seidel :
	(Subscript[T, G]) = (((D - L))^(-1)).U

Ejemplo 5º.Calcular las matrices de paso de los métodos de Jacobi y Gauss-Seidel de los ejemplos anteriores,y sus respectivos radios espectrales,comprobando que son menores que la unidad.

d=Table[A[[i,i]],{i,3}];diag=DiagonalMatrix[d]

L=-A;Do[If[j≥i,L[[i,j]]=0],{i,3},{j,3}];MatrixForm[L]

U=-A;Do[If[j<=i,U[[i,j]]=0],{i,3},{j,3}];MatrixForm[U]

Para el método de Jacobi T = (D^(-1)).(L+U):

Tj=Inverse[diag].(L+U);MatrixForm[Tj]

Max[Abs[Eigenvalues[Tj]]]

Para el método de Gauss-Seidel T = (((D-L))^(-1)).U:

Tg=Inverse[diag-L].U;MatrixForm[Tg]

Max[Abs[Eigenvalues[Tg]]]

Resolución de sistemas no lineales. El método de Newton. El comando FindRoot.

Sea el sistema {({{f1(x,y) = 0}, {f2(x,y) = 0}}) , y sea J(x,y) la matriz Jacobiana de f1 y f2 respecto a x e y. Sea (Subscript[P, k]) ((Subscript[x, k]), (Subscript[y, k])) la aproximación de la solución correspondiente a la etapa k. Los pasos correspondientes a la siguiente etapa son:
1º.- Evaluar F((Subscript[P, k])) = (({{f1((Subscript[P, k]))}, {f2((Subscript[P, k]))}})).
2º.- Calcular la matriz jacobiana en (Subscript[P, k]), J((Subscript[P, k])).
3º.- Resolver el sistema lineal :  J((Subscript[P, k])).Δp = - F( J((Subscript[P, k])).)
4º.- Calcular la siguiente aproximación : (Subscript[P, k+1]) = (Subscript[P, k]) + Δp .

Ejemplo 6º : Sea el sistema : 
	{({{(x^2) - y = 0.2}, {(y^2) - x = 0.3}}) , 
Representar gráficamente para obtener la aproximación inicial correspondiente a la raiz del primer cuadrante, y calcular 3 iteraciones del método de Newton .

f[x_,y_]=x^2-y-.2;g[x_,y_]=y^2-x-.3;

graf1=ContourPlot[f[x,y]==0,{x,-2,2},{y,-2,2}]

graf2=ContourPlot[g[x,y]==0,{x,-2,2},{y,-2,2}]

Show[graf1,graf2]

jacob[x_,y_]={{(∂_x) f[x,y],(∂_y) f[x,y]},{(∂_x) g[x,y],(∂_y) g[x,y]}};indep[x_,y_]={-f[x,y],-g[x,y]};

Print[MatrixForm[jacob[x,y]],"    ", MatrixForm[indep[x,y]]]

x0=1.2;y0=1.2;itermax=3;

Do[Print["Jacob= ",MatrixForm[jacob[x0,y0]],"  ","Term.In= ",MatrixForm[indep[x0,y0]]];sol=LinearSolve[jacob[x0,y0],indep[x0,y0]];Print["Delta= ",MatrixForm[sol]];dx=sol[[1]];dy=sol[[2]];x0=x0+dx;y0=y0+dy;Print["K= ",k," x= ",x0,"  y= ",y0],{k,itermax}]

El comando FindRoot

Como en el caso de una sola ec. no lineal, el comando FindRoot nos permite obtener una solución numérica (aproximada) de un sistema de ecuaciones no lineales partiendo de una aproximación inicial. Una de sus ociones nos permite especificar el nº máximo de iteraciones a realizar. Si se alcanza el nº máximo de iteraciones sin haber obtenido la solución el comando FindRoot nos devuelve la aproximación calculada en la última etapa que se puede utilizar como aproximación inicial de un nuevo comando FindRoot.
Es conveniente realizar una representación gráfica, cuando sea posible, para obtener una buena aproximación inicial.
El comando FindRoot puede resolver ec. no lineales en variable compleja. Para ello solo es preciso incluir la unidad imaginaria I en la aproximación inicial. Por ejemplo:

FindRoot[Sin[x] == 2, {x, 1}]

FindRoot[Sin[x] == 2, {x, I}]

Para resolver el sistema anterior utilizando el comando FindRoot:

FindRoot[{f[x,y]==0,g[x,y]==0},{x,1.2},{y,1.2}]

Ejercicios

1º.- El  sistema lineal :
0.04 x1 + 0.01 x2 - 0.01 x3 = 0.06
0.20 x1 + 0.50 x2 - 0.20 x3 = 0.30
       x1 +    2  x2 +    4  x3 =  11
tiene como solución exacta xe = {1.827586,0.6551724,1.965517} y como solución aproximada xa = {1.8,0.64, 1.9}. Calcular , para las normas (Subscript[l, 1]), (Subscript[l, 2]) y (Subscript[l, ∞]) : 
	a) la norma del vector error absoluto y del vector residuo. 
	b) el error relativo y el residuo relativo.
	c) la norma  y el número de condición de la matriz de coeficientes, comprobando que cualquiera de las normas es mayor que el radio espectral de la matriz.

mat={{.04,.01,-.01},{.2,.5,-.2},{1,2,4}};b={.06,.3,11};xe={1.827586,0.6551724,1.965517};xa={1.8,.64,1.9};

a)

b)

c)

2º.- Resolver el sistema : 
(({{-1, 11, -1, 3}, {2, -1, 10, -1}, {10, -1, 2, 0}, {0, 3, -1, 8}})). (({{x}, {y}, {z}, {t}}))= (({{25}, {-11}, {6}, {15}})) , partiendo de la aproximación inicial x0=(({{0}, {0}, {0}, {0}}))y
utilizando el método de Jacobi y con un error absoluto menor que (10^(-3)).

Reordenando el sistema para que la matriz de los coeficientes sea estrictamente diagonal dominante: (({{10.`, -1, 2, 0}, {-1, 11, -1, 3}, {2, -1, 10, -1}, {0, 3, -1, 8}})).(({{x}, {y}, {z}, {t}}))=(({{6}, {25}, {-11}, {15}})), y, partiendo de la aproximación inicial x0=(({{0}, {0}, {0}, {0}}))con el método de Jacobi y con un error absoluto menor que (10^(-3)).

3º.-Resolver el siguiente sistema :
(({{-2, 1, 0, 0}, {1, -2, 1, 0}, {0, 1, -2, 1}, {0, 0, 1, -2}})).(({{x}, {y}, {z}, {t}})) = (({{1}, {0}, {0}, {0}}))
mediante el método de Gauss-Seidel partiendo de la aproximación inicial : (({{0.7`}, {0.5`}, {0.3`}, {0.1`}})). Comprobar la convergencia del método.

4º- Resolver el siguiente sistema mediante los métodos de Jacobi y Gauss-Seidel para n=10, 50 y 100. Obtener las soluciones con un error menor que (10^(-3)) utilizando la norma del máximo. Comparar el nº de iteraciones requeridos por cada método.

(({{1, -(1/2), 0, 0, ., ., 0}, {-(1/2), 1, -(1/2), 0, ., ., 0}, {0, -(1/2), 1, -(1/2), 0, ., 0}, {., ., ., ., ., ., .}, {., ., ., ., ., ., .}, {0, ., ., 0, -(1/2), 1, -(1/2)}, {0, ., ., ., 0, -(1/2), 1}})).(({{(Subscript[x, 1])}, {(Subscript[x, 2])}, {(Subscript[x, 3])}, {.}, {.}, {.}, {.}, {(Subscript[x, n-1])}, {(Subscript[x, n])}}))= (({{(1/2)}, {0}, {0}, {.}, {.}, {.}, {.}, {0}, {0}}))

a) Jacobi

b) Gauss-Seidel