Вопрос о численном интегрировании и выборе верхнего предела

[alexey007]

Привет, всем!

Мне нужно численно интегрировать следующий интеграл:
$\int_{0}^{\infty} \ln\left(1 - 2 e^{-L\sqrt{M^2+x^2}} + e^{-2L\sqrt{M^2+x^2}}\right) dx $

$L,M$— некоторые константы. Я хочу интегрировать быстро и точно.

Я пока додумался до такого подхода:

1) Интегрировать до какого-то большого значения $x_{\max}$
⁡
2) Разбить интервал интегрирования на $n$ точек, используя формулы для оценки погрешности численного интегрирования. Там нужно будет найти значение второй производной. Пока кажется это не сложно.

У меня возникают трудности с первым пунктом: как найти такое значение $x_{\max}$, чтобы остаток интегрирования был меньше 0.1%. Например так $\frac{\int_{x_{\max}}^{\infty}F(x)dx}{\int_{0}^{\infty}F(x)dx}<0.001$

Можно ли как-то автоматизировать выбор верхнего предела, чтобы эта погрешность удовлетворялась? Мне нужно, чтобы оценка численного интегрирования тоже быстро делалась.
Или, может быть, есть другие методы, которые помогут быстро вычислить данный интеграл с заданной точностью быстро?

Заранее благодарен за любую помощь и советы!

[svv]

Давайте для начала немного упростим интеграл.
Во-первых,
$\ln(1-2e^{-y}+e^{-2y})=\ln (1-e^{-y})^2=2\ln (1-e^{-y})$
Во-вторых, сделаем замену $x=Mt$ (я считаю $M>0$).
Получится
$\int\limits_{0}^{\infty} \ln\left(1 - 2 e^{-L\sqrt{M^2+x^2}} + e^{-2L\sqrt{M^2+x^2}}\right) dx=2M\int\limits_{0}^{\infty} \ln\left(1 - e^{-P\sqrt{1+t^2}}\right) dt,$
где $P=LM$. Теперь интеграл зависит только от одного параметра $P$.
При желании можно переменную интегрирования опять обозначить $x$.

[svv]

Обозначим
$I=\int\limits_{0}^{\infty} \ln\left(1 - e^{-P\sqrt{1+t^2}}\right) dt$
Разложим логарифм в ряд Тейлора:
$\ln\left(1 - e^{-P\sqrt{1+t^2}}\right)=-\sum\limits_{n=1}^\infty\frac{1}{n} e^{-nP\sqrt{1+t^2}}$
Подставим результат в интеграл вместо логарифма и поменяем порядок суммирования и интегрирования:
$I=-\sum\limits_{n=1}^\infty\frac{1}{n} \int\limits_{0}^{\infty}e^{-nP\sqrt{1+t^2}} dt$
Этот интеграл есть в справочнике Градштейна и Рыжика (3.914, с.496), там надо взять $b=0,\gamma=1,\beta=nP$. Он равен $K_1(nP)$, где $K_1$ — модифицированная функция Бесселя второго рода. Получаем
$I=-\sum\limits_{n=1}^\infty\frac{1}{n} K_1(nP)$
Функция $K_1(nP)$ при $n\to+\infty$ стремится к нулю быстрее, чем $e^{-nP}$, поэтому ряд быстро сходится.

Программка на MATLAB, которая подтверждает, что интеграл $I$ найден верно:

Используется синтаксис Matlab M

P = 1.2; % для примера
f = @(t) log(1-exp(-P*sqrt(1+t.^2)));
I = integral(f,0,Inf)

S = 0;
for n=1:50
    S = S-besselk(1, n*P)/n;
end
S

Результат:
I = -0.484665986677563
S = -0.484665986677563
(это всё без множителя $2M$; потом не забыть умножить!)

Ряд лучше интеграла тем, что не нужно контролировать шаг. Недостаток в том, что нужна функция, вычисляющая $K_1$.

[мат-ламер]

alexey007 в сообщении #1642040 писал(а):

У меня возникают трудности с первым пунктом: как найти такое значение $x_{\max}$, чтобы остаток интегрирования был меньше 0.1%.

А асимптотику остатка (той части интеграла, что выше $x_{\max}$) вы пробовали найти?

[alexey007]

svv, спасибо большое, очень круто!
Буду использовать ряд, с помощью него действительно легче получить интеграл с заданной суммой.
Единственное, что данный интеграл это один из частных случаев, по которому я хотел оценивать предел и шаг интегрирования для более общих случаев. Но, даже для этого частного случая смогу вычислять интеграл быстрее.

мат-ламер, да пробовал. Т.е. теперь если я знаю полный интеграл посчитанный с помощью ряда, предложенного svv
Могу теперь еще посчитать интеграл от асимптотики при больших $t$:
$\int\limits_{t_{\max}}^{\infty}\ln(1 - e^{-P\sqrt{1+t^2}})\approx - \int\limits_{t_{\max}}^{\infty}e^{-P\sqrt{1+t^2}}$
и тут тупик... Не могу найти первооразную... Может если только дальше продолжить искать асимптотику и не отбросить 1 под корнем и свести уже вот к такому:
$\approx - \int\limits_{t_{\max}}^{\infty}e^{-Pt}$ Это будет не совсем уж?
т.е. если теперь я запишу соотношения для поиска $t_{\max}$
$\frac{\int\limits_{t_{\max}}^{\infty}e^{-Pt}}{\sum\limits_{n=1}^\infty\frac{1}{n} K_1(nP)}=0.001$
Решив данное данное уравнение, я же буду получать предел более консервативно $x_{\max}$, т.е. такой предел $x_{\max}$, который будет дальше находиться, чем если бы все считать без асимптотик и я же буду получать заведомо больше точность чем требуемая, то тогда хорошо. Ведь это будте так?
но если при таком упращении можно натнуться наоброт на большое несоответсие, то не очень тогда)))

[мат-ламер]

alexey007 в сообщении #1642096 писал(а):

мат-ламер, да пробовал. Т.е. теперь если я знаю полный интеграл посчитанный с помощью ряда, предложенного svv

Вот именно. С помощью ряда, предложенного svv . Там подынтегральная функция крайне быстро стремиться к нулю.

А в чём вообще смысл моего предложения? Вот вы писали:

alexey007 в сообщении #1642040 писал(а):

Я хочу интегрировать быстро и точно.

То есть интеграл будет считаться многократно. Допусти, интеграл выше $x_{max}$ считается у нас крайне быстро. А вот интеграл ниже $x_{max}$ можно очень быстро подсчитать каким-нибудь простым численным способом. Возможно мы тут получим ускорение по сравнению с методом, предложенным svv . Но тут надо провести эксперименты. Посмотреть, что получится.

[мат-ламер]

Тут есть интересный метод - адаптивного разбиения. Сам программировал. И он достаточно прост в реализации. Хотя требует знания, как программировать рекурсию. Для начала выбираем достаточно большой отрезок интегрирования. И к нему применяем метод Симпсона. Затем делим отрезок пополам и к двум полученным подотрезкам применяем метод Симпсона. Затем рекурсивно делим только те пары подотрезков, интеграл по которым достаточно сильно отличается от интеграла по объемлищему их одному отрезку. Так продолжаем, пока есть возможность деления.

[alexey007]

мат-ламер в сообщении #1642197 писал(а):

Тут есть интересный метод - адаптивного разбиения. Сам программировал. И он достаточно прост в реализации. Хотя требует знания, как программировать рекурсию. Для начала выбираем достаточно большой отрезок интегрирования. И к нему применяем метод Симпсона. Затем делим отрезок пополам и к двум полученным подотрезкам применяем метод Симпсона. Затем рекурсивно делим только те пары подотрезков, интеграл по которым достаточно сильно отличается от интеграла по объемлищему их одному отрезку. Так продолжаем, пока есть возможность деления.

Спасибо! А есть статья или название метода. Где можно по подробнее почитать?

[Markiyan Hirnyk]

Не открывая Америку, а применяя Математику 14:

ClearAll[M, L, x]; L = E/10; M = Pi; 
NIntegrate[ Log[1 - 2*Exp[-L*Sqrt[M^2 + x^2]] + Exp[-2*L*Sqrt[M^2 + x^2]]], {x,  0, Infinity}] // Timing

{0.015625, -5.74387}

NIntegrate[ Log[1 - 2*Exp[-L*Sqrt[M^2 + x^2]] + Exp[-2*L*Sqrt[M^2 + x^2]]], {x,    0, Infinity}, PrecisionGoal -> 20, AccuracyGoal -> 20,  WorkingPrecision -> 20] // Timing

{0.140625, -5.7438656740575619202}

[мат-ламер]

alexey007 в сообщении #1642199 писал(а):

Спасибо! А есть статья или название метода. Где можно по подробнее почитать?

Честно, не помню. Давно было. Видал в некоторых учебниках на английском. Но подробное описание метода не нужно. Вы легко сами по интуиции восстановите подробности. Если в классическом методе Симпсона у нас сетка разбиения равномерна, то тут получается сетка с неравномерным разбиением. Причём она сильнее разбита именно там, где нам нужно. А там, где функция стремится к нулю, останется крупная сетка и разбиений больше не будет. Ещё есть маленький необязательный нюанс. На следующем шаге разбиения можно использовать значения подынтегральной функции, вычисленные на предыдущем шаге.

[alexey007]

мат-ламер в сообщении #1642197 писал(а):

Тут есть интересный метод - адаптивного разбиения.

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

В итоге я написал библиотеку. Как только закончу тестирование, скину ссылку на репозиторий — возможно, кому-то она пригодится.

Прошу предложить интересные примеры двойных интегралов, желательно с бесконечными пределами. Хочу подготовить статью о том, как удалось обобщить метод Ромберга для интегралов от двух переменных и насколько быстро этот метод работает по сравнению со стандартными библиотеками. В SciPy нет функций для интегрирования методом Ромберга от двух переменных, а для интегралов от одной переменной я переписал функцию, сделав её более оптимальной, что также ускорило вычисления.

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

[ozheredov]

alexey007 в сообщении #1643399 писал(а):

Хочу подготовить статью о том, как удалось обобщить метод Ромберга для интегралов от двух переменных

Решил восполнить свой пробел в аброзовании и возгуглить метод Ромберга. Что нашел (например):

https://www.mathnet.ru/links/d154ca9978 ... vmmf35.pdf

Abstract писал(а):

Известный метод Ромберга, применяемый для повышения точности вычисления одномерных
интегралов, обобщен на случай кратных интегралов, если для их вычисления используется
произведение составных квадратурных формул. При определенных условиях коэффициенты
формулы Ромберга оказываются независимыми от кратности интеграла, что позволяет использовать простой алгоритм вычислений, развитый для одномерных интегралов. Даются
примеры расчетов по методу Ромберга для интегралов кратности от двух до шести и производится сравнение с некоторыми другими методами.

Оно или нет?