2014 dxdy logo

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

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




Начать новую тему Ответить на тему
 
 Вычисление интегралов от полиномов Лагерра.
Сообщение19.07.2010, 00:00 
Заслуженный участник


19/07/08
1266
В ходе решения некоторой задачи возникла задача посчитать интегралы вида
$$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 


09/06/06
367
nestoklon в сообщении #339833 писал(а):
В принципе, есть вариант выразить эти интегралы через
$$R_N(\gamma)=\int_0^{\infty} e^{-\gamma x - \frac12 x^2} x^N dx,$$

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

 Профиль  
                  
 
 Re: Вычисление интегралов от полиномов Лагерра.
Сообщение22.07.2010, 09:29 
Заслуженный участник


19/07/08
1266
ГАЗ-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 


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

 Профиль  
                  
 
 Re: Вычисление интегралов от полиномов Лагерра.
Сообщение23.07.2010, 11:09 
Заслуженный участник


19/07/08
1266
ГАЗ-67 в сообщении #340472 писал(а):
Не совсем понимаю - нужно посчитать любым способом или именно численным ?

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

 Профиль  
                  
 
 Re: Вычисление интегралов от полиномов Лагерра.
Сообщение26.07.2010, 09:38 


09/06/06
367
$$\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 
Заслуженный участник


19/07/08
1266
ГАЗ-67 в сообщении #340905 писал(а):
где-то так ...

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

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

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

 Профиль  
                  
 
 Re: Вычисление интегралов от полиномов Лагерра.
Сообщение26.07.2010, 11:23 
Заслуженный участник


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

 Профиль  
                  
 
 Re: Вычисление интегралов от полиномов Лагерра.
Сообщение26.07.2010, 14:54 


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

 Профиль  
                  
 
 Re: Вычисление интегралов от полиномов Лагерра.
Сообщение26.07.2010, 19:23 
Заслуженный участник


19/07/08
1266
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 ] 

Модераторы: Toucan, maxal, PAV, Karan, Супермодераторы



Кто сейчас на конференции

Сейчас этот форум просматривают: нет зарегистрированных пользователей


Вы не можете начинать темы
Вы не можете отвечать на сообщения
Вы не можете редактировать свои сообщения
Вы не можете удалять свои сообщения
Вы не можете добавлять вложения

Найти:
Powered by phpBB © 2000, 2002, 2005, 2007 phpBB Group