logo标签列表

使用Scipy进行累积辛普森积分

numpyscipynumerical-integration
5

我有一些代码，使用scipy.integration.cumtrapz来计算采样信号的反导数。我想使用辛普森规则(Simpson's rule)而不是梯形(Trapezoid)。然而，scipy.integration.simps似乎没有累积的对应函数... 我错过了什么吗？有没有简单的方法可以使用"scipy.integration.simps"实现累积积分？

3个回答
4

您始终可以自己编写：

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's rule函数要容易得多。当提供一个单一的数组时，选择间隔会很难处理，您可以选择：

  • 使用每个值的奇偶性作为simpson's rule的边缘值，其余值作为中心值？
  • 使用数组作为边界并插值计算中心值？

还有一些选项可以控制如何对间隔进行求和。这些复杂性可能是它没有在scipy中编码的原因。

0

你的问题早已得到解答，但我最近遇到了同样的问题。我编写了一些函数来计算等间距点的累积积分；代码可以在GitHub上找到。插值多项式的阶数从1（梯形法则）到7不等。正如丹尼尔在之前的回答中指出的那样，在如何对间隔进行求和方面必须做出一些选择，特别是在边界处；因此，使用的软件包可能会导致结果略有不同。还要注意，数值积分可能会因为高阶多项式的龙格现象（意外振荡）而受到影响。

这里是一个例子：

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]
'''
0
我会选择Daniel的解决方案。但是，如果您要集成的函数本身受到波动的影响，那么您需要小心。辛普森法则要求函数表现良好（在这种情况下，指连续函数）。
有一些技术可以使一个中等程度上不良行为的函数看起来比它实际上更良好（实际上是对函数的逼近形式），但在这种情况下，您必须确保该函数“足够”逼近您的函数。在这种情况下，您可能需要使区间非均匀以处理问题。
例如，在考虑一个场的流动时，长时间尺度上可以近似为良好的函数，但在短时间内却受到密度有限的随机波动的影响。
这并没有回答问题。一旦您拥有足够的声望，您将能够评论任何帖子；相反，提供不需要询问者澄清的答案。- 来自审核 - Trenton McKinney
网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接