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

Question

我有以下代码，我需要在其中求解表达式以找到根。需要求解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])

程序出错，我无法更正。错误是：TypeError: cannot create mpf from 0.005Iomega

编辑 1 - 我尝试根据 cmets 实施不同的方法。第一种方法是使用 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'

如何更正程序？请向我建议可以产生正确结果的方法。

Answer 1

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 放在一个整数周围的。在 Python 中除整数时，它会失去准确性，因为它使用浮点数。将其转换为 SymPy 对象可保持所有准确性。

既然我们知道我们有一个整数上的ImageSet，我们可以开始列出一些解决方案：

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)

凭借一些经验，您可以将其自动化，以便整个程序仅返回 1 个解决方案，无论 solveset 返回的输出类型如何。

方法 2 - SciPy

​​>

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

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)

这给了

    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)

看起来这是 SymPy 找到的解决方案之一。所以我们必须做正确的事。请注意，有许多初始值不会收敛，例如，sol = root(f, [1, 1])。