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

Question

Я решил определенный интеграл выражения, используя sympy, чтобы получить символическое выражение интеграла. Однако, когда я использую полученное выражение в функции, я не получаю тот же результат, который дает числовая оценка интеграла:

>> from sympy import *
>> x, y, a, b, c, d, k = symbols ('x y a b c d k', positive=True)
>> res = integrate(exp(-k*abs(x-y)), (x, a, b), (y, c, d))

>> res
(-exp(a*k) + exp(b*k))*exp(-b*k)*exp(-k*(a - d))/k**2 - (-exp(a*k) + exp(b*k))*exp(-b*k)*exp(-k*(a - c))/k**2

>> def integral_1(k1, a1, b1, c1, d1):
>>     return (-exp(a1*k1) + exp(b1*k1))*exp(-b1*k1)*exp(-k1*(a1 - d1))/k1**2 - (-exp(a1*k1) + exp(b1*k1))*exp(-b1*k1)*exp(-k1*(a1 - c1))/k1**2

>> integral_1(0.6, 0, 1, 0, 1)
1.0303623235681536

>> integrate(exp(-0.6*abs(x-y)), (x, 0, 1), (y, 0, 1))
0.826731311633480

Почему я получаю такую ​​разницу?

Answer 1

Символический вывод неверен из-за ошибки SymPy. В стороне: лучший способ продемонстрировать это res.subs({k: 0.6, a: 0, b: 1, c: 0, d: 1}), который возвращает 1.03 ..., что явно невозможно для интеграла функции, ограниченной 1 на единичном квадрате.

Проблема в том, что SymPy (начиная с версии 1.3) неправильно обрабатывает abs в этом интеграле. Его вывод res будет таким же, если вы вообще не включите abs. Более простой пример:

>>> integrate(exp(abs(x-y)), (x, a, b))
-exp(a - y) + exp(b - y)

Это не может быть правильно, потому что интеграл должен зависеть от того, как y размещен относительно a и b.

Средство защиты (до устранения ошибки) - применить rewrite(Piecewise) к подынтегральному выражению, заменив abs на кусочную функцию. Например, exp(abs(x-y)).rewrite(Piecewise) возвращает exp(Piecewise((x - y, x - y >= 0), (-x + y, True))), который интегрируется правильно. После

res = integrate(exp(-k*abs(x-y)).rewrite(Piecewise), (x, a, b), (y, c, d))

результат длинный и безобразный, но правильный:

>>> res.subs({k: 0.6, a: 0, b: 1, c: 0, d: 1})
0.826731311633481

Результат будет не таким длинным и безобразным, если вы скажете SymPy, что a >>> res.subs({a < b: True}) -2*Min(d, Max(a, c))/k + 2*Min(d, Max(a, b, c))/k + (-exp(a*k) + exp(b*k))*exp(-b*k)*exp(-k*(a - Min(d, Max(a, c))))/k**2 - (-exp(a*k) + exp(b*k))*exp(-b*k)*exp(-k*(a - Min(c, d)))/k**2 - (exp(a*k) - exp(b*k))*exp(-a*k)*exp(-k*(-a + Min(d, Max(a, b, c))))/k**2 + (exp(a*k) - exp(b*k))*exp(-a*k)*exp(-k*(-a + d))/k**2 - (-k**2*exp(a*k)*exp(k*(-a + Min(d, Max(a, c)))) + k**2*exp(b*k)*exp(-k*(-a + Min(d, Max(a, c)))))*exp(-b*k)/k**4 + (-k**2*exp(a*k)*exp(k*(-a + Min(d, Max(a, b, c)))) + k**2*exp(b*k)*exp(-k*(-a + Min(d, Max(a, b, c)))))*exp(-b*k)/k**4