Интеграл Mathematica со многими особенностями

Как лучше всего заставить Mathematica 7 или 8 выполнять интеграл?

NIntegrate[Exp[-x]/Sin[Pi x], {x, 0, 50}]

У каждого целого числа есть полюса, и нам нужно значение принципа Коши. Идея состоит в том, чтобы получить хорошее приближение для интеграла от 0 до бесконечности.

С Integrate есть вариант PrincipleValue -> True.

С помощью NIntegrate я могу дать ему вариант Exclusions -> (Sin[Pi x] == 0) или вручную указать полюса с помощью

NIntegrate[Exp[-x]/Sin[Pi x], Evaluate[{x, 0, Sequence@@Range[50], 50}]]

Исходная команда и два приведенных выше приема NIntegrate дают результат 60980 +/- 10. Но все выплевывают ошибки. Каков наилучший способ получить быстрый надежный результат для этого интеграла без того, чтобы Mathematica давала ошибки?


person Simon    schedule 07.04.2011    source источник
comment
К вашему сведению, если я попробую разные значения WorkPrecision, я получу совершенно разные результаты.   -  person Mr.Wizard    schedule 07.04.2011


Ответы (3)


Саймон, есть ли основания полагать, что ваш интеграл сходится?

In[52]:= f[k_Integer, eps_Real] := 
 NIntegrate[Exp[-x]/Sin[Pi x], {x, k + eps, k + 1 - eps}]

In[53]:= Sum[f[k, 1.0*10^-4], {k, 0, 50}]

Out[53]= 2.72613

In[54]:= Sum[f[k, 1.0*10^-5], {k, 0, 50}]

Out[54]= 3.45906

In[55]:= Sum[f[k, 1.0*10^-6], {k, 0, 50}]

Out[55]= 4.19199

Похоже, проблема в x==0. Разделение подынтегрального выражения k+eps на k+1-eps для целых значений k:

In[65]:= int = 
 Sum[(-1)^k Exp[-k ], {k, 0, Infinity}] Integrate[
   Exp[-x]/Sin[Pi x], {x, eps, 1 - eps}, Assumptions -> 0 < eps < 1/2]

Out[65]= (1/((1 + 
   E) (I + \[Pi])))E (2 E^(-1 + eps - I eps \[Pi])
     Hypergeometric2F1[1, (I + \[Pi])/(2 \[Pi]), 3/2 + I/(2 \[Pi]), 
     E^(-2 I eps \[Pi])] + 
   2 E^(I eps (I + \[Pi]))
     Hypergeometric2F1[1, (I + \[Pi])/(2 \[Pi]), 3/2 + I/(2 \[Pi]), 
     E^(2 I eps \[Pi])])

In[73]:= N[int /. eps -> 10^-6, 20]

Out[73]= 4.1919897038160855098 + 0.*10^-20 I

In[74]:= N[int /. eps -> 10^-4, 20]

Out[74]= 2.7261330651934049862 + 0.*10^-20 I

In[75]:= N[int /. eps -> 10^-5, 20]

Out[75]= 3.4590554287709991277 + 0.*10^-20 I

Как видите, имеется логарифмическая особенность.

In[79]:= ser = 
 Assuming[0 < eps < 1/32, FullSimplify[Series[int, {eps, 0, 1}]]]

Out[79]= SeriesData[eps, 0, {(I*(-1 + E)*Pi - 
     2*(1 + E)*HarmonicNumber[-(-I + Pi)/(2*Pi)] + 
          Log[1/(4*eps^2*Pi^2)] - 2*E*Log[2*eps*Pi])/(2*(1 + E)*Pi), 
     (-1 + E)/((1 + E)*Pi)}, 0, 2, 1]

In[80]:= Normal[
  ser] /. {{eps -> 1.*^-6}, {eps -> 0.00001}, {eps -> 0.0001}}

Out[80]= {4.191989703816426 - 7.603403526913691*^-17*I, 
 3.459055428805136 - 
     7.603403526913691*^-17*I, 
 2.726133068607085 - 7.603403526913691*^-17*I}

EDIT Out[79] приведенного выше кода дает расширение ряда для eps->0, и если объединить эти два логарифмических члена, мы получим

In[7]:= ser = SeriesData[eps, 0, 
       {(I*(-1 + E)*Pi - 2*(1 + E)*HarmonicNumber[-(-I + Pi)/(2*Pi)] + 
              Log[1/(4*eps^2*Pi^2)] - 2*E*Log[2*eps*Pi])/(2*(1 + E)*
       Pi), 
         (-1 + E)/((1 + E)*Pi)}, 0, 2, 1]; 

In[8]:= Collect[Normal[PowerExpand //@ (ser + O[eps])], 
 Log[eps], FullSimplify]

Out[8]= -(Log[eps]/\[Pi]) + (
 I (-1 + E) \[Pi] - 
  2 (1 + E) (HarmonicNumber[-((-I + \[Pi])/(2 \[Pi]))] + 
     Log[2 \[Pi]]))/(2 (1 + E) \[Pi])

Очевидно, что -Log[eps]/Pi исходит из полюса x==0. Итак, если вычесть это, как метод основного значения делает это для других полюсов, вы получите конечное значение:

In[9]:= % /. Log[eps] -> 0

Out[9]= (I (-1 + E) \[Pi] - 
 2 (1 + E) (HarmonicNumber[-((-I + \[Pi])/(2 \[Pi]))] + 
    Log[2 \[Pi]]))/(2 (1 + E) \[Pi])

In[10]:= N[%, 20]

Out[10]= -0.20562403655659928968 + 0.*10^-21 I

Конечно, этот результат трудно проверить численно, но вы можете знать больше, чем я, о вашей проблеме.

ИЗМЕНИТЬ 2

Это редактирование предназначено для подтверждения входных данных In[65], которые вычисляют исходный регуляризованный интеграл. Мы вычисляем

Sum[ Integrate[ Exp[-x]/Sin[Pi*x], {x, k+eps, k+1-eps}], {k, 0, Infinity}] ==  
  Sum[ Integrate[ Exp[-x-k]/Sin[Pi*(k+x)], {x, eps, 1-eps}], {k, 0, Infinity}] ==
  Sum[ (-1)^k*Exp[-k]*Integrate[ Exp[-x]/Sin[Pi*x], {x, eps, 1-eps}], 
       {k, 0, Infinity}] == 
  Sum[ (-1)^k*Exp[-k], {k, 0, Infinity}] * 
     Integrate[ Exp[-x]/Sin[Pi*x], {x, eps, 1-eps}]

В третьей строке использовалось значение Sin[Pi*(k+x)] == (-1)^k*Sin[Pi*x] для целого числа k.

person Sasha    schedule 07.04.2011
comment
Привет Саша, спасибо за всю вашу работу над этой проблемой. На самом деле я получил эту проблему от математика, который сказал, что использует Mma4, потому что у него есть пакет, который обеспечивает лучшую поддержку, чем более поздняя встроенная поддержка числовых интегралов PrincipleValue. Он привел мне этот пример, и теперь, когда вы заставили меня задуматься об этом, я согласен... это не очень хороший пример! Проблема при lim x->+0 заключается в том, что нет соответствующего lim x->-0, который отменял бы полюс. - person Simon; 09.04.2011
comment
У меня есть вопрос о вашем In[65], где вы определяете int. Sum кажется не связанным со всем остальным... Я думаю, что это приводит к проблемам в остальной части вашего анализа - но не к вашим выводам! - person Simon; 09.04.2011
comment
@Simon In[65] легко оправдать, разбив интеграл 0..inf на 0+..1- и 1+..2- и так далее. Затем оно становится суммой по k интегралов от Exp[-k-x]/Sin[Pi*(k+x)] для x от eps до 1-eps. Но Sin[Pi*(k+x)] == (-1)^kSin[Pix], и подынтегральная функция становится (-1)^kExp[-k]* Exp[-x]/Sin[Pix]. Таким образом, сумма и интеграл теперь могут быть вычислены независимо, как в [65]. Я считаю, что ваш вопрос стоит внести в редактирование, объясняющее этот шаг. - person Sasha; 09.04.2011

Саймон, я не уделял много времени твоему интегралу, но тебе следует попробовать взглянуть на аппроксимацию стационарной фазы. У вас есть плавная функция (эксп) и сильно колебательная функция (синус). Теперь работа заключается в том, чтобы превратить 1/sin(x) в форму exp(if(x)).

В качестве альтернативы вы можете использовать расширение серии cosecant (не действует на полюсах):

In[1]:=Series[Csc[x], {x, 0, 5}]
(formatted) Out[1]=1/x + x/6 + 7/360 x^3 + 31/15120 x^5 +O[x]^6

Обратите внимание, что для всех m>-1 у вас есть следующее:

In[2]:=Integrate[x^m Exp[-x], {x, 0, Infinity}, Assumptions -> m > -1]
Out[2]=Gamma[1+m]

Однако, суммируя ряд с коэффициентами косеканса (из википедии), не включая 1/x Exp[-x] случай, который не сходится на [0,Infinity].

c[m_] := (-1)^(m + 1) 2 (2^(2 m - 1) - 1) BernoulliB[2 m]/Factorial[2 m];
Sum[c[m] Gamma[1 + 2 m - 1], {m, 1, Infinity}]

тоже не сходится...

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

person abcd    schedule 07.04.2011
comment
Спасибо Р.М. Интеграл был плохим примером, данным мне кем-то, чьим математическим способностям я доверял... Есть полюс в нуле, где мы берем только односторонний предел. Это означает, что мы не можем получить PrincipleValue для этого полюса... - person Simon; 09.04.2011

Я должен согласиться с Сашей, интеграл не выглядит сходящийся. Однако если исключить x == 0 и разбить интеграл на части

Integrate[Exp[-x]/Sin[Pi x], {x, n + 1/2, n + 3/2}, PrincipalValue -> True]

где n >= 0 && Element[n, Integers], то вроде можно получить чередующийся ряд

I Sum[ (-1/E)^n, {n, 1, Infinity}] == - I / (1 + E )

Сейчас я отнес его только к n == 4, но выглядит разумно. Однако для приведенного выше интеграла с Assumptions -> Element[n, Integers] && n >= 0 Mathematica дает

If[ 2 n >= 1, - I / E, Integrate[ ... ] ]

что просто не соответствует отдельным случаям. В качестве дополнительного примечания, если полюс лежит на границе области интегрирования, то есть ваши пределы равны {x, n, n + 1}, вы получите только DirectedInfinitys. Быстрый взгляд на график подразумевает, что вы с пределами {x, n, n + 1} у вас есть только строго положительное или отрицательное подынтегральное выражение, поэтому бесконечное значение может быть связано с отсутствием компенсации, которую дает вам {x, n + 1/2, n + 3/2}. Однако проверка с помощью {x, n, n + 2} выдает только невычисленный интеграл.

person rcollyer    schedule 07.04.2011
comment
Я не получаю чередующийся ряд... скорее один со всеми положительными условиями (см. мой ответ выше). Интеграл расширяется до суммы Gamma функций (которая не сходится) плюс интеграл от Exp[-x]/x члена, который также не сходится из-за сингулярности в начале координат. - person abcd; 07.04.2011
comment
@Р. М., на самом деле мы используем разные регионы интеграции. Я исключаю часть x == 0, начиная с x == 1/2 и интегрируя интервалы шириной в 1 единицу; сумма которых выглядит как знакопеременный ряд, который сходится. - person rcollyer; 07.04.2011
comment
Спасибо rcollyer. Вы правы, точка x=0 — это односторонний предел, который не может быть уравновешен другой стороной. И это Integrate[Exp[-x]/Sin[Pi x], {x, n + 1/2, n + 3/2}, PrincipalValue -> True, Assumptions -> n \[Element] Integers] дает сложную константу, умноженную на (-1/E)^n. Но это не соответствует NIntegrate[Exp[-x]/Sin[Pi x], {x, n + 1/2, n, n + 3/2}, WorkingPrecision -> 100, MaxRecursion -> 40], что дает ответ, сильно зависящий от WorkingPrecision. В общем, это был плохой вопрос/пример. - person Simon; 09.04.2011