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.