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コードでもメモリエラーを発生させることなく実行できます）、32ビットまで管理しm=7ましn=7た

現在のタイミングを以下に要約できます (m=3 および n=3 で測定)。そこから、数値積分がボトルネックであることがわかります。

試行関数の構築 = 0%
微分方程式の評価 = 2%
k1 のラム化 = 22%
k1 の統合 = 74%
k2 のラム化と統合 = 2%
固有値の抽出 = 0%

関連する質問: lambdify について

4

2 に答える 2

5

計算の別の段階で数値評価に切り替えることで、ラムダ化時間を回避できると思います。

つまり、X が 5x5 の行列 (内部に微分演算がありますが、それは問題ではありません) であり、M が大きい 5xM であるという意味でk1、あなたの計算は対角線のように見えます。したがって。k2k = g^T X ggk[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 = Matrix( [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 = 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 )

テストはしませんでしたが (過去の午前)、アイデアはわかります。

ここで、すべてを遅くする巨大なシンボリック マトリックスを使用する代わりに、試行関数のインデックス用の 2 つのマトリックス インデックスと自由なパラメーターi1,j1を使用i2,j2し、それらの役割を果たし、最後に整数をそれらに代入する必要があります。

ラム化する行列は 5x5 のみであり、すべてのループの外側で 1 回だけラム化する必要があるため、ラム化と単純化のオーバーヘッドはなくなります。さらに、この問題は、大きな m、n の場合でも簡単にメモリに収まります。

統合はそれほど高速ではありませんが、式が非常に小さいため、たとえばFortran でダンプしたり、他のスマートなことを簡単に行うことができます。

于 2013-05-02T23:46:54.587 に答える