Python + Scipy + Integration : 스파이크가있는 함수의 정밀도 오류 처리

Question

scipy.integrate.quad를 사용하여 매우 큰 범위 (0..10,000)에 함수를 통합하려고합니다. 이 함수는 대부분의 범위에서 0이지만 매우 작은 범위 (예 : 1,602..1,618)의 스파이크를가집니다.

통합 할 때 나는 출력이 양의 값이 될 것으로 기대하지만 쿼드의 추측 알고리즘은 혼란스러워서 제로를 출력한다고 생각합니다. 제가 알고 싶은 것은,이를 극복 할 수있는 방법이 있습니까 (예 : 다른 알고리즘, 다른 매개 변수 등을 사용하는 것)? 나는 보통 스파이크가 어디로 가게 될지 알지 못한다. 그래서 나는 누군가가 그것을하는 법에 대해 좋은 아이디어가 없다면, 통합 범위를 나누고 그 부분들을 합산 할 수 없다.

감사합니다.

샘플 출력 : 당신은 같은 integrate.romberg() 방법으로 다른 통합 방법을 시도 할 수도

>>>scipy.integrate.quad(weighted_ftag_2, 0, 10000) 
(0.0, 0.0) 
>>>scipy.integrate.quad(weighted_ftag_2, 0, 1602) 
(0.0, 0.0) 
>>>scipy.integrate.quad(weighted_ftag_2, 1602, 1618) 
(3.2710994652983256, 3.6297354011338712e-014) 
>>>scipy.integrate.quad(weighted_ftag_2, 1618, 10000) 
(0.0, 0.0) 

Answer 1

.

weighted_ftag_2(x_samples).argmax()과 같이 함수가 큰 지점의 위치를 ​​얻은 다음 일부 휴리스틱 스를 사용하여 함수의 최대 값 (예 : x_samples[….argmax()])에서 통합 간격을 줄일 수 있습니다. 샘플링 된 가로 좌표 목록 (x_samples) : 함수가 최대 인 지역에 항상 포인트가 있어야합니다.

더 일반적으로 통합 할 함수에 대한 모든 특정 정보는 좋은 값을 얻을 수 있도록 도와줍니다. 그 통합에 대한 합리적인 분할과 함께 (Scipy가 제공 한 many methods 중 하나) 함수에 대해 잘 작동하는 방법을 결합 할 것입니다. 위에서 제안 된 선을 따라).

Answer 2

각 정수 범위 [x, x + 1], 에 대한 함수 f()를 계산하고 예를 들어. romb(), EOL에서 알 수 있듯이> 0 :

from __future__ import division 
import numpy as np 
from scipy.integrate import romb 

def romb_non0(f, a=0, b=10000, nromb=2**6+1, verbose=1): 
    """ sum romb() over the [x, x+1) where f != 0 """ 
    sum_romb = 0 
    for x in xrange(a, b): 
     y = f(np.arange(x, x+1, 1./nromb)) 
     if y.any(): 
      r = romb(y, 1./nromb) 
      sum_romb += r 
      if verbose: 
       print "info romb_non0: %d %.3g" % (x, r) # , y 
    return sum_romb 

#........................................................................... 
if __name__ == "__main__": 
    np.set_printoptions(2, threshold=100, suppress=True) # .2f 

    def f(x): 
     return x if (10 <= x).all() and (x <= 12).all() \ 
      else np.zeros_like(x) 

    romb_non0(f, verbose=1)