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

本文介绍了函数，SymPy和SciPy矩阵上的数值积分的处理方法，对大家解决问题具有一定的参考价值，需要的朋友们下面随着小编来一起学习吧！

问题描述

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

From my SymPy output I have the matrix shown below, which I must integrate in 2D. Currently I am doing it element-wise as shown below. This method works but it gets too slow (for both sympy.mpmath.quad and scipy.integrate.dblquad) for my real case (in which A and its functions are much bigger (see edit below):

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]

我当时正考虑一次拍摄，但是我不确定这是否可行:

I was considering doing it in one shot like, but I'm not sure if this is the way to go:

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人，只要他能做到以下几点:

The real case has been made available in this link. Just unzip and run shadmehri_2012.py (is the author from were this example was taken from: Shadmehri et al. 2012).I've started a bounty of 50 for the one who can do the following:

使其比提出的问题合理地
即使在代码中包含多个术语m=15和n=15的情况下，也可以在不给出内存错误的情况下运行)，我最多可以在32位中管理m=7和n=7

make it reasonably faster than the proposed question
manage to run without giving memory error even with a number of terms m=15 and n=15 in the code), I managed up to m=7 and n=7 in 32-bit

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

The current timing can be summarized below(measured with m=3 and n=3). From there it can be seen that the numerical integration is the bottleneck.

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

build trial functions = 0%
evaluating differential equations = 2%
lambdifying k1 = 22%
integrating k1 = 74%
lambdifying and integrating k2 = 2%
extracting eigenvalues = 0%

推荐答案

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

I think you can avoid the lambdification time by switching to numerical evaluation at a different stage of the calculation.

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

Namely, your calculation seems to be diagonal in the sense that k1 and k2 are both of the form k = g^T X g where X is some 5x5 matrix (with differential ops inside, but that doesn't matter), and g is 5xM, with M large. Therefore 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 )

(虽然是上午)没有测试，但是您知道了.

Didn't test (past am), but you get the idea.

现在，您有了两个用于试验函数索引的矩阵索引，以及自由参数i1,j1和i2,j2，它们发挥了作用，并且您应该在其中替换整数，而不是使用一个庞大的符号矩阵来使所有事情变慢.结束.

Now, instead of having a huge symbolic matrix which makes everything slow, you have two matrix indices for the trial function indices, and free parameters i1,j1 and i2,j2 which play their role and you should substitute integers into them in the end.

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

Since the matrix to lambdify is only 5x5, and needs to be lambdified only once outside all loops, the lambdification and simplification overhead is gone. Moreover, the problem fits easily into memory even for large m, n.

集成速度不是很快，但是由于表达式很小，因此您可以轻松地例如将其转储到Fortran中或执行其他一些巧妙的操作.

The integration is not any faster, but since the expressions are very small, you can easily e.g. dump them in Fortran or do something else smart.

这篇关于函数，SymPy和SciPy矩阵上的数值积分的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持！