7

从我的 SymPy 输出中，我得到了如下所示的矩阵，我必须将其集成为 2D。目前我正在按元素进行操作，如下所示。这种方法有效，但对于我的实际情况（其中和它的功能要大得多（见下面的编辑）来说，它变得太慢了（对于sympy.mpmath.quad和）：scipy.integrate.dblquadA

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 et al. 2012）。我已经为可以执行以下操作的人开始了 50 的赏金：

  • 使其比提出的问题合理地快
  • 即使有许多术语m=15和n=15代码也能在不出现内存错误的情况下运行），我管理到m=732n=7位

当前时间可以总结如下（用 m=3 和 n=3 测量）。由此可见，数值积分是瓶颈。

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

相关问题：关于lambdify

4

2 回答 2

5

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

也就是说，从某种意义上说，您的计算似乎是对角线的，k1并且k2都是k = g^T X gX 是某个 5x5 矩阵（内部带有微分运算，但这无关紧要）的形式，并且g是 5xM，M 大。因此k[i,j] = g.T[i,:] * X * g[:,j]。

所以你可以更换

对于 xrange(1,n+1) 中的 j：
    对于 xrange(1,m+1) 中的 i：
        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 = 矩阵（[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 = 矩阵（[g1，g2，g3，g4，g5]）

和

tmp = evaluateExpr(B*g)
k1 = r*tmp.transpose() * F * tmp
k2 = r*g.transpose()*evaluateExpr(Bc*g)
k2 = 评估表达式（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 = 评估表达式（k2）

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

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

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

集成并没有更快，但由于表达式非常小，您可以轻松地将它们转储到 Fortran 中或执行其他智能操作。

于 2013-05-02T23:46:54.587 回答
1

quadpy（我的一个项目）进行矢量化数值积分。这个

from numpy import sin, cos, pi
import quadpy


def f(X):
    x, t = X
    return [
        [(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)]
        ]


x1 = 30
x2 = 75
t1 = 0
t2 = 2*pi

sol = quadpy.quadrilateral.integrate(
        f,
        [[x1, t1], [x2, t1], [x2, t2], [x1, t2]],
        quadpy.quadrilateral.Product(quadpy.line_segment.GaussLegendre(5))
        )

print(sol)

给

[[ 1456.3701526 ]
 [ 2620.60490653]
 [ 5034.5831071 ]]

时间：

%timeit quadpy.quadrilateral.integrate(f, [[x1, t1], [x2, t1], [x2, t2], [x1, t2]], q)
1000 loops, best of 3: 219 µs per loop

这会在您的可下载示例中显着加快速度：

import numpy
array2mat = [{'ImmutableMatrix': numpy.array}, 'numpy']
k1_lambda = lambdify( (x,t), k1, modules=array2mat)
print 'Finished lambdifying k1:', time.clock()
import quadpy
sol = quadpy.quadrilateral.integrate(
    lambda X: k1_lambda(X[0], X[1]),
    [[x1, t1], [x2, t1], [x2, t2], [x1, t2]],
    quadpy.quadrilateral.Product(quadpy.line_segment.GaussLegendre(5))
    )

输出：

Start: 0.040001
Finished trial functions: 0.379929
Finished evaluating differential equations: 2.669536
Finished lambdifying k1: 29.961808
Finished integrating k1: 30.106988
Finished lambdifying and integrating k2: 34.229007
Finished calculating eigenvalues and eigenvectors: 34.229924

请注意，quadpy不进行自适应正交，因此请明智地选择您的方案。

于 2017-03-22T11:05:48.777 回答