Out[1]=
Tue 16 Sep 2014 17:45:31
Inicialización Subrutinas
Solución de ecuaciones algebraicas por eliminación de Gauss 1. Regla de Cramer vs eliminación de Gauss (EG) Tanto el método MNA como el STA, dan lugar a sistemas de ecuaciones algebraicas Ax=b. Su solución por la regla de Cramer, que es un método de interés didáctico y teórico es el que más multiplicaciones y divisiones requiere: Regla de Cramer : mult div n 1
1
donde n es el orden de la matriz. Con el método de eliminación de Gauss (EG), en cambio, se necesitan Eliminación de Gauss : mult div
1
n3
3
Para el caso de una matriz de 1010 se requerirán, con la regla de Cramer, 10!=39,916,800 multiplicaciones/divisiones, y unas 333 con EG. In[8]:= Out[8]=
10 1 39 916 800 1
In[9]:=
Out[9]=
103 N 3 333.333
EG es la base de los algoritmos para solución de sistemas de ecuaciones algebraicas lineales.
2
2
Eliminación de Gauss I.nb
2. EG semi-manual In[10]:=
A
1 2 1 2 2 3 ;x 1 3 0
x1 x2 ; b x3
0 3 ; 2
La matriz m es para guardar los pivotes m Table 0, i, 1, 3, j, 1, 3 ; Ab Join A, b, 2 ; Print "Ab", MatrixForm Ab ; Ab
In[14]:=
1 2 1 0 2 2 3 3 1 3 0 2
m 2, 1
Ab 2, 1
; Print "Pivote m21", m 2, 1 ; Ab 1, 1 Ab 2 Ab 2 m 2, 1 Ab 1 ; Print "Ab", MatrixForm Ab ; Pivote m212
Ab
In[17]:=
1 2 1 0 0 2 1 3 1 3 0 2
m 3, 1
Ab 3, 1
; Print "Pivote m31", m 3, 1 ; Ab 1, 1 Ab 3 Ab 3 m 3, 1 Ab 1 ; Print "Ab", MatrixForm Ab ; Pivote m31 1 1 2 1 0 Ab 0 2 1 3 0 1 1 2
In[20]:=
m 3, 2
Ab 3, 2 ; Print "Pivote m32", m 3, 2 ; Ab 2, 2
Ab 3 Ab 3 m 3, 2 Ab 2 ; Print "Ab", MatrixForm Ab ; 1 Pivote m32 2 1 2 1 0 Ab 0 2 1 3 1 1 0 0 2 2 In[23]:=
A Take Ab, 1, 3, 1, 3 ; b Take Ab, 1, 3, 4, 4 ; Print "U", MatrixForm A , " Print "L", MatrixForm m , "
b", MatrixForm b ; b", MatrixForm b ;
Eliminación de Gauss I.nb
1 2 1 U 0 2 1 1 0 0 0 2
L
1
0 b 3
2
1 2
0 0 0 0
0 b 3
1 2
0
1 2
Substitución hacia atrás In[27]:= Out[27]=
In[28]:= Out[28]=
In[29]:= Out[29]=
s3 Solve A 3 .x b 3 1 x3 1
s2 Solve A 2 .x b 2 . s3 1 x2 1
s1 Solve A 1 .x b 1 . s2 . s3 1 x1 1
El procedimiento de arriba es la llamada descomposición LU, con el que se obtiene 0 2
0 0 1 0 0 1 2 1 1 0 0 0 0 + 0 1 0 = 2 1 0 ; U= 0 2 1 , tales que LU=A, donde A es la L= 1 12 0 1 12 1 0 0 12 0 0 1 matriz inicial. Que LU=A se confirma como sigue: 1 2
0 0 1 2 1 1 0 . 0 2 1 1 1 1 2 1 0 0 2
In[30]:=
MatrixForm
Out[30]//MatrixForm=
1 2 1 2 2 3 1 3 0
2.5 Número de condición de una matriz Link: Planet Math Matrix Condition Number El número de condición Κ se define como: In[31]:=
Κ A : Norm A Norm Inverse A
El valor mínimo de Κ es 1, y corresponde a las matrices identidad.Son las matrices menos sensibles a perturbaciones numéricas. In[32]:= Out[32]=
Κ IdentityMatrix 20 1
En el caso de la matriz del apartado 2. anterior, Κ es In[33]:=
A
1 2 1 2 2 3 ; 1 3 0
3
4
Eliminación de Gauss I.nb
In[34]:= Out[34]=
Κ A N 62.2048
Conforme crece Κ, los resultados tienden a ser menos confiables. Sin embargo, no hay problemas de estabilidad en el apartado 2 porque se resolvió con la representación exacta de números racionales. La inestabilidad está asociada a la representación imperfecta en números reales.
2.6 Precisión Precision[x] da el número efectivo de dígitos de precisión del número x In[35]:= Out[35]=
Precision 1.2345 MachinePrecision
Mathematica utiliza la precisión de la máquina en que corre In[36]:= Out[36]=
$MachinePrecision 15.9546
Se puede cambiar el número de dígitos: In[37]:=
Out[40]=
In[41]:= Out[41]=
nDg 4; $MinPrecision nDg; $MaxPrecision nDg; 1 n4 N , nDg 3 0.3333 Precision n4 4.
In[42]:=
2 m4 N , nDg; 3
In[43]:=
Precision n4 m4
Out[43]=
4. n4
In[44]:=
Precision m4
Out[44]=
4.
Si se combinan dos números de precisión nDg (4), El resultado tiene la misma precisión. Sin embargo, si uno es real con precisión mayor que nDG, el número el resultado tendrá MachinePrecision In[45]:= Out[45]=
Precision nMach N 1.333333, nDg MachinePrecision
Si n4 se suma o se multiplica con nMach, la precisión es MachinePrecision
Eliminación de Gauss I.nb
In[46]:= Out[46]=
In[47]:= Out[47]=
In[48]:=
Precision nMach n4 MachinePrecision Precision nMach n4 MachinePrecision Regresa a precisión predetermina $MaxPrecision Infinity; $MinPrecision Infinity;
3. EG semi-manual con pivoteo Resultado exacto en fracciones racionales In[50]:=
A
6 2
2
2 3
1 3
2
1 2 1
; b
2 1 ; 0
La solución exacta en racionales es In[51]:=
LinearSolve A , b MatrixForm
Out[51]//MatrixForm=
13 5 19 5
5
In[52]:=
Print "Número de condición Κ", Κ A N ; Número de condición Κ31.5021
Como el número de condición Κ>>1, puede haber inestabilidad si se resuelve en reales.
Precisión fija sin pivoteo $MinPrecision=$MaxPrecision=n para precisión aritmética fija.
5
6
Eliminación de Gauss I.nb
In[53]:=
nDg 4; $MinPrecision nDg; $MaxPrecision nDg; 6 2 A N 2
2 3
2 1 3
, nDg;
1 2 1
x
x1 2 x2 ; b N 1 , nDg; x3 0
Print "A", MatrixForm A , "\n", "b", MatrixForm b ; 6.000 2.000 2.000 A 2.000 0.6667 0.3333 1.000 2.000 1.000 b
In[59]:=
2.000 1.000 0
m21
A 2, 1 ; A 1, 1
A 2 A 2 m21 A 1 ; b 2 b 2 m21 b 1 ; m31
A 3, 1 ; A 1, 1
A 3 A 3 m31 A 1 ; A 3, 1 0; b 3 b 3 m31 b 1 ; Print"A 1 ", MatrixForm A , "\n", "b1 ", MatrixForm b ; 6.000 2.000 2.000 A1 2.168 1019 5.421 1020 0.3333 0 1.667 1.333 b1
2.000 1.667 0.3333
Eliminación de Gauss I.nb
In[67]:=
m32
A 3, 2 A 2, 2
; Print "El pivote m32", m32, " es muy grande\n"
A 3 A 3 m32 A 2 ; b 3 b 3 m32 b 2 ; Print"A 2 ", MatrixForm A , "\n", "b2 ", MatrixForm b ; El pivote m32 3.074 1019 es muy grande 6.000 2.000 A2 2.168 1019 5.421 1020 6.667
b2
0
2.000 0.3333
1.025 1019
2.000 1.667 5.124 1019
m32 1.025 1019 hace que la última ecuación en x3 tenga coeficientes grandes y poco precisos, por lo que los resultados son incorrectos por completo. La inestabilidad resulta de la pobre representación en reales de los coeficientes de dicha ecuación. In[71]:= Out[71]=
In[72]:= Out[72]=
In[73]:= Out[73]=
s3 Solve A 3 .x b 3 1 x3 5.000 6.505 1019 x1 s2 Solve A 2 .x b 2 . s3 1 x2 0
s1 Solve A 1 .x b 1 . s2 . s3 1 x1 1.333
El resultado exacto obtenido en reales y convertido al final a reales es In[74]:= Out[74]=
In[75]:=
2.6`, 3.8`, 5.` 2.6 , 3.8 , 5.
$MaxPrecision Infinity; $MinPrecision Infinity;
7
8
Eliminación de Gauss I.nb
Precisión fija con pivoteo
Eliminación de Gauss I.nb
In[77]:=
con pivoteo nDg 4; $MinPrecision nDg; $MaxPrecision nDg; 6 2 A N 2
2 3
2 1 3
, nDg;
1 2 1
x
x1 2 x2 ; b N 1 , nDg; x3 0
Print "A", MatrixForm A , "\n", "b", MatrixForm b ; m21
A 2, 1 ; A 1, 1
A 2 A 2 m21 A 1 ; A 2, 1 0; b 2 b 2 m21 b 1 ; m31
A 3, 1 ; A 1, 1
A 3 A 3 m31 A 1 ; A 3, 1 0; b 3 b 3 m31 b 1 ; Print "Antes de intercambio A", MatrixForm A ; Print "intercambio de hileras" ; ax2 A 2 ; bx2 b 2 ; ax3 A 3 ; bx3 b 3 ; A 2 ax3; b 2 bx3; A 3 ax2; b 3 bx2; MatrixForm A Print"A 1 ", MatrixForm A , "\n", "b1 ", MatrixForm b ;
m32
A 3, 2 ; A 2, 2
A 3 A 3 m32 A 2 ; b 3 b 3 m32 b 2 ; Print"A 2 ", MatrixForm A , "\n", "b2 ", MatrixForm b ;
9
10
Eliminación de Gauss I.nb
6.000 2.000 2.000 A 2.000 0.6667 0.3333 1.000 2.000 1.000 b
2.000 1.000 0 6.000
Antes de intercambio A
0 0
2.000
2.000
5.421 1020 0.3333 1.667 1.333
intercambio de hileras Out[97]//MatrixForm=
6.000 2.000 2.000 0 1.667 1.333 0 5.421 1020 0.3333 A1
6.000 2.000 2.000 0 1.667 1.333 0 5.421 1020 0.3333
2.000 b1 0.3333 1.667 A2
6.000 2.000 2.000 0 1.667 1.333 0 0 0.3333
2.000 b2 0.3333 1.667 In[103]:= Out[103]=
In[104]:= Out[104]=
In[105]:= Out[105]=
s3 Solve A 3 .x b 3 1 x3 5.000
s2 Solve A 2 .x b 2 . s3 1 x2 3.800
s1 Solve A 1 .x b 1 . s2 . s3 1 x1 2.600
Ahora la solución coincide con la exacta In[106]:= Out[106]=
In[107]:=
2.6`, 3.8`, 5.` 2.6 , 3.8 , 5.
$MaxPrecision Infinity; $MinPrecision Infinity;
8. EG matriz de Κ=927: sensibilidad por b+∆b Soluciones MachinePrecision y exacta La componente 3 de b se utiliza como parámetro
Eliminación de Gauss I.nb
In[109]:=
12
b3
; 10 0.729 0.81 0.9 1 1 1 ;x 1.331 1.21 1.1
A
x1 x2 ; b x3
0.687 0.8338 ; b3
El número de condición Κ es grande In[111]:= Out[111]=
Κ A 927.186
La solución del sistema con MachinePrecision In[112]:= Out[112]=
SolExacta Solve A.x b x1 9.33212, x2 17.0264, x3 8.52804
No difiere de la solución exacta a continuación 729 1000 In[113]:=
Out[113]=
Solve
81 100
9 10
1
1
1
1331 1000
121 100
11 10
.x
687 1000 8338 10 000
N
b3
x1 9.33212, x2 17.0264, x3 8.52804
Κ no cambia si se intercambian hileras
In[114]:=
Out[114]=
Κ
729 1000 1331 1000
81 100 121 100
9 10 11 10
1
1
1
N
927.186
Solución con precisión fija de 3 dígitos In[115]:=
In[118]:=
nDg 3; $MinPrecision nDg; $MaxPrecision nDg; A N
729 1000
81 100
9 10
1
1
1
1331 1000
121 100
11 10
, nDg; x
x1 x2 ; b N x3
A MatrixForm Out[119]//MatrixForm=
0.729 0.810 0.900 1.00 1.00 1.00 1.33 1.21 1.10 In[120]:=
m es la matriz de pivotes m Table 0, i, 1, 3, j, 1, 3 ; Ab Join A, b, 2 ; Print "Ab", MatrixForm Ab ; Print "Operación con el pivote m21:" ;
687 1000 8338 10 000
b3
, nDg;
11
12
Eliminación de Gauss I.nb
Ab
0.729 0.810 0.900 0.687 1.00 1.00 1.00 0.834 1.33 1.21 1.10 1.20
Operación con el pivote m21:
In[124]:=
m 2, 1
Ab 2, 1 ; Print "m21", m 2, 1 Ab 1, 1
m211.37 In[125]:=
Resta de m21 Ab1 a la segunda hilera Ab 2 N Ab 2 m 2, 1 Ab 1 , nDg ; Print "Ab", MatrixForm Ab ; Ab
In[127]:=
0.729 0.810 0.900 0.687 0 0.111 0.235 0.109 1.33 1.21 1.10 1.20
m 3, 1
Ab 3, 1 ; Print "m31", m 3, 1 Ab 1, 1
m311.83 In[128]:=
Ab 3 Ab 3 m 3, 1 Ab 1 ; Print "Ab", MatrixForm Ab ; Ab
In[130]:=
0.729 0.810 0.900 0.687 0 0.111 0.235 0.109 0 0.269 0.543 0.0543
m 3, 2
Ab 3, 2 ; Print "m32", m 3, 2 Ab 2, 2
Ab 3 Ab 3 m 3, 2 Ab 2 ; Print "Ab", MatrixForm Ab ; m322.42
Ab
In[133]:=
A Take Ab, 1, 3, 1, 3 ; b Take Ab, 1, 3, 4, 4 ; Print "A", MatrixForm A , " A
In[136]:= Out[136]=
In[137]:= Out[137]=
In[138]:= Out[138]=
0.729 0.810 0.900 0.687 0 0.111 0.235 0.109 0 0 0.0244 0.208
0.729 0.810 0.900 0 0.111 0.235 0 0 0.0244
b", MatrixForm b ;
0.687 b 0.109 0.208
s3 Solve A 3 .x b 3 1 x3 8.5
s2 Solve A 2 .x b 2 . s3 1 x2 17.0
s1 Solve A 1 .x b 1 . s2 . s3 1 x1 9.3
Eliminación de Gauss I.nb
In[139]:= Out[139]=
13
SolExacta x1 9.33212, x2 17.0264, x3 8.52804
b x1 x2 x3 x1 exacta x2 exacta x3 exacta 0.9 4.3 8.9 3.74 4.3 8.9 3.74 ; 1 0.24 0.25 0.346 0.24 0.25 0.346 1.2 9.3 17.0 8.5 9.33 17.02 8.52
In[140]:=
De la tabla, Ax=b es muy sensible a cambios en b, una medida de los cuales es la expresión de abajo In[141]:=
Out[141]=
Norm 9.33, 17.02, 8.52 Norm 4.3, 8.9, 3.74 Norm 0.687, 0.8338, 1.2 Norm 0.687, 0.8338, 0.9 50.965
La norma del cambio en b provoca un cambio 51 veces mayor en el vector x. EG no da lugar a pivotes grandes, por lo que la solución con pocos dígitos no difiere tanto de la exacta como en el ejemplo anterior.l In[142]:=
$MaxPrecision Infinity; $MinPrecision Infinity;
7. Manipulación de matrices In[144]:=
a
a11, a12, a21, a22; MatrixForm a
Out[144]//MatrixForm=
In[145]:= Out[145]=
a11 a12 a21 a22
a 1 a 2 a11 a21, a12 a22
Define a second matrix. In[146]:=
b
b11, b12, b21, b22; MatrixForm b
Out[146]//MatrixForm=
In[147]:=
b11 b12 b21 b22
TreeForm a , TreeForm b List
Out[147]=
List
List
a11
,
List
a12
a21
a22
List
b11
List
b12
b21
b22
14
Eliminación de Gauss I.nb
This constructs a matrix by combining the columns of the two matrices. In[148]:=
jMJoin a, b,1 MatrixForm
Out[148]//MatrixForm=
a11 a21 b11 b21 In[149]:=
a12 a22 b12 b22
TreeForm jM
Out[149]//TreeForm=
List
List
a11
List
a12
a21
List
a22
b11
List
b12
b21
b22
A matrix can also be constructed by combining the rows of these matrices. In[150]:=
jM2 Join a, b, 2 MatrixForm
Out[150]//MatrixForm=
In[151]:=
a11 a12 b11 b12 a21 a22 b21 b22
TreeForm jM2
Out[151]//TreeForm=
List
List
a11
a12
b11
List
b12
a21
a22
b21
b22