python - 使用 SciPy 集成返回矩阵或数组的函数

Question

我有一个符号数组，可以表示为:

from sympy import lambdify, Matrix

g_sympy = Matrix([[   x,  2*x,  3*x,  4*x,  5*x,  6*x,  7*x,  8*x,   9*x,  10*x],
                  [x**2, x**3, x**4, x**5, x**6, x**7, x**8, x**9, x**10, x**11]])

g = lambdify( (x), g_sympy )

因此对于每个 x 我得到一个不同的矩阵:

g(1.) # matrix([[  1.,   2.,   3.,   4.,   5.,   6.,   7.,   8.,   9.,  10.],
      #         [  1.,   1.,   1.,   1.,   1.,   1.,   1.,   1.,   1.,   1.]])
g(2.) # matrix([[  2.00e+00,   4.00e+00,   6.00e+00,   8.00e+00,   1.00e+01, 1.20e+01,   1.40e+01,   1.60e+01,   1.80e+01,   2.00e+01],
      #         [  4.00e+00,   8.00e+00,   1.60e+01,   3.20e+01,   6.40e+01, 1.28e+02,   2.56e+02,   5.12e+02,   1.02e+03,   2.05e+03]])

等等……


我需要在 x 上对 g 进行数值积分，比如说 从 0. 到 100.(在实际情况下积分没有精确解) 并且在我目前的方法中，我必须 lambdify g 中的每个元素并将其单独集成。我正在使用 quad 进行元素集成，例如:

ans = np.zeros( g_sympy.shape )
for (i,j), func_sympy in ndenumerate(g_sympy):
    func = lambdify( (x), func_sympy)
    ans[i,j] = quad( func, 0., 100. )

这里有两个问题:1)多次使用lambdify和2)for循环；我相信第一个是瓶颈，因为 g_sympy 矩阵最多有 10000 个项(这对 for 循环来说不是什么大问题)。

如上所示，lambdify 允许对整个矩阵求值，所以我想:“有没有办法整合整个矩阵？”

scipy.integrate.quadrature 有一个参数 vec_func 这给了我希望。我期待的是这样的:

g_int = quadrature( g, x1, x2 )

得到完全集成的矩阵，但它给出了ValueError: matrix must be 2-dimensional

---

编辑:我正在尝试做什么 can apparently be done in Matlab使用 quadv 和 has already been discussed for SciPy

真实案例has been made available here .

要运行它，您需要:

• NumPy
• 科学
• matplotlib
• 同情

只需运行:“python curved_beam_mrs.py”。

您会看到该过程已经很慢，主要是因为集成，由文件 curved_beam.py 中的 TODO 指示。

如果您删除文件 curved_beam_mrs.py 中 TODO 之后指示的注释，速度会慢很多。

print.txt 文件中显示了集成的函数矩阵。

谢谢!

Answer 1

quad 或 quadrature 的第一个参数必须是可调用的。 quadrature 的 vec_func 参数是指此可调用对象的参数 是否为(可能是多维的)向量。从技术上讲，您可以向量化 quad 本身:

>>> from math import sin, cos, pi
>>> from scipy.integrate import quad
>>> from numpy import vectorize
>>> a = [sin, cos]
>>> vectorize(quad)(a, 0, pi)
(array([  2.00000000e+00,   4.92255263e-17]), array([  2.22044605e-14,   2.21022394e-14]))

但这只是等同于对 a 的元素进行显式循环。具体来说，如果您想要的话，它不会给您带来任何性能提升。所以，总而言之，问题是您为什么要在这里实现以及您究竟想要实现什么。