Sorry, you're right I made a mistake while deriving that equation in python the actual matrix that I computed using mathematica is
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, -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, 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, -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, -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], [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, 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, 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, -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, -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], [-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, 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, -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, 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, 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]]) Using the one I provided in mathematica yields an error right away. Sorry about it. The _rep approach seems to work nicely and definitely something I will use a lot from now on. The answer is a little convoluted and I'm trying to get it simplified just to compare it with mathematica I'll try some of the other suggestions. Thanks a lot Oscar El martes, 11 de enero de 2022 a las 20:36:56 UTC+1, [email protected] escribió: > Would Mathematic return a Penrose inverse, if det(M) =0? I do not know > Mathematica at all. > > On Tue 11. Jan 2022 at 20:08, Oscar Benjamin <[email protected]> > wrote: > >> On Tue, 11 Jan 2022 at 12:30, Gerardo Suarez <[email protected]> >> 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 [email protected]. >> > To view this discussion on the web visit >> https://groups.google.com/d/msgid/sympy/CAHVvXxTtngh4LHPu2aid9bqycn%2B7cvg0cZ8Ey0_r5ms4CYzQig%40mail.gmail.com >> . >> > -- > Best regards, > > Peter Stahlecker > -- 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 [email protected]. To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/dae89e4b-78cc-4caa-bd4e-01b1e34fbf23n%40googlegroups.com.
