Ошибка из-за экспоненциальной функции при использовании модулей 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: невозможно создать mpf из 0,005Iomega

Изменить 1. Я попытался реализовать другой подход на основе комментариев. Первый подход заключался в использовании модуля sympy.solveset. Второй подход заключался в использовании fsolve из scipy.optimise. Оба не дают надлежащего вывода.

Для ясности я копирую соответствующий код для каждого подхода вместе с выводом, который я получаю.

Подход 1 - Симпи


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 Сципи


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'

Как исправить программу? Пожалуйста, предложите мне подход, который даст правильные результаты.


python-3.x python sympy mpmath
person Jan    schedule 04.10.2020    source источник
comment
Если вы используете mpmath.findroot, почему вы используете scipy/numpy в lambdify?   -  person hpaulj    schedule 04.10.2020
comment
Вы говорите о двух разных ошибках, mpf from '0.005/omega' и exponential function. Я получаю только ошибку 'mpc' object has no attribute 'exp'.   -  person hpaulj    schedule 04.10.2020
comment
Просто придерживайтесь либо использования sympy, либо использования mpmath. Хотя sympy использует mpmath, если вы попытаетесь их смешать, вам будет сложно сделать это правильно. SymPy может решить это уравнение с помощью функцииsolve/solveset.   -  person Oscar Benjamin    schedule 04.10.2020
comment
@oscar benjamin Я пытался использовать набор решений из Sympy, но он возвращает пустой набор, если я решаю в реальном домене. Если я решаю в воображаемой области, я получаю эту ошибку. ImageSet(Lambda(_n, -200,0*I*(I*(2*_n*pi + pi) + 0,693147180559945)), целые числа)   -  person Jan    schedule 05.10.2020
comment
Это не правильный ответ? Обратите внимание, что solve вернет список выражений, если вы этого хотите.   -  person Oscar Benjamin    schedule 05.10.2020
comment
Нет, это должно дать какую-то ценность. Можете ли вы предложить какой-либо другой метод?   -  person Jan    schedule 05.10.2020

Ответы (1)


arrow_upward
0
arrow_downward

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)

Имея некоторый опыт, вы можете автоматизировать это так, чтобы вся программа возвращала только одно решение, независимо от типа вывода, возвращаемого 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]).

person Chris du Plessis    schedule 06.10.2020