与 scipy 的累积辛普森集成

Question

我有一些代码使用 scipy.integration.cumtrapz 来计算采样信号的反导数。我想使用辛普森规则而不是梯形。然而 scipy.integration.simps 似乎没有累积对应物......我错过了什么吗？有没有一种简单的方法来获得与“scipy.integration.simps”的累积集成？

Answer 1

你总是可以自己写：

def cumsimp(func,a,b,num):
    #Integrate func from a to b using num intervals.

    num*=2
    a=float(a)
    b=float(b)
    h=(b-a)/num

    output=4*func(a+h*np.arange(1,num,2))
    tmp=func(a+h*np.arange(2,num-1,2))
    output[1:]+=tmp
    output[:-1]+=tmp 
    output[0]+=func(a)
    output[-1]+=func(b)
    return np.cumsum(output*h/3)

def integ1(x):
    return x

def integ2(x):
    return x**2

def integ0(x):
    return np.ones(np.asarray(x).shape)*5

先看一个常数函数的求和和导数。

print cumsimp(integ0,0,10,5)
[ 10.  20.  30.  40.  50.]

print np.diff(cumsimp(integ0,0,10,5))
[ 10.  10.  10.  10.]

现在检查一些琐碎的例子：

print cumsimp(integ1,0,10,5)
[  2.   8.  18.  32.  50.]

print cumsimp(integ2,0,10,5)
[   2.66666667   21.33333333   72.          170.66666667  333.33333333]

在这里显式编写您的被积函数要容易得多，然后在这种情况下重现 scipy 的 simpson 规则函数。当提供单个数组时，选择间隔将很难做到，你是不是：

• 使用所有其他值作为辛普森规则的边缘，其余值作为中心？
• 使用数组作为边并插入中心值？

还有一些选项可用于您希望如何对间隔求和。这些并发症可能是它没有用 scipy 编码的原因。

Answer 2

您的问题很久以前就已经回答了，但我最近遇到了同样的问题。我写了一些函数来计算等距点的累积积分；代码可以在GitHub 上找到。插值多项式的顺序范围从 1（梯形规则）到 7。正如 Daniel 在上一个答案中指出的那样，必须对如何对区间求和做出一些选择，尤其是在边界处；因此，根据您使用的软件包，结果可能会略有不同。另请注意，对于高阶多项式，数值积分可能会受到龙格现象（意外振荡）的影响。

这是一个例子：

import numpy as np
from scipy import integrate as sp_integrate
from gradiompy import integrate as gp_integrate

# Definition of the function (polynomial of degree 7)
x = np.linspace(-3,3,num=15)
dx = x[1]-x[0]
y = 8*x + 3*x**2 + x**3 - 2*x**5 + x**6 - 1/5*x**7
y_int = 4*x**2 + x**3 + 1/4*x**4 - 1/3*x**6 + 1/7*x**7 - 1/40*x**8

# Cumulative integral using scipy
y_int_trapz = y_int [0] + sp_integrate.cumulative_trapezoid(y,dx=dx,initial=0)
print('Integration error using scipy.integrate:')
print('  trapezoid = %9.5f' % np.linalg.norm(y_int_trapz-y_int))

# Cumulative integral using gradiompy
y_int_trapz = gp_integrate.cumulative_trapezoid(y,dx=dx,initial=y_int[0])
y_int_simps = gp_integrate.cumulative_simpson(y,dx=dx,initial=y_int[0])
print('\nIntegration error using gradiompy.integrate:')
print('  trapezoid = %9.5f' % np.linalg.norm(y_int_trapz-y_int))
print('  simpson   = %9.5f' % np.linalg.norm(y_int_simps-y_int))

# Higher order cumulative integrals
for order in range(5,8,2):
    y_int_composite = gp_integrate.cumulative_composite(y,dx,order=order,initial=y_int[0])
    print('  order %i   = %9.5f' % (order,np.linalg.norm(y_int_composite-y_int)))
    
# Display the values of the cumulative integral
print('\nCumulative integral (with initial offset):\n',y_int_composite)

你应该得到以下结果：

'''
Integration error using scipy.integrate:
  trapezoid = 176.10502

Integration error using gradiompy.integrate:
  trapezoid = 176.10502
  simpson   =   2.52551
  order 5   =   0.48758
  order 7   =   0.00000

Cumulative integral (with initial offset):
 [-6.90203571e+02 -2.29979407e+02 -5.92267425e+01 -7.66415188e+00
  2.64794452e+00  2.25594840e+00  6.61937372e-01  1.14797061e-13
  8.20130517e-01  3.61254267e+00  8.55804341e+00  1.48428883e+01
  1.97293221e+01  1.64257877e+01 -1.13464286e+01]
'''

Answer 3

我会选择 Daniel 的解决方案。但是，如果您正在积分的功能本身会受到波动的影响，您需要小心。 Simpson's 要求函数表现良好（在这种情况下，意味着函数是连续的）。

有一些技术可以使表现中等的函数看起来比实际表现更好（实际上是函数的近似形式），但在这种情况下，您必须确保该函数“充分”近似于您的函数。在这种情况下，您可能会使间隔可能不均匀以处理问题。

一个例子可能是考虑一个场的流动，在较长的时间尺度上，它是由一个表现良好的函数近似的，但在较短的时间内，它的密度会受到有限的随机波动。