Использование квадроцикла SciPy для получения главного значения интеграла путем интегрирования чуть ниже и чуть выше особой точки

Я пытаюсь вычислить главное значение интеграла (по s) от 1/((s - q02)*(s - q2)) на [Ecut, inf] с q02 ‹ Ecut ‹ q2. Делая основное значение вручную (или Mathematica), можно получить общий результат

ln((q2-Ecut)/(Ecut-q02)) / (q02-q2)

В конкретном примере ниже это дает результат -1,58637*10^-11. Также можно получить тот же результат, разделив интеграл на два, проинтегрировав до q2 - eps, а затем начав с q2 + eps, а затем сложив два результата (расхождения должны сокращаться). Принимая eps все меньше и меньше, можно восстановить результат выше. Когда я реализую это в scipi с помощью quad, мой результат сходится к неправильному результату 6.04685e-11, как я показываю на графике eps против интегрального результата, который я включаю.
Почему quad делает это? даже если у меня eps = 0, это дает мне неправильный результат, когда я ожидал, что он выдаст мне ошибку, когда все взорвется...

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import quad


q02  = 485124412.
Ecut = 17909665929.
q2   = 90000000000.

def integrand(s):
    return 1/((s - q02)*(s - q2))

xx=[1.,0.1,0.01,0.001,0.0001,0.00001,0.000001,0.0000001,0.00000001,
    0.000000001,0.0000000001,0.00000000001,0.]

integral = [0*y for y in xx]
i=0
for eps in xx:

    ans1,err = quad(integrand, Ecut, q2 -eps )
    ans2,err= quad(integrand, q2 + eps, np.inf)

    integral[i] = ans1 + ans2
    i=i+1

plt.semilogx(xx,integral,marker='.')
plt.show()

eps против интегрального результата


scipy numerical-integration quad
person Rudyard    schedule 10.04.2017    source источник

Ответы (1)


arrow_upward
1
arrow_downward

Также можно получить тот же результат, разделив интеграл на два, проинтегрировав до q2 - eps, а затем начав с q2 + eps, а затем сложив два результата.

Только если расчеты были абсолютно точными. В численной практике то, что вы описали, в основном худшее, что можно было сделать. Вы получаете два больших интеграла противоположных знаков, которые при сложении почти компенсируют друг друга; то, что осталось, больше связано с ошибками интегрирования, чем с действительным значением интеграла.

Я заметил, что вы проигнорировали значения ошибок err в своем скрипте, даже не распечатав их. Плохая идея: они имеют размер 1e-10, что уже говорит вам о том, что окончательный результат с "что-то e-11" является мусором.

Вопрос вычислительной науки Числовая интеграция основных значений - как у Гильберта решает эту проблему. Один из подходов, который они указывают, состоит в том, чтобы сложить значения подынтегральной функции в точках, симметричных относительно особенности, перед попыткой ее интегрирования. Для этого необходимо взять интеграл по симметричному интервалу с центром в сингулярности q2 (то есть от Ecut до 2*q2-Ecut), а затем добавить вклад интеграла от 2*q2-Ecut к бесконечность. Это разделение в любом случае имеет смысл, потому что quad трактует бесконечные пределы совершенно по-другому (используя интегрирование Фурье), что является еще одной вещью, которая повлияет на способ устранения сингулярности.

Таким образом, реализация этого подхода будет

ans1, err = quad(lambda s: integrand(s) + integrand(2*q2-s), Ecut, q2)
ans2, err = quad(integrand, 2*q2-Ecut, np.inf)

Эпс не нужен. Тем не менее, результат по-прежнему отсутствует: это около -2.5e-11. Оказывается, виноват второй интеграл. К сожалению, интегральный подход Фурье не кажется здесь эффективным (или я не нашел способ заставить его работать). Оказывается, задание большого, но конечного значения в качестве верхнего предела приводит к лучшему результату, особенно если также используется опция epsabs, т.е. epsabs=1e-20.

А еще лучше прочтите документацию. quad и обратите внимание, что он напрямую поддерживает интегралы с весом Коши 1/(s-q2), выбирая для них соответствующий численный метод. Это по-прежнему требует конечного верхнего предела и небольшого значения epsabs, но результат довольно точен:

quad(lambda s: 1/(s - q02), Ecut, 1e9*q2, weight='cauchy', wvar=q2, epsabs=1e-20)

возвращает -1.5863735715967363e-11 по сравнению с точным значением -1.5863735704856253e-11. Обратите внимание, что множитель 1/(s-q2) не появляется в подынтегральном выражении выше, будучи отнесенным к вариантам веса.

person Community    schedule 11.04.2017