On Tue, 11 Jan 2022 at 12:30, Gerardo Suarez <gerardofingu...@gmail.com> wrote:
>
> Sure, sorry about the image, I was unaware of how easy is to generate python
> code with the print methods
>
> from sympy import Symbol, Matrix,I
> gamma_c = Symbol('gamma_c')
> n_c = Symbol('n_c')
> gamma_h = Symbol('gamma_h')
> n_h = Symbol('n_h')
> epsilon_c = Symbol('epsilon_c')
> g = Symbol('g')
> epsilon_h = Symbol('epsilon_h')
> e = Matrix([[-gamma_c*(n_c + 1) - gamma_h*(n_h + 1), 0, 0, 0, 0, gamma_c*n_c,
> 0, 0, 0, 0, gamma_h*n_h, 0, 0, 0, 0, 0], [0, -I*epsilon_c - gamma_c*n_c/2 -
> gamma_c*(n_c + 1)/2 - gamma_h*(n_h + 1), I*g, 0, 0, 0, 0, 0, 0, 0, 0,
> gamma_h*n_h, 0, 0, 0, 0], [0, I*g, -I*epsilon_h - gamma_c*(n_c + 1) -
> gamma_h*n_h/2 - gamma_h*(n_h + 1)/2, 0, 0, 0, 0, gamma_c*n_c, 0, 0, 0, 0, 0,
> 0, 0, 0], [0, 0, 0, -gamma_c*n_c/2 - gamma_c*(n_c + 1)/2 - gamma_h*n_h/2 -
> gamma_h*(n_h + 1)/2 - I*(epsilon_c + epsilon_h), 0, 0, 0, 0, 0, 0, 0, 0, 0,
> 0, 0, 0], [0, 0, 0, 0, -gamma_c*n_c/2 - gamma_c*(n_c + 1)/2 - gamma_h*(n_h +
> 1) - I*(epsilon_h - (epsilon_c + epsilon_h)), 0, 0, 0, -I*g, 0, 0, 0, 0, 0,
> gamma_h*n_h, 0], [gamma_c*(n_c + 1) - gamma_h*n_h, 0, 0, 0, 0, -gamma_c*n_c -
> gamma_h*n_h - gamma_h*(n_h + 1), I*g, 0, 0, -I*g, -gamma_h*n_h, 0, 0, 0, 0,
> 0], [0, 0, 0, 0, 0, I*g, -gamma_c*n_c/2 - gamma_c*(n_c + 1)/2 - gamma_h*n_h/2
> - gamma_h*(n_h + 1)/2 - I*(-epsilon_c + epsilon_h), 0, 0, 0, -I*g, 0, 0, 0,
> 0, 0], [0, 0, gamma_c*(n_c + 1), 0, 0, 0, 0, -I*epsilon_h - gamma_c*n_c -
> gamma_h*n_h/2 - gamma_h*(n_h + 1)/2, 0, 0, 0, -I*g, 0, 0, 0, 0], [0, 0, 0, 0,
> -I*g, 0, 0, 0, -gamma_c*(n_c + 1) - gamma_h*n_h/2 - gamma_h*(n_h + 1)/2 -
> I*(epsilon_c - (epsilon_c + epsilon_h)), 0, 0, 0, 0, gamma_c*n_c, 0, 0], [0,
> 0, 0, 0, 0, -I*g, 0, 0, 0, -gamma_c*n_c/2 - gamma_c*(n_c + 1)/2 -
> gamma_h*n_h/2 - gamma_h*(n_h + 1)/2 - I*(epsilon_c - epsilon_h), I*g, 0, 0,
> 0, 0, 0], [-gamma_c*n_c + gamma_h*(n_h + 1), 0, 0, 0, 0, -gamma_c*n_c, -I*g,
> 0, 0, I*g, -gamma_c*n_c - gamma_c*(n_c + 1) - gamma_h*n_h, 0, 0, 0, 0, 0],
> [0, gamma_h*(n_h + 1), 0, 0, 0, 0, 0, -I*g, 0, 0, 0, -I*epsilon_c -
> gamma_c*n_c/2 - gamma_c*(n_c + 1)/2 - gamma_h*n_h, 0, 0, 0, 0], [0, 0, 0, 0,
> 0, 0, 0, 0, 0, 0, 0, 0, -gamma_c*n_c/2 - gamma_c*(n_c + 1)/2 - gamma_h*n_h/2
> - gamma_h*(n_h + 1)/2 + I*(epsilon_c + epsilon_h), 0, 0, 0], [0, 0, 0, 0, 0,
> 0, 0, 0, gamma_c*(n_c + 1), 0, 0, 0, 0, I*epsilon_h - gamma_c*n_c -
> gamma_h*n_h/2 - gamma_h*(n_h + 1)/2, I*g, 0], [0, 0, 0, 0, gamma_h*(n_h + 1),
> 0, 0, 0, 0, 0, 0, 0, 0, I*g, I*epsilon_c - gamma_c*n_c/2 - gamma_c*(n_c +
> 1)/2 - gamma_h*n_h, 0], [gamma_c*n_c + gamma_h*n_h, 0, 0, 0, 0, gamma_c*n_c +
> gamma_h*n_h + gamma_h*(n_h + 1), 0, 0, 0, 0, gamma_c*n_c + gamma_c*(n_c + 1)
> + gamma_h*n_h, 0, 0, 0, 0, 0]])
>
> However, I care more about how to speed up matrix inversion in general, as it
> was terribly fast in Mathematica and it never finished in sympy. I guess
> there is a more efficient way than using
>
> e.inv()
The determinant of this matrix is zero so it isn't invertible. I
checked this using DomainMatrix:
In [87]: from sympy.polys.matrices import DomainMatrix
In [88]: dM = DomainMatrix.from_Matrix(e)
In [89]: %time ok = dM.charpoly()
CPU times: user 2min 36s, sys: 100 ms, total: 2min 36s
Wall time: 2min 36s
In [90]: ok[-1] # constant term of charpoly is the determinant
Out[90]: (0 + 0*I)
You can check more quickly just by substituting random values for the symbols:
In [91]: vals = dict(zip(e.free_symbols, randMatrix(7, 1)))
In [92]: e.subs(vals).det()
Out[92]: 0
In [93]: vals
Out[93]: {ε_c: 6, εₕ: 92, g: 44, γ_c: 84, γₕ: 63, n_c: 37, nₕ: 42}
Are you sure that this matrix is the same as the one that you used in
Mathematica?
Did Mathematica actually return an inverse or did it just say that the
matrix is not invertible?
Also I can see now how inv can be made much faster for cases like
this. The matrix has a number of different components:
In [101]: e.connected_components()
Out[101]: [[0, 5, 6, 10, 9, 15], [1, 2, 7, 11], [3], [4, 8, 13, 14], [12]]
The inverse can be computed from the submatrices generated by the
numbered rows and columns. The determinant is the product of their
determinants and the first submatrix has a zero determinant:
In [110]: ccs = e.connected_components()
In [111]: e[ccs[0],ccs[0]].det()
Out[111]: 0
In fact this is zero because its right column is zero which in turn is
because column 15 of the matrix is zero:
In [121]: list(e[:,15])
Out[121]: [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
Not sure why I didn't notice that straight away!
In [127]: e.print_nonzero()
[X X X ]
[ XX X ]
[ XX X ]
[ X ]
[ X X X ]
[X XX XX ]
[ XX X ]
[ X X X ]
[ X X X ]
[ X XX ]
[X XX XX ]
[ X X X ]
[ X ]
[ X XX ]
[ X XX ]
[X X X ]
That's definitely a case that could be handled more efficiently by the
various routines but it makes me wonder if you showed the right
matrix.
--
Oscar
--
You received this message because you are subscribed to the Google Groups
"sympy" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to sympy+unsubscr...@googlegroups.com.
To view this discussion on the web visit
https://groups.google.com/d/msgid/sympy/CAHVvXxTtngh4LHPu2aid9bqycn%2B7cvg0cZ8Ey0_r5ms4CYzQig%40mail.gmail.com.
