Numeryczna ocena całki silnie oscylacyjnej

11

W tym zaawansowanym kursie na temat zastosowań teorii funkcji złożonych w jednym punkcie ćwiczenia całka silnie oscylacyjna

I(λ)=-sałata(λsałatax)grzechxxrex

należy aproksymować dla dużych wartości λ stosując metodę punktu siodłowego w płaszczyźnie złożonej.

Ze względu na bardzo oscylacyjny charakter całka ta jest bardzo trudna do oceny przy użyciu większości innych metod. Są to dwa fragmenty wykresu całki dla λ=10 w różnych skalach:

cos (10 cos (x)) sin (x) / x

Wiodącym asymptotycznym przybliżeniem rzędu jest

ja1(λ)=sałata(λ-14π)2)πλ

a kolejne (znacznie mniejsze) udoskonalenie dodaje termin

ja2)(λ)=18grzech(λ-14π)2)πλ3)

Wykres przybliżonych wartości w funkcji wygląda następująco:λ

I (lambda) ok

Teraz pojawia się moje pytanie: aby wizualnie zobaczyć, jak dobre jest przybliżenie, chciałbym porównać je z „rzeczywistą wartością” całki, a dokładniej z dobrym przybliżeniem do tej samej całki przy użyciu niezależnego algorytmu. Ze względu na niewielką korektę subleadingową spodziewałbym się, że będzie to naprawdę bliskie.

Próbowałem oszacować całkę dla niektórych λ przy użyciu innych algorytmów, ale z bardzo małym powodzeniem: Mathematica i Matlab przy użyciu domyślnego integratora numerycznego nie są w stanie wygenerować znaczącej wartości (i podać to jawnie), mpmath przy użyciu zarówno podwójnie wykładniczej tanh(sinh) podstawienie i metoda Gaussa-Legendre'a daje bardzo głośne wyniki, chociaż ma niewielką tendencję do oscylowania wokół wartości, które daje metoda punktu siodłowego, ponieważ ten wykres może pokazywać:

mpmath ok

W końcu spróbowałem szczęścia z integratorem Monte-Carlo, wykorzystując próbkę ważności, którą zaimplementowałem, ale nie udało mi się też uzyskać stabilnych wyników.

Czy ktoś ma pojęcie o tym, jak tę całkę można niezależnie oszacować dla dowolnej stałej wartości λ>1 lub więcej?

doetoe
źródło
Czy funkcja jest parzysta?
nicoguaro
Tak, to jest nawet
doetoe
Czy próbowałeś zamienić całkę w ODE?
nicoguaro
1
Nie, różnicowanie wr a następnie rozwiązywanie równania różniczkowego numerycznie. x
nicoguaro
1
Twój pierwszy wykres wydaje się pełnić inną funkcję niż twój integrand. Mianowicie wydaje się, że zastąpione przez λ x . Tj. Wykres ma funkcję x ( cos ( λ x cos x ) sinc x )λλxx(sałata(λxsałatax)sincx)
Ruslan

Odpowiedzi:

12

Użyj twierdzenia Plancherela, aby ocenić tę całkę.

Podstawową ideą jest to, że dla dwóch funkcji fa,sol ,

ja=-fa(x)sol(x)rex=-fa(k)sol(k)rek

gdzie fa,sol są transformatami Fouriera fa,sol . Obie funkcje mają stosunkowo niewielkie wsparcie w dziedzinie spektralnej. Tutaj grzechx/xrect(k) i sałata(λsałatax) powinny mieć analityczną transformację Fouriera (lub szereg), podobnie jak rozszerzenie Jacobi-Anger . Możesz obciąć nieskończoną serię o wartości około λ ze względu na nadwykładniczy rozkład funkcji Bessela |jotn(x)|dla n>|x|. Mam nadzieję że to pomoże.

Edycja : Właściwie powinieneś użyć tutaj reprezentacji serii Fouriera zamiast transformacji. Ścieżka transformacji prowadzi do uzyskania asymptotycznej reprezentacji, którą już masz (okazuje się, że jest to tylko πjot0(λ) ). Powyższe twierdzenie Plancherela działa również dla szeregów Fouriera z domeną integracji [0,2)π] na ostatniej całce.

smh
źródło
Dzięki, to bardzo dobry pomysł!
doetoe
7

Kluczem do oceny całek oscylacyjnych jest obcięcie całki w odpowiednim punkcie. W tym przykładzie musisz wybrać górną granicę formy

πN.+π2)
Zanim wyjaśnię, dlaczego to powinno działać, pozwól mi najpierw pokazać, że faktycznie przynosi dobre wyniki.

Asymptotyki

Łatwo zgadnąć, że szereg asymptotyczny ma postać

ja(λ)2)πλ[sałata(λ-π4)+do1grzech(λ-π4)λ+do2)sałata(λ-π4)λ2)+do3)grzech(λ-π4)λ3)+]
Aby sprawdzić liczbowo, żedo1=18 wystarczy wykreślić różnicę między całkową i wiodącą ekspresją asymptotyczną.

int := NIntegrate[Cos[l*Cos[x]]*Sinc[x], {x, 0, 20.5*Pi}]; 
Plot[{l*(Sqrt[2*l/Pi]*int - Cos[l-Pi/4]), Sin[l-Pi/4]/8}, {l, Pi/4, 20}]

Jako wynik otrzymujesz całkiem niezły sinus, który pokrywa się z tym, który wyprowadziłeś powyżej.

18

Jeśli chcesz znaleźć następujące współczynniki, w razie potrzeby nieco bardziej zaawansowany fragment kodu. Poniższy kod ma na celu przyjęcie kilku wysokich górnych wartości granicznych i „uśrednienie” ich wyników.

J[l_?NumericQ] := Block[{n=500},
  f[k_] := NIntegrate[Cos[l*Cos[x]]*Sinc[x], {x,0,(n+k)*Pi+Pi/2},
  Method->{"DoubleExponential"}, AccuracyGoal->14, MaxRecursion->100];
  1/2*((f[0]+f[1])/2+(f[1]+f[2])/2)
]
t = Table[{l, l^2*(Sqrt[2*l/Pi]*J[l] - Cos[l-Pi/4] - 1/8*Sin[l-Pi/4]/l)}, 
    {l, 4*Pi+Pi/4, 12*Pi+Pi/4, Pi/36}];
Fit[t, Table[Cos[l-Pi/4+Pi/2*n]/l^n, {n, 0, 10}], l]

c2=9128,c3=751024,c4=367532768,

Wyjaśnienie

Prosty przykład

S(x)=0xsin(y)ydy.
S()=π2

sinus

S(x)

SN=n=1N(1)nn.
SSN+12(1)N+1N+1.
S(x)0πN+π2sinxxdx
max|S(x)|

Twój problem

jax0(λ)=2)0x0sałata(λsałata(x))sinc(x)rex
x0=πN.+π2)λ=12π

tab = Table[{x0, 2*NIntegrate[Cos[12*Pi*Cos[x]]*Sinc[x], {x, 0, x0}, 
     Method->{"DoubleExponential"}, AccuracyGoal->12, MaxRecursion->100]},
    {x0, 10*Pi+Pi/2, 30*Pi+Pi/2, Pi}];
tab1 = Table[(tab[[i]] + tab[[i+1]])/2, {i,1,Length[tab]-1}];
ListPlot[{tab, tab1}]

zgodnie z

S.N.=12)(S.N.+S.N.+1)
S.N.

David Saykin
źródło
Ładny! Czy instruktorzy kursu są prawdziwymi profesorami? Ich trasa jest fantastyczna, choć bardzo trudna i szybka
doetoe
@doetoe tak, jestem studentem Konstantin. Udostępnił mi link do twojego pytania.
David Saykin
6

Działa tutaj metoda Ooury dla całek sinusoidalnych Fouriera, patrz:

Ooura, Takuya i Masatake Mori, Solidna formuła podwójnie wykładnicza dla całek typu Fouriera. Journal of obliczeniowej i stosowanej matematyki 112.1-2 (1999): 229-241.

Napisałem implementację tego algorytmu, ale nigdy nie wkładałem pracy, aby uzyskać go szybko (powiedzmy, buforowanie węzłów / wag), ale mimo to otrzymuję spójne wyniki we wszystkim poza precyzją zmiennoprzecinkową:

float     = 0.0154244
double    = 0.943661538060268
long dbl  = 0.943661538066058702
quad      = 0.943661538066060288748574485677942
oct       = 0.943661538066060288748574485680878906503533004997613278231689169604876
asymptotic= 0.944029734

Oto kod:

#include <iostream>
#include <boost/math/quadrature/ooura_fourier_integrals.hpp>
#include <boost/math/constants/constants.hpp>
#include <boost/multiprecision/float128.hpp>
#include <boost/multiprecision/cpp_bin_float.hpp>

template<class Real>
Real asymptotic(Real lambda) {
    using std::sin;
    using std::cos;
    using boost::math::constants::pi;
    Real I1 = cos(lambda - pi<Real>()/4)*sqrt(2*pi<Real>()/lambda);
    Real I2 = sin(lambda - pi<Real>()/4)*sqrt(2*pi<Real>()/(lambda*lambda*lambda))/8;
    return I1 + I2;
}

template<class Real>
Real osc(Real lambda) {
    using std::cos;
    using boost::math::quadrature::ooura_fourier_sin;
    auto f = [&lambda](Real x)->Real { return cos(lambda*cos(x))/x; };
    Real omega = 1;
    Real Is = ooura_fourier_sin<decltype(f), Real>(f, omega);
    return 2*Is;
}

template<class Real>
void print(Real val) {
   std::cout << std::defaultfloat;
   std::cout << std::setprecision(std::numeric_limits<Real>::digits10);
   std::cout <<  val <<  " = " << std::hexfloat << val;
}

int main() {
    using boost::multiprecision::float128;
    float128  s = 7;
    std::cout << "Asymptotic = " << std::setprecision(std::numeric_limits<float128>::digits10) << asymptotic(s) << "\n";
    std::cout << "float precision = ";
    print(osc(7.0f));
    std::cout << "\n";
    std::cout << "double precision= ";
    print(osc(7.0));
    std::cout << "\n";
    std::cout << "long double     = ";
    print(osc(7.0L));
    std::cout << "\n";
    print(osc(s));

    print(osc(boost::multiprecision::cpp_bin_float_oct(7)));
    std::cout << "\n";
}

λ0wprowadź opis zdjęcia tutaj

użytkownik14717
źródło
Dzięki, to jest naprawdę miłe! Jeszcze go nie uruchomiłem, moja instalacja doładowania nie jest kompatybilna, ale teraz pobieram najnowszą wersję.
doetoe
Dla pewności: w 23 masz cos (lambda * cos (x)) / x, bez współczynnika sin (x) z całki. Czy to ooura_fourier_sin, który przyjmuje ten czynnik sin (x), aby pomnożyć przekazaną mu całkę?
doetoe
Mam to działa. Wydaje się, że to i jego zależności są tylko nagłówkami, więc nawet nie musiałem instalować ani kompilować (z wyjątkiem plików wykonywalnych). Mam nadzieję, że zostanie włączony do ulepszenia!
doetoe
grzech(x)
@doetoe: Scalono w Boost 1.71. Interfejs API różni się nieco od tej odpowiedzi.
user14717,