Вычисление интегралов от полиномов Лагерра.

[nestoklon]

В ходе решения некоторой задачи возникла задача посчитать интегралы вида

$$P_N^M(\gamma)=\int_0^{\infty} e^{-\gamma x - \frac12 x^2} x^M L_{N-M}^{M} (x^2) dx,$$

где $L_N^M$ -- полиномы Лагерра. Причём хочется делать это точно и по мере возможности быстро.
Начал с простого -- для небольших $N$ численно интегрирую, для больших заменяю полиномы на функцию Бесселя и пишу ответ. Хотелось бы сделать что-то получше.

В принципе, есть вариант выразить эти интегралы через

$$R_N(\gamma)=\int_0^{\infty} e^{-\gamma x - \frac12 x^2} x^N dx,$$

но тут опять засада. Для этих интегралов можно получить рекуррентное соотношение

$$R_N(\gamma) = \sqrt2 S_N(\gamma) + \left[ \sqrt{\frac{\pi}2} e^{\frac12 \gamma^2} \mathrm{erfc}\left( \frac{\gamma}{\sqrt2} \right) \right] T_N(\gamma),$$

$$S_{N+1}(\gamma) = T_N(\gamma) - \frac{\partial S_N(\gamma)}{\partial \gamma} ,\;\;\;\; T_{N+1}(\gamma) = - \gamma T_N(\gamma) - \frac{\partial T_N(\gamma)}{\partial \gamma}$$

Плюс $S_0=0$ $T_0=1$.
Если попытаться посмотреть на несколько первых $N$, оказывается что так считать точно и быстро не получится, потому что по отдельности оба слагаемых довольно велики, то есть для вычисления $R_N$ с разумной точностью придётся сделать что-то весьма нетривиальное.

Собственно, в этом и вопрос. Что такое нетривиальное можно сделать чтобы достаточно точно посчитать $R_N$ или сразу $P_N^M$? Может, кто-то сталкивался с подобной задачей?

Да, сразу оговорюсь. Вопрос задаю в Computer Science потому, что в принципе $R_N$ можно выразить через конечную сумму функций параболического цилиндра, т.е. получить замкнутый ответ. Вот только посчитать для аргументов больше 1 не получится, так как в какой-то момент придётся умножать большое значение экспоненты на маленькое значение функции параболического цилиндра, что приводит к ещё большим проблемам, несмотря на "аналитическое" решение.
Вопрос именно в том, как это можно посчитать.

Заранее спасибо.

[ГАЗ-67]

nestoklon в сообщении #339833 писал(а):

В принципе, есть вариант выразить эти интегралы через

$$R_N(\gamma)=\int_0^{\infty} e^{-\gamma x - \frac12 x^2} x^N dx,$$

В показателе экспоненты выделите полный квадрат и сделайте соответствующую замену переменной .
Далее рассматривайте случаи чётных и нечётных $N$.

[nestoklon]

ГАЗ-67 в сообщении #340320 писал(а):

В показателе экспоненты выделите полный квадрат и сделайте соответствующую замену переменной .
Далее рассматривайте случаи чётных и нечётных $N$.

Вопрос не в том, как получить аналитический ответ. У меня есть пара вариантов рекурсивного вычисления этой функции, например такой (для простоты привожу ответы для $N=0$):
$P_0^0=\sqrt{\frac{\pi}2}e^{\frac{x^2}2}\operatorname{erfc}\left( \frac{x}{\sqrt2}\right)$
$P_1^0=x-x^2 P_0^0$
$P_N^0=\frac{(2N-1)P_{N-1}^0-\frac{\partial^2}{\partial x^2} P_{N-1}^0 - (N-1)P_{N-2} }N$
Можно и без дифференцирования -- если сначала найти функцию, потом считать, никакой разницы.

Проблема в том, что результат -- выражение вида $P_N^0(x)=S_N(x)+T_N(x)P_0^0(x)$, где $S$ и $T$ -- полиномы степени $N$. Для $N$ начная где-то с 10 и $x$ начиная где-то с 8 это выражение становится численно неустойчиво -- нельзя в машинной точности вычислять величину как разницу двух больших чисел. Да ещё и $\operatorname{erfc}$ не самая хорошая функция.

Если интересно посмотреть на проблему, привожу маленькую програмку на питоне которая рисует картинку, которая показывает что с этой функцией не так. Обратите внимание, что сначала аналитически вычисляется функция, потом в неё просто подставляется число. Делать что-то надо с этой функцией -- представлять $\operatorname{erfc}$ в виде ряда и как-то тасовать порядок суммирования.

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

#!/usr/bin/python

import operator
import numpy 
import matplotlib.pyplot as plt
from scipy import special
from scipy import integrate
import sympy

def factorial_ratio(M,N):
   if M==0: return 1
   if M>=N: return 1
   return reduce(operator.mul,xrange(N+1-M,N+1))
def func_in_P(M,N,x):
   if M==0:
      return lambda xi: numpy.exp(-x*xi-xi*xi/2)*special.laguerre(N)(xi*xi)
   else:
      return lambda xi: 1.0/math.sqrt(float(factorial_ratio(M,N)))*numpy.exp(-x*xi-0.5*xi*xi)*numpy.power(xi,M)*special.genlaguerre(N,M)(xi*xi)
def P_int(M,N,x):
   integral=integrate.quad(func_in_P(M,N,x),0,10.0/x+10.0,limit=20) # we might use inf but we know that the difference is small
   return integral[0]

def xP(N,x):
   P0=sympy.sqrt(sympy.pi/2)*sympy.exp(x*x/2)*(1-sympy.erf(x/sympy.sqrt(2)))
   P1=x-x*x*P0
   yield P0
   P,Pprev,n= P1,P0,1
   while n<N:
      yield P
      P,Pprev,n = ((2*n+1)*P-sympy.diff(P,x,2) - n*Pprev )/(n+1) , P,n+1
      
x=sympy.Symbol("x")
Ns,M=range(10),0
P_anal_f=[P for P in xP(11,x)]
xxx=numpy.linspace(0.01,20.01,10) 
fig1=plt.figure()
ax1=fig1.gca()
fig2=plt.figure()
ax2=fig2.gca()
for N in Ns:
  P_ex  =numpy.array( [a*P_int(M,N,a) for a in xxx ] )
  P_anal=numpy.array( [a*P_anal_f[N]({x:a}).evalf() for a in xxx] )
  ax1.plot(xxx,numpy.abs(P_ex-P_anal) ,label="N=%i"%N)
  ax2.plot(xxx,P_ex,label="ex N=%i"%N)
  ax2.plot(xxx,P_anal,":",label="anal N=%i"%N)
ax1.legend()
plt.show()

[ГАЗ-67]

Не совсем понимаю - нужно посчитать любым способом или именно численным ?

[nestoklon]

ГАЗ-67 в сообщении #340472 писал(а):

Не совсем понимаю - нужно посчитать любым способом или именно численным ?

Нужно посчитать. Как угодно.
Под "посчитать" я понимаю написать функцию $P_N^M(x)$, которая давала бы численный ответ с разумной ($10^{-3}$ меня бы устроило) точностью. При достаточно больших (хотя бы до 30) $N$ и не очень маленьких (хотя бы до 10) $x$.
Аналитическое решение "в лоб" не подходит, так как оно численно неустойчиво (см. выше).

[ГАЗ-67]

$$\int_0^{\infty} e^{- px^2} x^N dx=\frac{(2n-1)!!}{2*(2p)^n}sqrt{\frac{\pi}{p}}$$

для чётных степеней

$$\int_0^{\infty} e^{- px^2} x^{2N+1} dx=\frac{N!}{2*p^{n+1}}}$$

для нечётных степеней
где-то так ...

[nestoklon]

ГАЗ-67 в сообщении #340905 писал(а):

где-то так ...

Простите, это шутка или вы не обратили внимание на такую мелочь как то, что параметр в экспоненте стоит при переменной интегрирования в первой степени? В результате чего интеграл не выражается в элементраных функциях. Что можно проверить, посчитав его при $N=0$ и получив функцию ошибок? Если вы не получили функцию ошибок, значит вы взяли интеграл неправильно.
Или там у Вас в знаменателе свёртка и я чего-то недопонимаю?.. Хотя свёртка с константой функцию не поменяет... Так что наверное что-то недопонимаю не я.

ГАЗ-67 в сообщении #340320 писал(а):

В показателе экспоненты выделите полный квадрат и сделайте соответствующую замену переменной .

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

[ewert]

А чем не устраивает просто численное интегрирование?... Тем более если Вам достаточно $10^{-3}$ -- пары десятков шагов и хватит. Где находится экстремум (для чистого $x^N$) -- нам известно, так что настройка вычислительной процедуры не должна вызвать затруднений...

[ГАЗ-67]

Прошу прощения , мозги от жары совсем засохли
На такую существенную деталь не обратил внимания ...

[nestoklon]

ewert в сообщении #340930 писал(а):

А чем не устраивает просто численное интегрирование?... Тем более если Вам достаточно $10^{-3}$ -- пары десятков шагов и хватит. Где находится экстремум (для чистого $x^N$) -- нам известно, так что настройка вычислительной процедуры не должна вызвать затруднений...

Спасибо, что заставили меня над этим задуматься. В итоге я понял, что через $\int_0^{\infty} e^{-\gamma x - \frac12 x^2} x^N dx$ просто нельзя считать. Следите за руками. Если мне нужна точность порядка $10^{-3}$ для интеграла с лагеррами, который порядка 1, это значит что мне надо считать все слагаемые с абсолютной точностью до $10^{-3}$. Это значит что интеграл с $x^N$ для $N$ равного, скажем, 10, мне тоже надо считать с точностью до $10^{-3}$. А максимум этой функции зависит от $N$ чуть медленнее чем $N^{\frac{N}2}$ (ну, может я ошибся и просто экспоненциально, не суть). То есть для довольно небольших $N$ относительная точность вычисления интеграла должна быть запредельной, чтобы гарантировать требуемую абсолютную точность.
Видимо, надо копаться в численно устойчивых алгоритмах вычисления лагерров и функции ошибок. Может, пойму как надо делать тут -- идеологически должно быть похоже.
А на первый взгляд такая тривиальная задача была...

(Оффтоп)

ГАЗ-67 в сообщении #340961 писал(а):

На такую существенную деталь не обратил внимания ...

Ничего, со всеми бывает. Я уж испугался что проворонил что-то совсем тривиальное.