scipy.integrate.quad точность по большим числам

Я считаю, что проблема связана с тем, что np.exp(-x) быстро становится очень маленьким с ростом x , что приводит к оценке как нулю из-за ограниченной числовой точности. Например, даже для x , np.exp(-x) x=10**2* , np.exp(-x) оценивается до 3.72007597602e-44 , тогда как значения x порядка 10**3 или выше приводят к 0 .

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

Один из способов избежать этой проблемы состоит в том, чтобы ограничить верхнюю границу интеграции до значения, выше которого числовая функция становится очень малой (и, следовательно, незначительно влияет на интегральное значение). Из вашего кода snipet значение 10**4 представляется хорошим выбором, однако значение 10**2 также приводит к точной оценке интеграла.

Другим подходом к решению проблем с числовой точностью является использование модуля, который выполняет вычисления в арифметике произвольной точности, такой как mpmath . Например, для x=10**5 mpmath вычисляет exp(-x) следующим образом (с использованием экспоненциальной функции нативного mpmath )

 import mpmath as mp print(mp.exp(-10**5)) 

3.56294956530937e-43430

Обратите внимание, насколько мало это значение. При стандартной аппаратной числовой точности (используется numpy ) это значение становится равным 0 .

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

 import mpmath as mp print(mp.quad(lambda x : .5 * mp.exp(-.5 * x), [0, mp.inf])) print(mp.quad(lambda x : .5 * mp.exp(-.5 * x), [0, 10**13])) print(mp.quad(lambda x : .5 * mp.exp(-.5 * x), [0, 10**8])) print(mp.quad(lambda x : .5 * mp.exp(-.5 * x), [0, 10**5])) 

 1.0 0.999999650469474 0.999999999996516 0.999999999999997 

Мы можем также получить еще более точные оценки, увеличивая точность до, скажем, 50 десятичных точек (от 15 которая является стандартной точностью)

 mp.mp.dps = 50; print(mp.quad(lambda x : .5 * mp.exp(-.5 * x), [0, mp.inf])) print(mp.quad(lambda x : .5 * mp.exp(-.5 * x), [0, 10**13])) print(mp.quad(lambda x : .5 * mp.exp(-.5 * x), [0, 10**8])) print(mp.quad(lambda x : .5 * mp.exp(-.5 * x), [0, 10**5])) 

 1.0 0.99999999999999999999999999999999999999999829880262 0.99999999999999999999999999999999999999999999997463 0.99999999999999999999999999999999999999999999999998 

В общем, стоимость получения этой точности – это увеличенное время вычисления.

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

Используйте аргумент points чтобы сообщить алгоритму, где поддержка вашей функции примерно такова:

 import numpy as np from scipy.integrate import quad def g(x): return .5 * np.exp(-.5 * x) print quad(g, a=0., b=10**3, points=[1, 100]) print quad(g, a=0., b=10**6, points=[1, 100]) print quad(g, a=0., b=10**9, points=[1, 100]) print quad(g, a=0., b=10**12, points=[1, 100])