在 Python 中重写 Mathematica 输出

Question

我想将 Mathematica 中计算的分段函数输出转换为 Python。从 this 页面获取 Mathematica->Python 转换和 this 页面编写分段函数的灵感，我有

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()在内部或者我们可以直接使用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版。