scipy.integrate.quad не срабатывает (иногда), когда интегрируемая функция также является интегральной - PullRequest
Новая
       71

scipy.integrate.quad не срабатывает (иногда), когда интегрируемая функция также является интегральной

0 голосов
/ 25 октября 2018

Я использую sqipy.integrate.quad для вычисления двойного интеграла.В основном я пытаюсь вычислить интеграл по exp [-mu_wx_par], где mu_wx_par также является интегралом.

Мой код в основном работает.Тем не менее, для некоторых значений он терпит неудачу, то есть он возвращает неправильные значения.

import numpy as np
from scipy import integrate


def mu_wx_par(x, year, par):
    """ First function to be integrated """
    m = max(par['alfa'], 0) + par['beta'] * 10 ** (par['gamma'] * x)
    w = np.minimum(par['frem_a'] + par['frem_b'] * x + par['frem_c'] * x**2, 0)
    return m * (1 + w/100)**(year - par['innf_aar'])


def tpx_wx(x, t, par, year):
    """ Second function to be integrated (which contains an integral itself)"""
    result, _ = integrate.quad(lambda s: mu_wx_par(x + s, year + s, par), 0, t)
    return np.exp(-result)


def est_lifetime(x, year, par):
    """ Integral of second function. """
    result, _ = integrate.quad(lambda s: tpx_wx(x, s, par, year), 0, 125 - x)

    return result


# Test variables
par = {'alfa': 0.00019244401470947973,
       'beta': 2.420260552210541e-06,
       'gamma': 0.0525500987420195,
       'frem_a': 0.3244611019518985,
       'frem_b': -0.12382978382606026,
       'frem_c': 0.0011901237463116591,
       'innf_aar': 2018
       }


year = 2018
estimate_42 = est_lifetime(42, year, par)
estimate_43 = est_lifetime(43, year, par)
rough_estimate_42 = sum([tpx_wx(42, t, par, year) for t in range(0, 100)])

print(estimate_42)
print(estimate_43)
print(rough_estimate_42)
3.1184634065887544
46.25925442287578
47.86323490659588

Значение estimate_42 неверно.Это должно быть примерно то же значение, что и rough_estimate_42.Обратите внимание, что estimate_43 выглядит хорошо.Что здесь происходит?

Я использую scipy v1.1.0 и numpy v1.15.1 и Windows.

Было высказано предположение, что функция близка к нулю почти везде, как в этомpost scipy integrate.quad возвращает неправильное значение .Это не так, поскольку простой график tpx_wx для x=42 от a=0 до b=125-42 ясно показывает 

from matplotlib import pyplot as plt

plt.plot(range(125-42), [tpx_wx(42, t, par, year) for t in range(0, 125-42)])
plt.show()

imagetpx_wx over integration range">

1 Ответ

0 голосов
/ 28 октября 2018

Похоже, это известная проблема , связанная с тем, что некоторый код Fortran, стоящий за quad, скомпилирован для Windows: рекурсивный вызов этого кода в некоторых случаях может привести к сбою.См. Также Большая ошибка интеграции с integrate.nquad .

За исключением перекомпиляции SciPy с лучшими флагами, кажется, что следует избегать вложенной интеграции с quad в Windows.Одним из способов решения этой проблемы является использование метода romberg для одного из этапов интеграции.Замена quad в est_lifetime на 

integrate.romberg(lambda s: tpx_wx(x, s, par, year), 0, 125 - x, divmax=20)

приводит к 47.3631754795 для estimate_42, что соответствует тому, что quad возвращает в Linux.

Одним из способов визуализации процесса интеграции является объявление глобального списка eval_points и вставка eval_points.append(t) в tpx_wx.С той же версией SciPy (0.19.1 в этом тесте) результаты plt.plot(eval_points, '.') выглядят иначе.

В Linux:

linux

В Windows:

Windows

Повторное деление пополам окрестности хитрой точки около 60 преждевременно завершается в Windows, и кажется, что полученный результат является неким частным интегралом по подынтервалу.

...