使用 SciPy 集成返回矩阵或数组的函数答案

【问题标题】：using SciPy to integrate a function that returns a matrix or array使用 SciPy 集成返回矩阵或数组的函数
【发布时间】：2013-06-08 20:25:06
【问题描述】：

我有一个符号数组，可以表示为：

from sympy import lambdify, Matrix

g_sympy = Matrix([[   x,  2*x,  3*x,  4*x,  5*x,  6*x,  7*x,  8*x,   9*x,  10*x],
                  [x**2, x**3, x**4, x**5, x**6, x**7, x**8, x**9, x**10, x**11]])

g = lambdify( (x), g_sympy )

所以对于每个x，我都会得到一个不同的矩阵：

g(1.) # matrix([[  1.,   2.,   3.,   4.,   5.,   6.,   7.,   8.,   9.,  10.],
      #         [  1.,   1.,   1.,   1.,   1.,   1.,   1.,   1.,   1.,   1.]])
g(2.) # matrix([[  2.00e+00,   4.00e+00,   6.00e+00,   8.00e+00,   1.00e+01, 1.20e+01,   1.40e+01,   1.60e+01,   1.80e+01,   2.00e+01],
      #         [  4.00e+00,   8.00e+00,   1.60e+01,   3.20e+01,   6.40e+01, 1.28e+02,   2.56e+02,   5.12e+02,   1.02e+03,   2.05e+03]])

等等...

我需要对g 与x 进行数值积分，比如from 0. to 100.（在实际情况下，积分没有精确解），在我目前的方法中，我必须lambdifyg 中的每个元素并积分它单独。我正在使用quad 进行元素集成，例如：

ans = np.zeros( g_sympy.shape )
for (i,j), func_sympy in ndenumerate(g_sympy):
    func = lambdify( (x), func_sympy)
    ans[i,j] = quad( func, 0., 100. )

这里有两个问题：1）lambdify用了很多次和2）for循环；我相信第一个是瓶颈，因为g_sympy 矩阵最多有 10000 个术语（这对 for 循环来说没什么大不了的）。

如上图lambdify允许对整个矩阵求值，于是我想：“有没有办法对整个矩阵求积分？”

scipy.integrate.quadrature 有一个参数vec_func 给了我希望。我期待的是这样的：

g_int = quadrature( g, x1, x2 )

要获得完全集成的矩阵，但它给出的ValueError: 矩阵必须是二维的

编辑：我正在尝试使用can apparently be done in Matlab 和quadv 和has already been discussed for SciPy 做的事情

真实案例has been made available here。

要运行它，您需要：

numpy
scipy
matplotlib
同情

只需运行："python curved_beam_mrs.py"。

你会看到这个过程已经很慢了，主要是因为集成，文件curved_beam.py中的TODO表示。

如果您删除文件curved_beam_mrs.py 中TODO 之后指示的注释，它会慢得多。

集成的函数矩阵显示在print.txt文件中。

谢谢！

【问题讨论】：

你能用数学术语来说明积分矩阵是什么意思吗？结果应该是什么，是矩阵还是标量或其他什么？
我将 x 从 0. 更改为 100.，并且我想整合矩阵中的每个元素，就像做 quad( g[i,j], 0., 100. ) for (i,j),v in ndenumerate(g) 一样，但这是我试图避免的元素方式。 ..
代码在 cmets 中看起来很尴尬，所以添加了一个答案 --- 实际上更多的是一个非答案:-)。
您在这里寻找什么样的准确度？
@flebool 问题已更新！

标签： python matrix numpy scipy numerical-integration

【解决方案1】：

quad 或 quadrature 的第一个参数必须是可调用的。 quadrature 的 vec_func 参数是指此可调用对象的参数是否是（可能是多维的）向量。从技术上讲，您可以 vectorize quad 本身：

>>> from math import sin, cos, pi
>>> from scipy.integrate import quad
>>> from numpy import vectorize
>>> a = [sin, cos]
>>> vectorize(quad)(a, 0, pi)
(array([  2.00000000e+00,   4.92255263e-17]), array([  2.22044605e-14,   2.21022394e-14]))

但这仅相当于显式循环a 的元素。具体来说，如果这就是您所追求的，它不会给您带来任何性能提升。所以，总而言之，问题是你为什么以及你想要在这里实现什么。

【讨论】：

感谢您的回答。我正在寻找性能提升。但是您的解决方案比我现在使用的 for 循环更优雅...
如果我将a = [sin, cos] 替换为a = lambda t: [sin(t), cos(t)]，则代码不再起作用。我应该如何调整代码以使其适用于a = lambda t: [sin(t), cos(t)]？

【解决方案2】：

向量化梯形和辛普森积分规则。 Trapezoidal 只是从另一个使用 logspace 而不是 linspace 的项目复制和粘贴，以便它可以利用非均匀网格。

def trap(func,a,b,num):
    xlinear=np.linspace(a,b,num)
    slicevol=np.diff(xlinear)
    output=integrand(xlinear)
    output=output[:,:-1]+output[:,1:]
    return np.dot(output,slicevol)/2

def simpson(func,a,b,num):
    a=float(a)
    b=float(b)
    h=(b-a)/num

    output=4*np.sum(integrand(a+h*np.arange(1,num,2)),axis=1)
    output+=2*np.sum(integrand(a+h*np.arange(2,num-1,2)),axis=1)
    output+=np.sum(integrand(b),axis=1)
    output+=np.sum(integrand(a),axis=1)
    return output*h/3

def integrand(rlin):
    first=np.arange(1,11)[:,None]
    second=np.arange(2,12)[:,None]
    return np.vstack((rlin*first,np.power(rlin,second)))

检查梯形和辛普森规则累积相对误差：

b=float(100)
first=np.arange(1,11)*(b**2)/2
second=np.power(b,np.arange(3,13))/np.arange(3,13)
exact=np.vstack((first,second))

for x in range(3):
    num=x*100+100
    tr=trap(integrand,0,b,num).reshape(2,-1)
    si=simpson(integrand,0,b,num).reshape(2,-1)
    rel_trap=np.sum(abs((tr-exact)/exact))*100
    rel_simp=np.sum(abs((si-exact)/exact))*100
    print 'Number of points:',num,'Trap Rel',round(rel_trap,6),'Simp Rel',round(rel_simp,6)

Number of points: 100 Trap Rel 0.4846   Simp Rel 0.000171
Number of points: 200 Trap Rel 0.119944 Simp Rel 1.1e-05
Number of points: 300 Trap Rel 0.053131 Simp Rel 2e-06

时间。请注意，两个梯形规则都使用 200 点，而基于上述收敛性，辛普森一家的时间仅为 100。对不起，我没有同情：

s="""
import numpy as np
from scipy.integrate import trapz

def integrand(rlin):
    first=np.arange(1,11)[:,None]
    second=np.arange(2,12)[:,None]
    return np.vstack((rlin*first,np.power(rlin,second)))

def trap(func,a,b,num):
    xlinear=np.linspace(a,b,num)
    slicevol=np.diff(xlinear)
    output=integrand(xlinear)
    output=output[:,:-1]+output[:,1:]
    return np.dot(output,slicevol)/2

def simpson(func,a,b,num):
    a=float(a)
    b=float(b)
    h=(b-a)/num

    output=4*np.sum(integrand(a+h*np.arange(1,num,2)),axis=1)
    output+=2*np.sum(integrand(a+h*np.arange(2,num-1,2)),axis=1)
    output+=np.sum(integrand(b),axis=1)
    output+=np.sum(integrand(a),axis=1)
    return output*h/3

def simpson2(func,a,b,num):
    a=float(a)
    b=float(b)
    h=(b-a)/num
    p1=a+h*np.arange(1,num,2)
    p2=a+h*np.arange(2,num-1,2)
    points=np.hstack((p1,p2,a,b))
    mult=np.hstack((np.repeat(4,p1.shape[0]),np.repeat(2,p2.shape[0]),1,1))
    return np.dot(integrand(points),mult)*h/3

def x2(x):
    return x**2
def x3(x):
    return x**3
def x4(x):
    return x**4
def x5(x):
    return x**5
def x5(x):
    return x**5
def x6(x):
    return x**6
def x7(x):
    return x**7
def x8(x):
    return x**8
def x9(x):
    return x**9
def x10(x):
    return x**10
def x11(x):
    return x**11
def xt5(x):
    return 5*x
"""

zhenya="""
a=[xt5,xt5,xt5,xt5,xt5,xt5,xt5,xt5,xt5,xt5,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11]
vectorize(quad)(a, 0, 100)
"""

usethedeathstar="""
g=lambda x: np.array([[x,2*x,3*x,4*x,5*x,6*x,7*x,8*x,9*x,10*x],[x**2,x**3,x**4,x**5,x**6,x**7,x**8,x**9,x**10,x**11]])
xv=np.linspace(0,100,200)
trapz(g(xv))
"""

vectrap="""
trap(integrand,0,100,200)
"""

vecsimp="""
simpson(integrand,0,100,100)
"""

vecsimp2="""
simpson2(integrand,0,100,100)
"""

print 'zhenya took',timeit.timeit(zhenya,setup=s,number=100),'seconds.'
print 'usethedeathstar took',timeit.timeit(usethedeathstar,setup=s,number=100),'seconds.'
print 'vectrap took',timeit.timeit(vectrap,setup=s,number=100),'seconds.'
print 'vecsimp took',timeit.timeit(vecsimp,setup=s,number=100),'seconds.'
print 'vecsimp2 took',timeit.timeit(vecsimp2,setup=s,number=100),'seconds.'

结果：

zhenya took 0.0500509738922 seconds.
usethedeathstar took 0.109386920929 seconds.
vectrap took 0.041011095047 seconds.
vecsimp took 0.0376999378204 seconds.
vecsimp2 took 0.0311458110809 seconds.

在时间上需要指出的是，真雅的答案应该更准确。我相信一切都是正确的，如果需要更改，请告诉我。

如果您提供您将使用的功能和范围，我可能会为您的系统提供更好的东西。您是否有兴趣使用额外的核心/节点？

【讨论】：

非常感谢您的回答。我已经用我正在整合的真实案例更新了这个问题......

【解决方案3】：

在实际情况下，积分没有精确解，您是指奇点吗？您能否更精确地了解它以及您希望整合的矩阵的大小。我不得不承认 sympy 在某些事情上非常慢（不确定集成是否是其中的一部分，但我更喜欢远离 sympy 并坚持使用 numpy 解决方案）。您想通过矩阵还是更快的矩阵来获得更优雅的解决方案？

-注意：显然我不能在你的帖子中添加评论来询问这些东西，所以我不得不发布这个作为答案，也许这是因为我没有足够的声誉左右？-

编辑：像这样？

    import numpy
    from scipy.integrate import trapz
    g=lambda x: numpy.array([[x,2*x,3*x],[x**2,x**3,x**4]])
    xv=numpy.linspace(0,100,200)
    print trapz(g(xv))

已经看到你想要整合 sum(a*sin(bx+c)^n*cos(dx+e)^m) 之类的东西，用于 a、b、c、d、e 的不同系数， m,n，我建议分析性地做所有这些。（应该有一些公式，因为您可以将 sin 重写为复指数

我在检查这些函数时注意到的另一件事是，sin(a*x+pi/2) 和 sin(a*x+pi) 以及类似的东西可以以某种方式重写为 cos 或 sin删除 pi/2 或 pi。另外我看到的是，只需查看函数矩阵中的第一个元素：

a*sin(bx+c)^2+d*cos(bx+c)^2 = a*(sin^2+cos^2)+(d-a)*cos(bx+c)^2 = a+(d-a)*cos(bx+c)^2

这也简化了计算。如果您以不涉及大量 txtfile 左右的方式获得公式，请检查您需要集成的最通用公式是什么，但我猜它类似于 a*sin^n(bx+c)*cos^ m(dx+e)，其中 m 和 n 为 0 1 或 2，这些东西可以简化为可以解析积分的东西。所以如果你找到你得到的最通用的分析函数，你可以很容易地做出类似的东西

f=lambda x: [[s1(x),s2(x)],[s3(x),s4(x)]]
res=f(x2)-f(x1)

其中 s1(x) 等只是您的函数的分析集成版本？

（不是真的打算检查整个代码以查看其余代码的作用，但它只是将这些函数集成到 txt 文件中从 a 到 b 或类似的东西？或者有类似的东西你拿每个函数的平方或任何可能破坏分析的可能性的东西？）

我猜这应该可以简化你的积分？

first integral 和：第二个

嗯，第二个链接不起作用，但我猜你从第一个链接中得到了想法

编辑，因为您不想要分析解决方案：改进仍然在于摆脱 sympy：

from sympy import sin as SIN
from numpy import sin as SIN2
from scipy.integrate import trapz
import time
import numpy as np

def integrand(rlin):
    first=np.arange(1,11)[:,None]
    second=np.arange(2,12)[:,None]
    return np.vstack((rlin*first,np.power(rlin,second)))

def simpson2(func,a,b,num):
    a=float(a)
    b=float(b)
    h=(b-a)/num
    p1=a+h*np.arange(1,num,2)
    p2=a+h*np.arange(2,num-1,2)
    points=np.hstack((p1,p2,a,b))
    mult=np.hstack((np.repeat(4,p1.shape[0]),np.repeat(2,p2.shape[0]),1,1))
    return np.dot(integrand(points),mult)*h/3


A=np.linspace(0,100.,200)

B=lambda x: SIN(x)
C=lambda x: SIN2(x)

t0=time.time()
D=simpson2(B,0,100.,200)
print time.time()-t0
t1=time.time()
E=trapz(C(A))
print time.time()-t1

    t2=time.time()
    F=simpson2(C,0,100.,200)
    print time.time()-t2

结果：

0.000764131546021 sec for the faster method, but when using sympy

7.58171081543e-05 sec for my slower method, but which uses numpy

0.000519037246704 sec for the faster method, when using numpy,

结论：使用 numpy，ditch sympy，（在这种情况下，我较慢的 numpy 方法实际上更快，因为在此示例中，我只在一个 sin 函数上尝试过，而不是在其中的一个 ndarray 上尝试过，但要放弃比较快速方法的 numpy 版本与快速方法的 sympy 版本的时间时，sympy 仍然存在）

【讨论】：

是的，您需要 50 声望才能添加 cmets，但非常感谢您的评论。我只想集成一个返回矩阵/数组的函数，包含 sympy 示例以提供一些背景
但是你的函数“足够温和”吗？（它里面有奇点或其他讨厌的东西吗？或者它只是一些函数，比如一些无法解析解决但仍然“足够温和，可以在数值上积分”的东西的产物？）
那么，您是否在寻找性能提升？或者更优雅的代码？
我建议一个分析解决方案 ;-) 因为你想要整合的东西都是 asin(bx+c)*cos(e*x+f)，或类似的（假设快速浏览后，您只想将您在该 txt 文件中获得的功能集成），（如果可以进行分析，那么就这样做），并摆脱 sympy，因为它很慢。通过 numpy.sin 等更改所有 sympy.sin，这应该会大大提高你的性能
@sgpc，澄清一下，您的意思是您对分析解决方案不感兴趣？

【解决方案4】：

我可能已经找到了一些有趣的方法，但代价是为矩阵 g_symp 定义不同的符号：

import numpy as np
from scipy.integrate import quad
import sympy as sy

@np.vectorize
def vec_lambdify(var, expr, *args, **kw):
    return sy.lambdify(var, expr, *args, **kw)

@np.vectorize
def vec_quad(f, a, b, *args, **kw):
    return quad(f, a, b, *args, **kw)[0]

Y = sy.symbols("y1:11")
x = sy.symbols("x")
mul_x = [y.subs(y,x*(i+1)) for (i,y) in enumerate(Y)]
pow_x = [y.subs(y,x**(i+1)) for (i,y) in enumerate(Y)]

g_sympy = np.array(mul_x + pow_x).reshape((2,10))
X = x*np.ones_like(g_sympy)
G = vec_lambdify(X, g_sympy)
I = vec_quad(G, 0, 100)
print(I)

结果：

[[  5.00000000e+03   1.00000000e+04   1.50000000e+04   2.00000000e+04
    2.50000000e+04   3.00000000e+04   3.50000000e+04   4.00000000e+04
    4.50000000e+04   5.00000000e+04]
[  5.00000000e+03   3.33333333e+05   2.50000000e+07   2.00000000e+09
   1.66666667e+11   1.42857143e+13   1.25000000e+15   1.11111111e+17
   1.00000000e+19   9.09090909e+20]]

并使用 ipython magic%timeit vec_quad(G,0,100) 我得到了

1000 loops, best of 3: 527 µs per loop

我认为这种方法更简洁一些，尽管要使用符号。

【讨论】：

【解决方案5】：

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

import numpy
import quadpy


def f(x):
    return [
        [k * x for k in range(2, 11)],
        [x ** k for k in range(2, 11)],
    ]


sol, err = quadpy.quad(f, 0.0, 100.0, epsabs=numpy.inf, epsrel=1.0e-10)

print(sol)
print(err)

给予

[[1.00000000e+04 1.50000000e+04 2.00000000e+04 2.50000000e+04
  3.00000000e+04 3.50000000e+04 4.00000000e+04 4.50000000e+04
  5.00000000e+04]
 [3.33333333e+05 2.50000000e+07 2.00000000e+09 1.66666667e+11
  1.42857143e+13 1.25000000e+15 1.11111111e+17 1.00000000e+19
  9.09090909e+20]]
[[5.11783704e-16 4.17869644e-16 1.02356741e-15 9.15506521e-16
  8.35739289e-16 1.19125717e-15 2.04713482e-15 1.93005721e-15
  1.83101304e-15]
 [6.69117036e-14 9.26814751e-12 1.05290634e-09 1.12081237e-07
  1.09966583e-05 1.09356156e-03 1.00722052e-01 9.31052614e+00
  9.09545305e+02]]

时间安排：

%timeit quadpy.quad(f, 0.0, 100.0, epsabs=numpy.inf, epsrel=1.0e-10)    
904 µs ± 3.02 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

【讨论】：