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