symbolic complex expression simplification

Although my question is related to a specific problem, I would like to approach it in more general terms. I would like to simplify a fractional complex expression obtained by multiplying symbolic matrices using the sympy package. What I get is a fraction with real parameters and many complex exponential terms (phase terms) like exp(-jd), exp(-2jd) and also exp(-4j*d). I get the correct result, but when I try to calculate the ||**2, which is a real expression, sympy.simplify() is not able to manage the phase terms and I obtain a huge expression I have to reduce by hand. My test procedure, being T, M, M_inv, F and T, 2×2 symbolic matrices is:

import simply as sym
Mf = T @ M_inv @ F @ T
Mf = sym.simplify(Mf)
eig = Mf.eigenvals() 
for key,value in eig.items():
    val = key
eig_cc = np.conj(val)
final = sym.simplify(val*eig_cc)

I wonder if a specific tool exists to improve the calculation.

After a fast research I found two possibilities in the post:

Exponential to trigonometric conversion in SymPy while simplifying – a stubborn expression

1. a.trigsimp() for a in expr.rewrite(cos).expand().as_real_imag()

2. expr.rewrite(cos).simplify().trigsimp()

While the first provides the correct result for easier expressions (not for the one I am trying to solve), none of them works in my case.

Second question, eig output has a dict shape, so to take the value I used the for loop.

for key,value in eig.items():
    val = key

Is there a better, faster, way to achieve the same result?

The code I am using is:

import numpy as np
import sympy as sym
from sympy import *
t = sym.Symbol('t', real=True)
r = sym.Symbol('r', real=True)
u = sym.Symbol('u', real=True)
d = sym.Symbol('d', real=True)
a = sym.Symbol('a', real=True)


I = np.identity(2)

T = sym.Matrix ([[t,0],[0,t]])
R = sym.Matrix ([[r,0],[0,r]])
F = sym.Matrix ([[sym.sqrt(u)*sym.exp(-1j*d),0],[0,sym.sqrt(u)*sym.exp(-1j*d)]])



sym.init_printing()
M = (I - F @ R @ F @ R)
M_inv = M.inv() 

Mf = T @ M_inv @ F @ T

Mf = sym.simplify(Mf)

eig = Mf.eigenvals()



for key,value in eig.items():
    val = key

eig_cc = np.conj(val)
expr = finale = sym.simplify(val*eig_cc)

"""
Test1
"""

result = [a.trigsimp() for a in expr.rewrite(sin).expand().as_real_imag()]

"""
Test2
"""

result1 = expr.rewrite(sin).simplify().trigsimp()

Thank you in advance

Answer

A few points:

1. Don’t mix up numpy and sympy like this unless you know exactly what you are doing. There is no need to use numpy at all here so use e.g. sym.eye(2) and sym.conjugate(val)

2. Don’t use floats unless you have a good reason – use sym.I instead of 1j. Using numpy can potentially introduce floats so don’t do that unless you know what you are doing.

3. Although eigenvals returns a dict in this case you only care about the values of the dict so you can just do list(M.eigenvals()).

4. Although you declare all symbols as real you are using sqrt(u) which is real only if u is positive. Unless you intend for sqrt(u) to be potentially imaginary then u should be declared as positive.

With the above changes your code looks like this:

import sympy as sym
from sympy import *
t = sym.Symbol('t', real=True)
r = sym.Symbol('r', real=True)
u = sym.Symbol('u', positive=True)
d = sym.Symbol('d', real=True)
a = sym.Symbol('a', real=True)

I = sym.eye(2)

T = sym.Matrix ([[t,0],[0,t]])
R = sym.Matrix ([[r,0],[0,r]])
F = sym.Matrix ([[sym.sqrt(u)*sym.exp(-sym.I*d),0],[0,sym.sqrt(u)*sym.exp(-sym.I*d)]])

sym.init_printing()
M = (I - F @ R @ F @ R)
M_inv = M.inv()

Mf = T @ M_inv @ F @ T

Mf = sym.simplify(Mf)

[eig] = Mf.eigenvals()

expr = eig * sym.conjugate(eig)

Now the final expression expr is this:

In [8]: expr
Out[8]: 
               4                
              t ⋅u              
────────────────────────────────
⎛ 2      -2⋅ⅈ⋅d⎞ ⎛ 2      2⋅ⅈ⋅d⎞
⎝r ⋅u - ℯ      ⎠⋅⎝r ⋅u - ℯ     ⎠

In [9]: expr.expand()
Out[9]: 
                  4                   
                 t ⋅u                 
──────────────────────────────────────
 4  2    2    2⋅ⅈ⋅d    2    -2⋅ⅈ⋅d    
r ⋅u  - r ⋅u⋅ℯ      - r ⋅u⋅ℯ       + 1

In [10]: expr.expand().rewrite(sin)
Out[10]: 
                                   4                                    
                                  t ⋅u                                  
────────────────────────────────────────────────────────────────────────
 4  2    2                               2                              
r ⋅u  - r ⋅u⋅(-ⅈ⋅sin(2⋅d) + cos(2⋅d)) - r ⋅u⋅(ⅈ⋅sin(2⋅d) + cos(2⋅d)) + 1

In [11]: expr.expand().rewrite(sin).expand()
Out[11]: 
             4             
            t ⋅u           
───────────────────────────
 4  2      2               
r ⋅u  - 2⋅r ⋅u⋅cos(2⋅d) + 1

I presume that the final result there is the one that you wanted.