Evidentemente se ocasiona el problema cuando se usan numeros de tipo flotante en las entradas de la matríz. No pude entender más que eso y no pude resolver el problema ingresando los números flotantes. El único remedio que encontré es introducir los números de tipo racional. El tema es el siguiente, si uno pone 2/100, para python eso es cero porque el dice "tengo un conciente de dos enteros devuelvo un entero que represente el cociente", cómo el cociente no es entero, devuelve el número redondeado. Para introducir el racionales hay que usar la función Rational de sympy.
>> Rational(2,100)
Da el
$\frac{1}{2}$
Les digo que es más práctico multiplicar por 100 la matríz, como hacían uds.
In [1]:
import sympy as sym
sym.init_printing()
Como es poco práctico escribir "Rational(x,y)" cada vez que quiero escribir una entrada de la matríz. Defino una función (R), que es esencialmente Rational pero con un nombre más breve. Uso el comando lambda que permite definir funciones rapidamente.
In [2]:
R=lambda n,d: sym.Rational(n,d)
R(1,2)
Out[2]:
In [3]:
A=sym.Matrix([[R(2,100), 0,R(119,10),0], [R(5,100),R(12,100),R(57,10),0],\
[0,R(14,100),R(21,100),R(32,100)],[0,0,R(43,100),R(11,100)]])
A
Out[3]:
In [4]:
lambd=sym.symbols('lambda')
p=sym.det(A-lambd*sym.eye(4))
p
Out[4]:
In [5]:
sym.roots(p)
Out[5]:
La precisión infinita de sympy lo lleva a escribir expresiones monstruosas (pero exactas). Vamos a usar el comando N de sympy que convierte las expresiones en complejos-flotantes. El comando N se aplica a números y no a diccionarios (resultado de roots(p)), por eso lo tengo que aplicar uno por uno a las claves del diccionario, pero se hace facil.
In [6]:
autovalores=[sym.N(a) for a in sym.roots(p)]
autovalores
Out[6]:
Los autovalores, si bien complejos, tienen una parte imaginaria ridiculamente chica, de un orden inferior a los errores de redondeo por aritmética de punto flotante. Desestimemos la parte imaginaria. El comando re de sympy calcula la parte real de un complejo.
In [29]:
AutSym=[ sym.re(a) for a in autovalores]
AutSym
Out[29]:
In [8]:
Aut=list(sym.roots(p))
Aut
Out[8]:
In [9]:
#sistema=A-Aut[0]*sym.eye(4),sym.zeros(4,1)
#x=sym.symbols('x:4')
#sym.linsolve(sistema,x)
In [11]:
import numpy as np
from numpy import linalg as LA
In [12]:
M=np.array([[.02, 0, 11.9, 0], [.05, .12, 5.7, 0], [0, .14, .21, .32], [0, 0, .43, .11]])
In [23]:
P=LA.eig(M)[1]
In [26]:
sym.Matrix(P)
Out[26]:
In [36]:
v=sym.Matrix([1,1,1,1])
Vect=sym.N((A-AutSym[0]*sym.eye(4))**(-100)*v/(max((A-AutSym[0]*sym.eye(4))**(-100)*v)))
Vect
Out[36]:
In [37]:
0.398998983743746*Vect
Out[37]:
In [ ]: