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

Question

我一直在为一个数学项目探索 scipy 和 scipy 的核心包。

我需要对一些方程式进行微积分运算...所以为了学习 scipy，我决定测试一个简单的方程式 (PDF of a normal random variable)。我需要在微积分运算期间保持常量不变……而不是为其赋值。

我能够使用 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？

Answer 1

此简化代码可能是您实现目标的一步:

它需要一个函数，这里是一个二次函数以尽可能简单，并计算它的封闭形式积分(使用 sympy)

然后它使用 scipy 对 0 到 3 之间的整数值进行数值计算。

并且闭式积分使用与 sympy 相同的极限。

它有点像原型(prototype)，但似乎可以满足您的要求。

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