如何解决使用mpmath和sympy模块时由于指数函数而出错

我有以下代码需要解析表达式以查找根的代码。该表达式需要针对欧米茄进行求解。

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])

程序出现错误，我无法纠正。错误是：TypeError：无法从0.005 I omega

创建mpf

编辑1-我尝试根据评论实施不同的方法。第一种方法是使用sympy.solveset模块。第二种方法是使用scipy.optimise中的fsolve。两者都没有给出正确的输出。

为清楚起见，我将相关代码与获得的输出一起复制到每种方法。

方法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)

方法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'

如何纠正程序？请给我建议可以产生适当结果的方法。

解决方法

SymPy是一个计算机代数系统，可以像人类一样求解方程。 SciPy使用数值优化。如果您需要所有解决方案，建议您使用SymPy。如果您想要一种解决方案，建议您使用SciPy。

方法1-SymPy

SymPy提供的解决方案对于您（作为开发人员）将更具“交互性”。但这几乎总是正确的。

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)

给予

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

这是真正的数学解决方案。

请注意在除以S之前，如何将ImageSet放在整数周围。在Python中对整数进行除法时，由于使用浮点数，因此准确性下降。将其转换为SymPy对象可以保持所有准确性。

由于我们知道整数上有for n in range(-3,3): print(complex(sol.lamda(n)))，因此我们可以开始列出一些解决方案：

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

哪个给

solveset

根据一些经验，您可以将其自动化，以便无论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)返回的输出类型如何，整个程序仅返回1个解决方案。

方法2-SciPy

SciPy提供的解决方案将更加自动化。您将永远找不到完美的答案，并且初始条件的不同选择可能不会一直收敛。

    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)

哪个给

sol = root(f,[1,1])

看起来这是SymPy找到的解决方案之一。因此，我们必须做正确的事。请注意，有许多不收敛的初始值，例如<div class="content">Content</div> <div class="content">Content</div> .content { padding: 5rem 2rem; }。

使用mpmath和sympy模块时由于指数函数而出错