Metoda całkowania numerycznego trudnej całki oscylacyjnej

Muszę ocenić liczbowo całkę poniżej:

gdzie ,x∈R+iλ,κ,ν>0. TutajKjest zmodyfikowaną funkcją Bessela drugiego rodzaju. W moim szczególnym przypadku mamλ=0,00313,κ=0,00825iν=0,33.$E(r) = r^4 (\lambda\sqrt{\kappa^2+r^2})^{-\nu-5/2} K_{-\nu-5/2}(\lambda\sqrt{\kappa^2+r^2})$$x \in \mathbb{R}_+$$\lambda, \kappa, \nu >0$$K$$\lambda = 0.00313$$\kappa = 0.00825$$\nu = 0.33$

Korzystam z MATLAB i wypróbowałem wbudowane funkcje integrali quadgk, co daje mi wiele błędów (patrz poniżej). Naturalnie próbowałem również wielu innych rzeczy, takich jak całkowanie przez części i sumowanie całek od do ( k + 1 ) x π .$kx\pi$$(k+1)x\pi$

Czy masz jakieś sugestie, którą metodę powinienem wypróbować?

AKTUALIZACJA (dodane pytania)
Przeczytałem artykuł z linkiem @Pedro i nie sądzę, aby było to zbyt trudne do zrozumienia. Mam jednak kilka pytań:

• Byłoby dobrze, aby skorzystać jako podstawa-elements * F k , w jednoczynnikowej Levin metody opisane?$x^k$$\psi_k$
• Czy zamiast tego mogę po prostu użyć metody Filon, ponieważ częstotliwość drgań jest stała?

Przykładowy kod
>> integral(@(r) sin(x*r).*sqrt(E(r)),0,Inf)
Warning: Reached the limit on the maximum number of intervals in use. Approximate
bound on error is 1.6e+07. The integral may not exist, or it may be difficult to
approximate numerically to the requested accuracy.
> In funfun\private\integralCalc>iterateScalarValued at 372
In funfun\private\integralCalc>vadapt at 133
In funfun\private\integralCalc at 84
In integral at 89

ans =

3.3197e+06

Napisałem własnego integratora, quadccktóry radzi sobie znacznie lepiej niż integratory Matlab z osobliwościami i zapewnia bardziej wiarygodne oszacowanie błędów.

Aby użyć go do rozwiązania problemu, wykonałem następujące czynności:

>> lambda = 0.00313; kappa = 0.00825; nu = 0.33;
>> x = 10;
>> E = @(r) r.^4.*(lambda*sqrt(kappa^2 + r.^2)).^(-nu-5/2) .* besselk(-nu-5/2,lambda*sqrt(kappa^2 + r.^2));
>> sincp = @(x) cos(x)./x - sin(x)./x.^2;
>> f = @(r) sincp(x*r) .* r .* sqrt( E(r) );

Ta funkcja fjest teraz twoją integracją. Zauważ, że właśnie przypisałem dowolną starą wartość x.

Aby zintegrować w nieskończonej domenie, stosuję podstawienie zmiennych:

>> g = @(x) f ( tan ( pi / 2 * x ) ) .* ( 1 + tan ( pi * x / 2 ).^2 ) * pi / 2;

gf$\infty$g

Następnie dzwonię do mojego integratora, quadccktóry może poradzić sobie z NaNs na obu końcach:

>> [ int , err , npoints ] = quadcc( g , 0 , 1 , 1e-6 )
int =
  -1.9552e+06
err =
   1.6933e+07
npoints =
       20761

Zauważ, że oszacowanie błędu jest ogromne, tzn. quadccNie ma większego zaufania do wyników. Patrząc na tę funkcję, nie jest to zaskakujące, ponieważ oscyluje ona przy wartościach o trzy rzędy wielkości powyżej rzeczywistej całki. Ponownie użycie innej transformacji przedziałowej może dać lepsze wyniki.

Możesz także przyjrzeć się bardziej szczegółowym metodom, takim jak to . Jest to nieco bardziej zaangażowane, ale zdecydowanie odpowiednia metoda dla tego rodzaju problemów.

Jak zauważa Pedro, metody typu Levina są najlepiej ustalonymi metodami dla tego rodzaju problemów.

Czy masz dostęp do Mathematica? W przypadku tego problemu Mathematica domyślnie je wykrywa i wykorzystuje:

In[1]:= e[r_] := 
 r^4 (l Sqrt[k^2 + r^2])^(-v - 5/2) BesselK[-v - 5/2, l Sqrt[k^2 + r^2]]

In[2]:= {l, k, v} = {0.00313, 0.00825, 0.33};

In[3]:= Block[{x = 10}, 
 NIntegrate[Sinc'[x r] r Sqrt[e[r]], {r, 0, \[Infinity]}, 
  PrecisionGoal -> 3]]

Out[3]= -112494.

Oto wykres w zakresie wartości x:

In[4]:= ListLinePlot[
 Table[NIntegrate[Sinc'[x r] r Sqrt[e[r]], {r, 0, \[Infinity]}, 
   PrecisionGoal -> 3], {x, .5, 10, 0.1}]]

Możesz również ręcznie określić konkretną metodę typu Levina, która ma zostać zastosowana, co w tym przypadku może przynieść nieznaczną poprawę wydajności:

In[5]:= method = {"LevinRule", "Kernel" -> {Cos[r x], Sin[r x]}, 
   "DifferentialMatrix" -> {{0, -x}, {x, 0}}, 
   "Amplitude" -> {(
     3497.878840962873` Sqrt[(
      r^4 BesselK[-2.17`, 
        0.00313` Sqrt[
         0.00006806250000000001` + r^2]])/(0.00006806250000000001` + 
        r^2)^1.415`])/
     x, -((3497.878840962873` Sqrt[(
       r^4 BesselK[-2.17`, 
         0.00313` Sqrt[
          0.00006806250000000001` + r^2]])/(0.00006806250000000001` + 
         r^2)^1.415`])/(r x^2))}, "AdditiveTerm" -> 0};

In[6]:= Block[{x = 10}, 
 NIntegrate[Sinc'[x r] r Sqrt[e[r]], {r, 0, \[Infinity]}, 
  PrecisionGoal -> 3, Method -> method]]

Out[6]= -112495.

Zobacz dokumentację, aby uzyskać szczegółowe informacje na temat metod typu Levin w Mathematica .

Jeśli nie masz dostępu do Mathematica, możesz napisać w Matlabie metodę Levina (lub inną specjalistyczną metodę oscylacyjną), jak sugeruje Pedro.

Czy korzystasz z biblioteki chebfun dla Matlaba? Właśnie dowiedziałem się, że zawiera implementację podstawowej metody Levin typu tutaj . Implementacja została napisana przez Olvera (jednego z ekspertów w dziedzinie kwadratur oscylacyjnych). Nie dotyczy osobliwości, adaptacyjnego podziału itp., Ale może być właśnie tym, czego potrzebujesz, aby zacząć.

Transformacja zalecana przez Pedro to świetny pomysł. Czy próbowałeś pobawić się parametrami w funkcji „quadgk” Matlaba? Na przykład, używając transformacji Pedro, możesz wykonać następujące czynności:
quadgk(f, 0.0+eps, 1.0-eps, 'AbsTol', eps, 'MaxIntervalCount', 100000)
Użycie tego daje mi rozwiązanie:
-2184689.50220729
i zajmuje tylko 0,8 sekundy (używając wartości wymienionych powyżej: x = 10)
Walter Gander i Walter Gautschi mają artykuł na temat kwadratury adaptacyjnej z Matlabem kod, którego możesz również użyć (link tutaj )