2014 dxdy logo

Научный форум dxdy

Математика, Физика, Computer Science, Machine Learning, LaTeX, Механика и Техника, Химия,
Биология и Медицина, Экономика и Финансовая Математика, Гуманитарные науки




 
 Вычисление интегралов от полиномов Лагерра.
Сообщение19.07.2010, 00:00 
В ходе решения некоторой задачи возникла задача посчитать интегралы вида
$$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 не получится, так как в какой-то момент придётся умножать большое значение экспоненты на маленькое значение функции параболического цилиндра, что приводит к ещё большим проблемам, несмотря на "аналитическое" решение.
Вопрос именно в том, как это можно посчитать.

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

 
 
 
 Re: Вычисление интегралов от полиномов Лагерра.
Сообщение22.07.2010, 07:32 
nestoklon в сообщении #339833 писал(а):
В принципе, есть вариант выразить эти интегралы через
$$R_N(\gamma)=\int_0^{\infty} e^{-\gamma x - \frac12 x^2} x^N dx,$$

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

 
 
 
 Re: Вычисление интегралов от полиномов Лагерра.
Сообщение22.07.2010, 09:29 
ГАЗ-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()

 
 
 
 Re: Вычисление интегралов от полиномов Лагерра.
Сообщение23.07.2010, 07:19 
Не совсем понимаю - нужно посчитать любым способом или именно численным ?

 
 
 
 Re: Вычисление интегралов от полиномов Лагерра.
Сообщение23.07.2010, 11:09 
ГАЗ-67 в сообщении #340472 писал(а):
Не совсем понимаю - нужно посчитать любым способом или именно численным ?

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

 
 
 
 Re: Вычисление интегралов от полиномов Лагерра.
Сообщение26.07.2010, 09:38 
$$\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}}}$$ для нечётных степеней
где-то так ...

 
 
 
 Re: Вычисление интегралов от полиномов Лагерра.
Сообщение26.07.2010, 10:32 
ГАЗ-67 в сообщении #340905 писал(а):
где-то так ...

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

ГАЗ-67 в сообщении #340320 писал(а):
В показателе экспоненты выделите полный квадрат и сделайте соответствующую замену переменной .

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

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

 
 
 
 Re: Вычисление интегралов от полиномов Лагерра.
Сообщение26.07.2010, 14:54 
Прошу прощения , мозги от жары совсем засохли :oops:
На такую существенную деталь не обратил внимания ... :oops:

 
 
 
 Re: Вычисление интегралов от полиномов Лагерра.
Сообщение26.07.2010, 19:23 
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 писал(а):
На такую существенную деталь не обратил внимания ... :oops:
Ничего, со всеми бывает. Я уж испугался что проворонил что-то совсем тривиальное.

 
 
 [ Сообщений: 10 ] 


Powered by phpBB © 2000, 2002, 2005, 2007 phpBB Group