python - 函数，SymPy和SciPy矩阵上的数值积分

从我的SymPy输出中，我得到如下所示的矩阵，必须将其集成为2D格式。目前，我正在按元素进行操作，如下所示。此方法有效，但对于我的实际情况(其中sympy.mpmath.quad及其功能要大得多(请参见下面的编辑))，它变得太慢(对于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]

编辑:
实际情况已在this link中提供。只需解压缩并运行shadmehri_2012.py(此示例的作者来自Shadmehri et al. 2012)。
我已经开始奖励50人，只要他能做到以下几点:

使它比建议的问题

更快
即使在代码中包含多个术语m=15和n=15的情况下，

仍能在不给内存错误的情况下运行)，我在32位

中最多管理了m=7和n=7
当前时序可以总结如下(以m = 3和n = 3测量)。从那里可以看出，数值积分是瓶颈。

建立试用功能= 0％
评估微分方程= 2％
lambdifying k1 = 22％
积分k1 = 74％
Lambdify和积分k2 = 2％
提取特征值= 0％

相关问题:about lambdify

最佳答案

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

也就是说，从k1和k2都是k = g^T X g的形式来看，您的计算似乎是对角线的，其中X是一些5x5矩阵(内部有差分运算，但这没关系)，而g是5xM，M大。因此k[i,j] = g.T[i,:] * X * g[:,j]。

所以你可以更换

对于xrange(1，n + 1)中的j:
对于我在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 =矩阵([g1，g2，g3，g4，g5])

和

i1 =符号('i1')
j1 =符号('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 =符号('i2')
j2 =符号('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])

和

tmp =运算值expr(B * g)
k1 = r * tmp.transpose()* F * tmp
k2 = r * g.transpose()* evaluateExpr(Bc * g)
k2 =运算值expr(k2)

经过

tmp_right = EvaluationExpr(B * g_right)
tmp_left = EvaluationExpr(B * g_left)
k1 = r * tmp_left.transpose()* F * tmp_right
k2 = r * g_left.transpose()* evaluateExpr(Bc * g_right)
k2 =运算值expr(k2)

没有测试(上午)，但是您知道了。

现在，不再有一个庞大的符号矩阵，它使所有事情变慢，而是有两个用于试验函数索引的矩阵索引，以及自由参数i1,j1和i2,j2发挥作用，最后应将整数替换为它们。

由于要进行lambdify的矩阵仅为5x5，并且只需要在所有循环外部进行一次lambdify，就消除了lambdification和简化的开销。此外，即使对于较大的m，n，该问题也很容易放入内存中。

集成速度不是很快，但是由于表达式很小，因此您可以轻松地例如dump them in Fortran或做其他聪明的事。

关于python - 函数，SymPy和SciPy矩阵上的数值积分，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/16295140/