scipy dblquad дает неверный результат в простом двойном интеграле - PullRequest
Новая
       8

scipy dblquad дает неверный результат в простом двойном интеграле

0 голосов
/ 12 ноября 2018

Я пытаюсь вычислить простой двойной определенный интеграл в Python: функция Max (0, (4-12x) + (6-12y)) в квадрате [0,1] x [0,1].

Мы можем сделать это с Mathematica и получить точный результат: 

Integrate[Max[0, 4-12*u1 + 6-12*u2], {u1, 0, 1}, {u2, 0,1}] = 125/108.

С помощью простого моделирования Монте-Карло я могу подтвердить этот результат.Однако, используя scipy.integrate.dblquad, я получаю значение 0,0005772072907971 с ошибкой 0,0000000000031299 

from scipy.integrate import dblquad

def integ(u1, u2):
    return max(0, (4 - 12*u1) + (6 - 12*u2))
sol_int, err = dblquad(integ, 0, 1, lambda _:0, lambda _:1, epsabs=1E-12, epsrel=1E-12)
print("dblquad: %0.16f. Error: %0.16f" % (sol_int, err) )

Согласна, что функция не выводима, но она непрерывна, и я не вижу причин, по которым этот конкретный интеграл может быть проблематичным.

Я подумал, что, возможно, dblquad имеет аргумент 'options', в котором я могу попробовать разные численные методы, но я ничего подобного не нашел.

Итак, что я делаю не так?

1 Ответ

0 голосов
/ 13 ноября 2018

попробуйте разные числовые методы

Это то, что я хотел бы предложить, учитывая проблему, которая повторяется quad в Windows.Изменив его на явный двухэтапный процесс, вы можете заменить один из quad другим методом, romberg кажется мне лучшей альтернативой.

from scipy.integrate import quad, romberg

def integ(u1, u2):
    return max(0, (4 - 12*u1) + (6 - 12*u2))

sol_int = romberg(lambda u1: quad(lambda u2: integ(u1, u2), 0, 1)[0], 0, 1)    
print("romberg-quad: %0.16f " % sol_int)

Это напечатает 1.1574073959987758 на моем компьютере, и, надеюсь, вы получите то же самое.

...