函数矩阵、SymPy 和 SciPy 上的数值积分

从我的 SymPy 输出中，我得到了如下所示的矩阵，我必须将其积分为 2D。目前我正在按元素进行操作，如下所示。此方法有效，但速度太慢（对于sympy.mpmath.quad and scipy.integrate.dblquad）对于我的真实案例（其中A它的功能要大得多（见下面的编辑）：

from sympy import Matrix, sin, cos
import sympy
import scipy
sympy.var( 'x, t' )
A = Matrix([[(sin(2-0.1*x)*sin(t)*x+cos(2-0.1*x)*cos(t)*x)*cos(3-0.1*x)*cos(t)],
            [(cos(2-0.1*x)*sin(t)*x+sin(2-0.1*x)*cos(t)*x)*sin(3-0.1*x)*cos(t)],
            [(cos(2-0.1*x)*sin(t)*x+cos(2-0.1*x)*sin(t)*x)*sin(3-0.1*x)*sin(t)]])

# integration intervals
x1,x2,t1,t2 = (30, 75, 0, 2*scipy.pi)

# element-wise integration
from sympy.utilities import lambdify
from sympy.mpmath import quad
from scipy.integrate import dblquad
A_int1 = scipy.zeros( A.shape, dtype=float )
A_int2 = scipy.zeros( A.shape, dtype=float )
for (i,j), expr in scipy.ndenumerate(A):
    tmp = lambdify( (x,t), expr, 'math' )
    A_int1[i,j] = quad( tmp, (x1, x2), (t1, t2) )
    # or (in scipy)
    A_int2[i,j] = dblquad( tmp, t1, t2, lambda x:x1, lambda x:x2 )[0]

我正在考虑一次性完成它，但我不确定这是否是正确的方法：

A_eval = lambdify( (x,t), A, 'math' )
A_int1 = sympy.quad( A_eval, (x1, x2), (t1, t2)                 
# or (in scipy)
A_int2 = scipy.integrate.dblquad( A_eval, t1, t2, lambda x: x1, lambda x: x2 )[0]

编辑： 真实案例已在这个链接。只需解压并运行shadmehri_2012.py（这个例子的作者来自：沙德梅里等人。 2012年）。 我已经开始悬赏 50 美元奖励那些能够做到以下几点的人：

• 使其比提出的问题更快
  
• 即使有很多术语，也能成功运行而不会出现内存错误m=15 and n=15在代码中），我设法达到m=7 and n=732 位

当前时序可总结如下（以 m=3 和 n=3 测量）。从这里可以看出数值积分是瓶颈。

构建试用功能 = 0%
计算微分方程 = 2%
羔羊化 k1 = 22%
积分 k1 = 74%
羔羊化和积分 k2 = 2%
提取特征值 = 0%

相关问题：关于 lambdify

我认为您可以通过在计算的不同阶段切换到数值计算来避免羔羊化时间。

也就是说，您的计算似乎是对角的，因为k1 and k2都是以下形式k = g^T X g其中 X 是某个 5x5 矩阵（内部有差分运算，但这并不重要），并且g是 5xM，其中 M 大。所以k[i,j] = g.T[i,:] * X * g[:,j].

所以你可以直接替换

for j in xrange(1,n+1):
    for i in xrange(1,m+1):
        g1 += [uu(i,j,x,t),          0,          0,          0,          0]
        g2 += [          0,vv(i,j,x,t),          0,          0,          0]
        g3 += [          0,          0,ww(i,j,x,t),          0,          0]
        g4 += [          0,          0,          0,bx(i,j,x,t),          0]
        g5 += [          0,          0,          0,          0,bt(i,j,x,t)]
g = Matrix( [g1, g2, g3, g4, g5] )
  

with

i1 = Symbol('i1')
j1 = Symbol('j1')
g1 = [uu(i1,j1,x,t),          0,          0,          0,          0]
g2 = [          0,vv(i1,j1,x,t),          0,          0,          0]
g3 = [          0,          0,ww(i1,j1,x,t),          0,          0]
g4 = [          0,          0,          0,bx(i1,j1,x,t),          0]
g5 = [          0,          0,          0,          0,bt(i1,j1,x,t)]
g_right = Matrix( [g1, g2, g3, g4, g5] )

i2 = Symbol('i2')
j2 = Symbol('j2')
g1 = [uu(i2,j2,x,t),          0,          0,          0,          0]
g2 = [          0,vv(i2,j2,x,t),          0,          0,          0]
g3 = [          0,          0,ww(i2,j2,x,t),          0,          0]
g4 = [          0,          0,          0,bx(i2,j2,x,t),          0]
g5 = [          0,          0,          0,          0,bt(i2,j2,x,t)]
g_left = Matrix( [g1, g2, g3, g4, g5] )
  

and

tmp = evaluateExpr( B*g )
k1 = r*tmp.transpose() * F * tmp
k2 = r*g.transpose()*evaluateExpr(Bc*g)
k2 = evaluateExpr( k2 )
  

by

tmp_right = evaluateExpr( B*g_right )
tmp_left = evaluateExpr( B*g_left )
k1 = r*tmp_left.transpose() * F * tmp_right
k2 = r*g_left.transpose()*evaluateExpr(Bc*g_right)
k2 = evaluateExpr( k2 )
  

没有测试（上午过去），但你明白了。

现在，不再有一个巨大的符号矩阵，这会让一切变慢，而是有两个用于试验函数索引的矩阵索引和自由参数i1,j1 and i2,j2它们发挥了作用，最后您应该将整数替换为它们。

由于要进行羔羊化的矩阵只有 5x5，并且只需要在所有循环之外进行一次羔羊化，因此羔羊化和简化的开销就消失了。此外，即使对于较大的 m、n，该问题也很容易装入内存。

集成并没有更快，但由于表达式非常小，您可以轻松地例如将它们转储到 Fortran 中或者做一些其他聪明的事情。