python - 使用 mpmath 和 sympy 模块时因指数函数而出错
问题描述
我有以下代码,我需要在其中求解表达式以找到根。需要求解 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.005 I omega
编辑 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
请注意我在除法之前如何放置一个整数。在 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])
.
推荐阅读
- clips - 两个不同定义模板的联合规则
- javascript - 从 api 链接获取 json 的值
- mysql - MYSQL 查询有待优化
- javascript - 如果一个函数在另一个函数内部执行,它会立即成为回调函数吗?
- ssh - 连接失败。用尽了可用的身份验证方法
- asp.net - 检测到 Microsoft.AspNetCore.Http.Abstractions 的版本冲突
- wordpress - htaccess 规则以避免重复以 index.php 结尾
- php - 如何在 c++ 进程中通过 exec 运行 php 脚本并从该脚本中获取 std :: out?
- tfs - 如何还原所有更改并转到源代码的先前版本?获取此特定版本与回滚?
- sql - SQL 时间戳查询:检查一个时间戳是否早于另一个