Metoda całkowania numerycznego trudnej całki oscylacyjnej

25

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

0sinc(xr)rE(r)dr

gdzie ,xR+iλ,κ,ν>0. TutajKjest zmodyfikowaną funkcją Bessela drugiego rodzaju. W moim szczególnym przypadku mamλ=0,00313,κ=0,00825iν=0,33.E(r)=r4(λκ2+r2)ν5/2Kν5/2(λκ2+r2)xR+λ,κ,ν>0Kλ=0.00313κ=0.00825ν=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π(k+1)xπ

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?xkψ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

torbonde
źródło
Co to jest w twojej całce? x
Pedro
Każda dodatnia, rzeczywista liczba. Właśnie zaktualizowałem swój post.
torbonde
Jeśli potrafisz pokazać kod i błędy, prawdopodobnie nie jest zbyt trudno rozwiązać większość z nich. Oczywiście, proszę najpierw przeczytaj uważnie błąd i sprawdź, czy możesz sam go usunąć.
Dennis Jaheruddin
Później dodam komentarz z kodem i błędami. Albo jutro.
torbonde
Okej, więc zapomniałem. Ale teraz zaktualizowałem swój post o przykład (podzieliłem całkę na dwie części, obliczając jawnie ). sinc
torbonde

Odpowiedzi:

12

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;

gfg

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.

Pedro
źródło
Dziękuję Ci bardzo. Przyjrzę się różnym metodom. Dla moich celów błąd nie musi być tak mały, jak jest to standard w eq integral(chyba 1e-10), ale 1,7e + 07 jest nadal naprawdę bardzo duży. Być może kolejna transformacja przyniesie dobre efekty, jak wspominasz.
torbonde
@ cimrg.joe: Należy pamiętać, że oszacowanie błędu jest oszacowaniem błędu bezwzględnego opartego między innymi na maksymalnych wartościach bezwzględnych całki. W niektórych ekstremalnych przypadkach zwracana wartość może być całkiem dobra. Jeśli szukasz dziesięciu cyfr dokładności, zdecydowanie zalecam stosowanie metod typu Levin, o których wspomniałem na końcu mojego postu.
Pedro
Może nie potrzebuję dziesięciu cyfr dokładności, ale myślę, że potrzebuję co najmniej pięciu. Czy twoja metoda może to wytworzyć?
torbonde
Ta metoda nie gwarantuje tego rodzaju precyzji dla całki, ponieważ wartości na prawym końcu przedziału są o kilka rzędów wielkości większe niż sama całka.
Pedro
11

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}]]

Wykres od x = 0,5 do x = 10

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 .

Andrew Moylan
źródło
Niestety nie mam dostępu do Mathematica - tylko MATLAB. Po prostu zaktualizuję moje pytanie, dodając kilka dodatkowych pytań dotyczących artykułu @Pedro, z którym jest powiązane.
torbonde
OK, jak mówisz, będziesz musiał zadowolić się Matlabem. Dodam kolejną odpowiedź na ten temat.
Andrew Moylan,
5

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ąć.

Andrew Moylan
źródło
Myślałem o wdrożeniu metody Levina, ale nie jestem pewien, czy jestem gotowy na to wyzwanie. Myślę, że muszę trochę lepiej zrozumieć tę metodę. Może mógłbym porozmawiać o tym z moim doradcą. W każdym razie powodem, dla którego pytałem o metody Filon, jest to, że wydają się one łatwiejsze do wdrożenia. A ponieważ nie potrzebuję bardzo wysokiej dokładności, ale jest to część mojej pracy magisterskiej, trudność waży.
torbonde
Przejrzałem bibliotekę chebfun (imponującą) i przykład integracji z Levinem. Ale nie mogę tego uruchomić. I rzeczywiście pisał pytanie dotyczące go tutaj .
torbonde
0

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 )

Joel Vroom
źródło