函数矩阵,SymPy和SciPy的数值积分

Question

从我的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]

---

编辑:真实案例已在此链接中提供.只需解压缩并运行shadmehri_2012.py(这个例子的作者来自:Shadmehri等人,2012).我已经为可以执行以下操作的人开始了50美元的奖励:

• 使它比提出的问题合理地快
  
• 设法不给内存错误即使有一些术语的运行m=15,并n=15在代码中),我设法达到m=7和n=732位

当前时序可以总结如下(用m = 3和n = 3测量).从那里可以看出,数值积分是瓶颈.

构建试验函数= 0%
评估微分方程= 2%
lambdifying k1 = 22%
积分k1 = 74%
lambdifying和积分k2 = 2%
提取特征值= 0%

---

相关问题:关于lambdify

Answer 1

我认为你可以通过在计算的不同阶段切换到数值评估来避免lambdification时间.

也就是说,你的计算在某种意义上似乎是对角的,k1并且k2都是k = g^T X gX是一些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] )

同

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] )

和

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

通过

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和i2,j2它们扮演它们的角色,你最后应该将整数替换成它们.

由于lambdify的矩阵只有5x5,并且只需要在所有循环之外进行一次lambdified,因此lambdification和简化开销就消失了.而且,即使对于大的m,n,该问题也容易适合存储器.

集成并不是更快,但由于表达式非常小,您可以轻松地将它们转储到Fortran或者做其他聪明的事情.