如何在scipy中创建数学表达式？

Question

我一直在探索Scipy和其核心包，用于一个数学项目。

我需要对一些方程进行微积分运算……因此为了学习Scipy，我决定测试一个简单的方程（正态随机变量的概率密度函数）。在微积分运算中，我需要保持常数不变，而不是将其赋值。

我使用Sympy成功地创建了它。以下是代码：

from sympy import *

x = Symbol('x')
mu = Symbol('mu')
sigma = Symbol('sigma')

def normpdfeqn():
    y = exp(-(((x-mu)**2)) / (2*(sigma**2))) / (sqrt(2*pi*(sigma**2)))
    return y

print(integrate(normpdfeqn(), (x)))

然后我得到了适当的输出：

sigma*erf(sqrt(2)*(-mu + x)/(2*sigma))/(2*sqrt(sigma**2))

接着我尝试使用scipy做同样的事情。

我一直在阅读http://docs.scipy.org/doc/scipy/reference/tutorial/integrate.html，但是无法弄清如何为其创建方程。以下是我到目前为止尝试过的内容（几乎与上面的代码相同）：

from sympy import exp
from sympy import sqrt
from sympy import pi
from scipy.integrate import quad
from sympy import Symbol
x = Symbol('x')
mu = Symbol('mu')
sigma = Symbol('sigma')

def integrand():
    y = exp(-((x-mu)**2) / (2*(sigma**2))) / (sqrt(2*pi*(sigma**2)))
    return y


I = quad(integrand(), 0, 1,)
print(I)

这段代码可能远未完成，我不知道如何让它工作。

如果我总是使用上述展示的方程式，是否需要学习scipy integrate？还是继续使用sympy和numpy？

- Jeet Parekh

---

如果你想进行“符号”操作，你应该使用sympy。Scipy用于进行“数值”操作。 - BrenBarn

@BrenBarn 所以，我使用sympy创建方程式，然后将它们lambdify（或者是其他一些函数，我记不清了），以便可以在numpy或scipy中使用该方程式...这个工作流程可以吗？ - Jeet Parekh

非常有趣。大多数人将这两个努力分开，因为他们的方程很少改变。但是你打算进行分析性输入，然后让代码首先进行分析性分析（无监督），然后计算数值结果？这可能是雄心勃勃的。我喜欢它。 - roadrunner66

@roadrunner66 我正在和我的数学教授一起研究一个项目，该项目通过改变一些因素来研究信号的变化。这就是为什么我的方程式会经常变化的原因。 - Jeet Parekh

@vicky96：我认为你可以做到，没错。 - BrenBarn

@BrenBarn 好的，非常感谢。 - Jeet Parekh

Answer 1

这段简化的代码可能是实现你目标的一步：

它使用一个函数，这里是二次方程以尽可能简单的方式，并使用sympy计算其闭合形式积分。

然后使用scipy在0和3之间数值地评估积分值。

再使用相同的限制条件使用sympy计算闭合形式积分。

这是一个小型的原型，但似乎可以满足你的要求。

from sympy import *
from scipy.integrate import quad

x = Symbol('x')

def func():
    y = x**2 + 2*x +3
    return y

def integrand(expr):
    y = lambda x: eval(expr)
    return y 


y = repr(func())             # a str representation of the original function

res = integrate(func(), (x)) # a sympy representation of the closed form integral 
                             # of the original function

I = quad(integrand(y), 0, 3) # a scipy numerical integration of the original function

print(y)        # sympy repr of original function
print(I)        # scipy numerical integration of the original function
print()
print(res)      # sympy repr of closed form integral of original function
print(integrand(repr(res))(3) - integrand(repr(res))(0)) # closed form integer evaluation

结果输出：

x**2 + 2*x + 3
(27.0, 2.9976021664879227e-13)

x**3/3 + x**2 + 3*x
27.0

- Reblochon Masque