Интеграция Python + Scipy +: работа с ошибками точности в функциях с пиками - PullRequest
Новая
       45

Интеграция Python + Scipy +: работа с ошибками точности в функциях с пиками

2 голосов
/ 06 июля 2010

Я пытаюсь использовать scipy.integrate.quad для интеграции функции в очень большом диапазоне (0..10,000).Функция равна нулю в большей части своего диапазона, но имеет всплеск в очень маленьком диапазоне (например, 1 602,11818).

При интегрировании я ожидаю, что результат будет положительным, но я предполагаю, что каким-то образом квадалгоритм гадания запутывается и выдает ноль.Я хотел бы знать, есть ли способ преодолеть это (например, с помощью другого алгоритма, какого-то другого параметра и т. Д.)?Я обычно не знаю, где будет всплеск, поэтому я не могу просто разделить диапазон интеграции и суммировать части (если у кого-то нет хорошей идеи о том, как это сделать).

Спасибо!

Пример вывода: 

>>>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)

Ответы [ 2 ]

2 голосов
/ 06 июля 2010

Возможно, вы захотите попробовать другие методы интеграции, такие как метод integrate.romberg().

В качестве альтернативы, вы можете получить местоположение точки, где ваша функция велика, с помощью weighted_ftag_2(x_samples).argmax(), а затем использоватьнекоторая эвристика, чтобы сократить интервал интеграции вокруг максимума вашей функции (который находится в x_samples[….argmax()]. Вы должны указать список отобранных абсцисс (x_samples) для вашей проблемы: он всегда должен содержать точки, которые находятся в той области, гдеваша функция максимальна.

В целом, любая конкретная информация об интегрируемой функции может помочь вам получить хорошее значение для ее интеграла. Я бы комбинировал метод, который хорошо работает для вашей функции (один из много методов , предложенных Сципи) с разумным разбиением интервала интегрирования (например, по линиям, предложенным выше).

1 голос
/ 07 июля 2010

Как насчет вычисления вашей функции f () в каждом целочисленном диапазоне [x, x + 1) и сложения, например, 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 )
...