Avaliação numérica da integral altamente oscilatória

Em deste curso avançado sobre aplicações da teoria de função complexa em um ponto em um exercício integral altamente oscilatório

$$I(\lambda)=\int_{-\infty}^{\infty} \cos (\lambda \cos x) \frac{\sin x}{x} d x$$

deve ser aproximado para grandes valores de usando o método do ponto de sela no plano complexo.$\lambda$

Devido à sua natureza altamente oscilatória, é muito difícil avaliar essa integral usando a maioria dos outros métodos. Estes são dois fragmentos do gráfico do integrando para em escalas diferentes:$\lambda = 10$

Uma aproximação assintótica de ordem principal é

$$I_{1}(\lambda) = \cos \left(\lambda-\frac{1}{4} \pi\right) \sqrt{\frac{2 \pi}{\lambda}}$$

e um refinamento adicional (muito menor) adiciona o termo

$$I_2(\lambda)=\frac{1}{8} \sin \left(\lambda-\frac{1}{4} \pi\right) \sqrt{\frac{2 \pi}{\lambda^{3}}}$$

Um gráfico dos valores aproximados em função de tem a seguinte aparência:$\lambda$

Agora vem a minha pergunta: para ver visualmente quão boa é a aproximação, eu gostaria de compará-la com o "valor real" da integral ou, mais precisamente, com uma boa aproximação à mesma integral usando um algoritmo independente. Devido à pequenez da correção subliminar, eu esperaria que isso fosse bem próximo.

Tentei avaliar a integral para alguns usando outros algoritmos, mas com muito pouco sucesso: o Mathematica e o Matlab, usando o integrador numérico padrão, não conseguem produzir um valor significativo (e relatam isso explicitamente), mpmath usando os dois pontos de exponencial substituição e o método Gauss-Legendre produzem resultados muito ruidosos, embora tenha uma leve tendência a oscilar em torno dos valores que o método do ponto de sela fornece, pois esse gráfico pode mostrar:$\lambda$$\tanh(\sinh)$

Finalmente, tentei a sorte com um integrador de Monte-Carlo usando um exemplo de importância que implementei, mas também não consegui obter resultados estáveis.

Alguém tem uma idéia de como essa integral pode ser avaliada independentemente para qualquer valor fixo de ou menos?$\lambda > 1$

Use o teorema de Plancherel para avaliar esta integral.

A idéia básica é que, para duas funções ,$f,g$

$$I=\int_{-\infty}^{\infty} f(x) g^*(x)dx = \int_{-\infty}^{\infty} F(k) G^*(k) dk$$

onde são as transformadas de Fourier de . Suas funções têm suporte relativamente pequeno no domínio espectral. Aqui, e devem ter uma transformação analítica de Fourier (ou série), como a expansão Jacobi-Anger . Você pode truncar a série infinita em termos aproximadamente devido ao decaimento função Besselpara. Espero que isto ajude.$F,G$$f,g$$\sin x / x \rightarrow \text{rect}(k)$$\cos(\lambda\cos x)$$\lambda$$|J_n(x)|$$n>|x|$

Editar : Na verdade, você deve usar as representações da série Fourier aqui em vez das transformações. O caminho da transformação leva a derivar a representação assintótica que você já possui (acontece que é apenas ). O teorema de Plancherel acima também funciona para as séries de Fourier com um domínio de integração de na última integral.$\pi J_0(\lambda)$$[0,2\pi]$

A chave para a avaliação das integrais oscilatórias é truncar a integral no ponto certo. Neste exemplo, você precisa escolher o limite superior do formulário Antes de explicar por que deve funcionar, deixe-me mostrar primeiro que ele realmente produz bons resultados.

$$\pi\mathbb{N}+\frac{\pi}{2}$$

Assintóticos

É fácil adivinhar que séries assintóticas têm a forma Para verificar numericamente se é suficiente traçar a diferença entre uma expressão assintótica integral e principal.

$$I(\lambda) \sim \sqrt{\frac{2\pi}{\lambda}}\left[\cos\left(\lambda-\frac{\pi}{4}\right) + c_1 \frac{\sin\left(\lambda-\frac{\pi}{4}\right)}{\lambda} + c_2 \frac{\cos\left(\lambda-\frac{\pi}{4}\right)}{\lambda^2} + c_3 \frac{\sin\left(\lambda-\frac{\pi}{4}\right)}{\lambda^3} + \dots \right]$$ $c_1 = \frac{1}{8}$

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

Como resultado, você obtém um seno bastante agradável, que coincide com o que você derivou acima.

Se você deseja encontrar os seguintes coeficientes, um pedaço de código um pouco mais sofisticado, se necessário. A idéia do código abaixo é usar vários valores-limite superiores altos e "mediar" seus resultados.

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]

Isso produz a seguinte resposta.

$$\boxed{ \boxed{ c_2 = -\frac{9}{128}, \qquad c_3 = -\frac{75}{1024}, \qquad c_4 = \frac{3675}{32768}, \quad\dots }}$$

Explicação

Exemplo simples

Para a ilustração, vou usar um exemplo mais simples da integral senoidal Deixe-me imaginar que estou interessado no valor , mas não o conheço.

$$S(x) = \int_{0}^x \frac{\sin(y)}{y} dy.$$ $S(\infty)=\frac{\pi}{2}$

Você vê que oscila em torno de seu valor limite de forma semelhante à soma parcial de alternância na série de sinais oscila com o corte superior. Quando você deseja estimar essa soma, de acordo com o método de aceleração da série Euler, você deve tomar Ou em termos de função seno-integral, você deve integrá-lo até o ponto entre o máximo e o mínimo de oscilações. Como é claramente visto no gráfico, esse ponto é dado aproximadamente por para grandes valores de argumento. De maneira mais geral, esse ponto é aquele em queocorre.$S(x)$

$$S_N = \sum_{n=1}^N \frac{(-1)^n}{n}.$$ $$S \approx S_N + \frac{1}{2}\frac{(-1)^{N+1}}{N+1}.$$ $$S(x) \approx \int_0^{\pi\mathbb{N}+\frac{\pi}{2}} \frac{\sin x}{x}dx$$ $\max |S'(x)|$

Seu problema

Voltando à integral do curso de Konstantin e Yaroslav, você pode ver que ela se comporta exatamente da mesma maneira que o seno - integral como uma função do limite superior. Isso significa que você só precisa calcular os valores com . Abaixo está o gráfico de vários desses valores com .

$$I_{x_0}(\lambda) = 2\int_{0}^{x_0}\cos\left(\lambda\cos(x)\right)\operatorname{sinc}(x)dx$$ $x_0 = \pi\mathbb{N}+\frac{\pi}{2}$$\lambda = 12\pi$

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

Aqui você pode ver o resultado de outro método de aceleração. Reorganizo somas parciais da seguinte maneira e obtenho uma nova sequência que converge muito mais rapidamente. Esse truque também é útil se você deseja avaliar a integral com alta precisão.

$$S_N' = \frac{1}{2}(S_N+S_{N+1})$$ $S_N'$

O método de Ooura para integrais senoidais de Fourier funciona aqui, veja:

Ooura, Takuya e Masatake Mori, uma fórmula exponencial dupla e robusta para integrais do tipo Fourier. Jornal de matemática computacional e aplicada 112.1-2 (1999): 229-241.

Eu escrevi uma implementação desse algoritmo, mas nunca trabalhei para obtê-lo rapidamente (digamos, nós / pesos de cache), mas mesmo assim, estou obtendo resultados consistentes em tudo além da precisão de flutuação:

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

Aqui está o código:

#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";
}

Você realmente não consegue ver a diferença entre a quadratura e o assintótico, porque eles ficam um em cima do outro, exceto como : $\lambda\to 0$