Переписать вывод Mathematica в Python

Question

Я хочу преобразовать кусочно-функциональный вывод вычисления в Mathematica в Python.Черпая вдохновение из этой страницы для преобразования Mathematica-> Python и этой страницы для написания кусочных функций, у меня есть

from numpy import linspace, vectorize, array
from numpy import arctan, log
import matplotlib.pyplot as plt
from sympy.parsing.mathematica import parse

def fun(x,a,b,c):
    # string inside parse('') is my mathematica output
    if a == 0:
        out = parse('a (-I b + x) ArcTan[(b - a)/c]')
    else:
        out = parse('4 c a^2 Log[c^2 + (x - a)^2]')
    return out

a = 0.17
b = 0.44
c = 0.29
x = linspace(0,50,int(1e3))

vfun = vectorize(fun)
y = vfun(x,a,b,c)
yp = 4*c*a**2*log(c**2 + (x - a)**2)

plt.figure()
plt.plot(x,y,'.-')
plt.title('mathematica -> python conversion')

plt.figure()
plt.plot(x,yp,'.-')
plt.title('expected')

plt.show()

Сюжет выглядит так:

тогда как должно быть:

Я сделал что-то не такпри преобразовании Mathematica в Python?Или есть проблема при присвоении числовых значений a, b и c?(Обратите внимание, что это MWE, и код Mathematica, который я хочу преобразовать, намного длиннее и сложнее, чем показано выше.)

Answer 1

Гораздо более простым решением является использование eval.Теперь я знаю, что eval очень опасно, но большую часть времени, когда он требует ввода от пользователя, но здесь ввод уже определен для нас.

Теперь на ответ.
Вы не получаетеожидаемый результат, потому что ваш векторизованный массив не содержит никаких чисел с плавающей запятой, но содержит только выходные данные из синтаксического анализатора mathematica, который возвращает неоцененную строку, поэтому мы можем использовать .evalf() в качестве ответа, данного @David, lambdify(), который использует eval() для внутреннего использования или мы можем напрямуюuse eval().
Вот документация об используемых методах: https://docs.sympy.org/latest/tutorial/basic_operations.html

from numpy import linspace, vectorize, array
from numpy import arctan, log
import matplotlib.pyplot as plt
from sympy.parsing.mathematica import MathematicaParser

def fun(x,a,b,c):
    obj = MathematicaParser()
    if a == 0:
        out = obj.parse('a (-I b + x) ArcTan[(b - a)/c]')
    else:
        out = obj.parse('4 c a^2 Log[c^2 + (x - a)^2]')

    return out

a = 0.17
b = 0.44
c = 0.29
x = linspace(0,50,int(1e3))
yp = 4*c*a**2*log(c**2 + (x - a)**2)
vfun = vectorize(fun)
y = vfun(x,a,b,c)
#Here y is a type <class 'numpy.ndarray'> containing 1000 <class 'numpy.str_'>
#['4*c*a**2*log(c**2+(x-a)**2)' '4*c*a**2*log(c**2+(x-a)**2)'
#'4*c*a**2*log(c**2+(x-a)**2)' '4*c*a**2*log(c**2+(x-a)**2)'
#'4*c*a**2*log(c**2+(x-a)**2)' '4*c*a**2*log(c**2+(x-a)**2)'
#'4*c*a**2*log(c**2+(x-a)**2)' '4*c*a**2*log(c**2+(x-a)**2)'
#....
y = eval(y[0]) #y[0] is '4*c*a**2*log(c**2+(x-a)**2)'
#When we evaluate y[0] we again get y as <class 'numpy.ndarray'> conatining 1000 <class 'numpy.float64'>
#Because x is <class 'numpy.ndarray'> it evaluates the first string over
#all the points in x.
#[-0.07309464 -0.07770262 -0.08110382 -0.0828403  -0.08263539 -0.08052339
#-0.07683235 -0.07203573 -0.06659307 -0.06086366 -0.05509179 -0.04942739
#-0.04395413 -0.03871304 -0.03371924 -0.0289728  -0.02446552 -0.02018495
#.....
plt.figure()
plt.plot(x, y,'.-')
plt.title('mathematica -> python conversion')

plt.figure()
plt.plot(x,yp,'.-')
plt.title('expected')

plt.show()

Выход:

Answer 2

Это грязное решение:

import numpy as np
import matplotlib.pyplot as plt
from sympy.parsing.mathematica import mathematica
from sympy import symbols

def fun(x, a, b, c):
    # string inside parse('') is my mathematica output
    if a == 0:
        out = mathematica('a (-I b + x) ArcTan[(b - a)/c]')
    else:
        out = mathematica('4 c a^2 Log[c^2 + (x - a)^2]')
    return out

aa = 0.17
bb = 0.44
cc = 0.29
II = 1
xx = np.linspace(0, 50, int(1e3))

x, a, b, c, I = symbols('x, a, b, c, I')

fun1 = fun(x, a, b, c)
fun2 = fun1.subs({a: aa, c: cc})
print(fun2.evalf())

y_list = []

for item in xx:
    y_list.append(fun2.subs({x:item}).evalf())

print(y_list[:10])

plt.figure()
plt.plot(xx, y_list,'.-')
plt.title('mathematica -> python conversion')
plt.show()

Я объясню позже, как все это работает.

Обновление

Как вы можете видеть, когда a == 0, функция имеет значение 0, вам нужно проверить это.

Когда вы используете mathematica, тип fun1 является функцией sympy (sympy.core.mul.Mul), поэтому вы должны использовать symbols и fun1.subs() и fun2.evalf(), другими словами, вам нужно знать, кто использует основные функции sympy.

То, как вы оцениваете функцию для ее построения, ну, вы используете список:

y_list = []

for item in xx:
    y_list.append(fun2.subs({x:item}).evalf())

Кстати, я использую sympy версию 1.4.