sympy: simplify a larger expression with Binomial formula and quadratic complement

Question

I get some eigenvalues of the matrix

import sys
import mpmath
from sympy import *

X,Y,Z = symbols("X,Y,Z")
Rxy,Rxz, Ry,Ryx,Ryz, Rz,Rzy,Rzz = symbols("Rxy,Rxz,  Ry,Ryx,Ryz, Rz,Rzy,Rzz")


J = Matrix([
       [      -1,        0,        0],
       [       0,    -Ry*Y, Ry*Rzy*Y],
       [Rxz*Rz*Z, Ryz*Rz*Z,    -Rz*Z]])

which are the following:

{-Ry*Y/2 - Rz*Z/2 + sqrt(Ry**2*Y**2 + 4*Ry*Ryz*Rz*Rzy*Y*Z - 2*Ry*Rz*Y*Z + Rz**2*Z**2)/2: 1,
 -Ry*Y/2 - Rz*Z/2 - sqrt(Ry**2*Y**2 + 4*Ry*Ryz*Rz*Rzy*Y*Z - 2*Ry*Rz*Y*Z + Rz**2*Z**2)/2: 1,
 -1: 1}

lets just look at eigenvalue one:

In [25]: J.eigenvals().keys()[0]
Out[25]: -Ry*Y/2 - Rz*Z/2 + sqrt(Ry**2*Y**2 + 4*Ry*Ryz*Rz*Rzy*Y*Z - 2*Ry*Rz*Y*Z + Rz**2*Z**2)/2

I want to simplify this term the following: factor out 1/2 and (this is important) the radicant.

I can transform the radicant as following by adding the quadratic complement

Ry**2*Y**2 + 4*Ry*Ryz*Rz*Rzy*Y*Z - 2*Ry*Rz*Y*Z + Rz**2*Z**2     | + 4*Ry*Rz*Y*Z -4*Ry*Rz*Y*Z

which leads to

Ry**2*Y**2 + Rz**2*Z**2 + 2*Ry*Rz*Y*Z - 4*Ry*Rz*Y*Z + 4*Ry*Ryz*Rz*Rzy*Y*Z

which can be factorized to

(Ry*Y + Rz*Z)**2 - 4*Ry*Rz*Y*Z*(1 - Ryz*Rzy)

with these evaluations the complete eigenvalue should look like this

-1/2*(Ry*Y + Rz*Z - sqrt((Ry*Y + Rz*Z)**2 - 4*Ry*Rz*Y*Z*(1 - Ryz*Rzy)))

This calculation is really important for me because I have to evaluate if the eigenvalue is <0. Which is much easier in the last form.

Let me show you what I did until now.

In [24]: J.eigenvals().keys()[0]
Out[24]: -Ry*Y/2 - Rz*Z/2 + sqrt(Ry**2*Y**2 + 4*Ry*Ryz*Rz*Rzy*Y*Z - 2*Ry*Rz*Y*Z + Rz**2*Z**2)/2

In [25]: J.eigenvals().keys()[0].factor() 
Out[25]: -(Ry*Y + Rz*Z - sqrt(Ry**2*Y**2 + 4*Ry*Ryz*Rz*Rzy*Y*Z - 2*Ry*Rz*Y*Z + Rz**2*Z**2))/2

In [26]: J.eigenvals().keys()[0].simplify()
Out[26]: -Ry*Y/2 - Rz*Z/2 + sqrt(Ry**2*Y**2 + 4*Ry*Ryz*Rz*Rzy*Y*Z - 2*Ry*Rz*Y*Z + Rz**2*Z**2)/2

So .simplify() doesnt changes the result at all. .factor() just factors out the -1/2. If I remember correct, I can pass an argument to .factor() like Y or Z, which variable should be factorized. But I get a lot of slightly different eigenvalues as output and I don't want to specify each argument of factor() by hand (if this solution even works).

I also tried to calculate the eigenvalues by myself by calculating the determinant and solve determinat==0... I also used determinat.factor() and solved it afterwards but the best result of this approach was the same as J.eigenvals().keys()[0].factor().

Do you have any idea how to solve this problem?

Thank you in advance

Alex

Answer 1

This sort of thing is asked for a lot (see also, for instance, this question: Expression simplification in SymPy), but there isn't really a good way in SymPy to do it. The problem is that such "partial" factorizations are not unique (there may be multiple ways to convert a polynomial into a sum of products).

I opened this issue about it in the SymPy issue tracker. I showed there a way that you can get close (here a is the term under the square root)

In [92]: collect(expand(a.subs(Ry*Y, x - Rz*Z)), x, func=factor).subs(x, Ry*Y + Rz*Z)
Out[92]:
      2  2                                                                   2
- 4⋅Rz ⋅Z ⋅(Ryz⋅Rzy - 1) + 4⋅Rz⋅Z⋅(Ry⋅Y + Rz⋅Z)⋅(Ryz⋅Rzy - 1) + (Ry⋅Y + Rz⋅Z)

Here I am temporarily replacing Ry*Y + Rz*Z with a variable x so that I can get the squared term that you want.

I couldn't figure out a way to get closer to exactly what you want (i.e., factor the Ryz*Rzy - 1 out of the remaining terms).