Error because of exponential function when using mpmath and sympy modules

Question

I have the following code where I need to solve an expression to find the roots. The expression needs to be solved for omega.

import numpy as np
from sympy import Symbol,lambdify
import scipy
from mpmath import findroot, exp
eta = 1.5 
tau = 5 /1000
omega = Symbol("omega")
Tf = exp(1j * omega * tau)
symFun = 1 + Tf * (eta - 1) 
denom = lambdify((omega), symFun, "scipy")
Tf_high = 1j * 2 * np.pi * 1000 * tau
sol = findroot(denom, [0+1j,Tf_high])

The program gives an error and I am not able to correct. The error is : TypeError: cannot create mpf from 0.005Iomega

Edit 1 - I have tried to implement different approach based on comments. First approach was to use the sympy.solveset module. Second approach was to use fsolve from scipy.optimise. Both are not giving proper output.

For clarity, I am copying the relevant code to each approach along with the output I am getting.

Approach 1 - Sympy

import numpy as np
from sympy import Symbol,exp
from sympy.solvers.solveset import solveset,solveset_real,solveset_complex
import matplotlib.pyplot as plt 

def denominator(eta,Tf):
    
    return 1 + Tf * (eta - 1)

if __name__ == "__main__":
    eta = 1.5 
    tau = 5 /1000
    omega = Symbol("omega")
    n = 1 
    Tf = exp(1j * omega * tau)
    denom = 1 + Tf * (eta - 1)
    symFun = denominator(eta,Tf)
    sol = solveset_real(denom,omega)
    sol1 = solveset_complex(denom,omega)
    print('In real domain', sol)
    print('In imaginary domain',sol1)

Output: 
In real domain EmptySet
In imaginary domain ImageSet(Lambda(_n, -200.0*I*(I*(2*_n*pi + pi) + 0.693147180559945)), Integers)

Approach 2 Scipy

import numpy as np
from scipy.optimize import fsolve, root

def denominator(eta,tau,n, omega):
    
    Tf = n * np.exo(1j * omega * tau)
    return 1 + Tf * (eta - 1)

if __name__ == "__main__":
    eta = 1.5 
    tau = 5 /1000
    n = 1 
    func = lambda omega :  1 + (eta - 1) * (n * np.exp( 1j * omega * tau))
    sol = fsolve(func,10)
    print(sol)

Output: 
Cannot cast array data from dtype('complex128') to dtype('float64') according to the rule 'safe'

How do I correct the program? Please suggest me the approach that will give proper results.

Answer 1

SymPy is a computer algebra system and solves the equation like a human would. SciPy uses numeric optimization. If you want ALL the solutions, I suggest going with SymPy. If you want one solution, I suggest going with SciPy.

Approach 1 - SymPy

The solutions SymPy gives will be more "interactive" for you as the developer. But it will be perfectly correct almost all the time.

from sympy import *

eta = S(3)/2
tau = S(5) / 1000
omega = Symbol("omega")
n = 1
Tf = exp(I * omega * tau)
denom = 1 + Tf * (eta - 1)
sol = solveset(denom, omega)
print(sol)

Giving

ImageSet(Lambda(_n, -200*I*(I*(2*_n*pi + pi) + log(2))), Integers)

This is the true mathematical solution.

Notice how I put S around an integer before dividing it. When dividing integers in Python, it loses accuracy because it uses floating point numbers. Converting it to SymPy objects keep all the accuracy.

Since we know we have an ImageSet over integers, we can start listing a few solutions:

for n in range(-3, 3):
    print(complex(sol.lamda(n)))

Which gives

(-3141.5926535897934-138.62943611198907j)
(-1884.9555921538758-138.62943611198907j)
(-628.3185307179587-138.62943611198907j)
(628.3185307179587-138.62943611198907j)
(1884.9555921538758-138.62943611198907j)
(3141.5926535897934-138.62943611198907j)

With some experience, you could automate it so that the whole program only returns 1 solution no matter on the type of output returned by solveset.

Approach 2 - SciPy

The solutions SciPy gives will be more automated. You will never have a perfect answer and different choices of the initial conditions may not converge all the time.

import numpy as np
from scipy.optimize import root

eta = 1.5
tau = 5 / 1000
n = 1
def f(omega: Tuple):
    omega_real, omega_imag = omega
    omega: complex = omega_real + omega_imag*1j
    result: complex = 1 + (eta - 1) * (n * np.exp(1j * omega * tau))
    return result.real, result.imag
sol = root(f, [100, 100])
print(sol)
print(sol.x[0]+sol.x[1]*1j)

Which gives

    fjac: array([[ 0.00932264,  0.99995654],
       [-0.99995654,  0.00932264]])
     fun: array([-2.13074003e-12, -8.86389816e-12])
 message: 'The solution converged.'
    nfev: 30
     qtf: array([ 2.96274855e-09, -6.82780898e-10])
       r: array([-0.00520194, -0.00085702, -0.00479143])
  status: 1
 success: True
       x: array([ 628.31853072, -138.62943611])

(628.3185307197314-138.62943611241522j)

Looks like that's one of the solutions SymPy found. So we must be doing something right. Note that there are many initial values that don't converge, for example, sol = root(f, [1, 1]).