Evaluación numérica de integral altamente oscilatoria

11

En este curso avanzado sobre aplicaciones de la teoría de funciones complejas en un punto de un ejercicio, la integral altamente oscilatoria

I(λ)=cos(λcosx)sinxxdx

tiene que aproximarse para valores grandes de λ utilizando el método del punto de silla de montar en el plano complejo.

Debido a su naturaleza altamente oscilatoria, esta integral es muy difícil de evaluar utilizando la mayoría de los otros métodos. Estos son dos fragmentos de la gráfica del integrando para λ=10 a diferentes escalas:

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

Una aproximación asintótica de orden principal es

I1(λ)=cos(λ14π)2πλ

y un refinamiento adicional (mucho más pequeño) agrega el término

I2(λ)=18sin(λ14π)2πλ3

Un gráfico de los valores aproximados en función de λ ve de la siguiente manera:

I (lambda) aprox.

Ahora viene mi pregunta: para ver visualmente qué tan buena es la aproximación, me gustaría compararla con el "valor real" de la integral, o más precisamente con una buena aproximación a la misma integral usando un algoritmo independiente. Debido a lo pequeño de la corrección de subtítulos, esperaría que esto sea muy cercano.

λtanh(sinh)

mpmath aprox

Finalmente probé suerte con un integrador de Montecarlo usando una muestra de importancia que implementé, pero tampoco pude obtener resultados estables.

λ>1

doetoe
fuente
¿Es la función par?
nicoguaro
Sí, es par
doetoe
¿Has intentado convertir tu integral en una ODE?
nicoguaro
1
x
1
λλxx(cos(λxcosx)sincx)

Respuestas:

12

Use el teorema de Plancherel para evaluar esta integral.

f,g

I=f(x)g(x)dx=F(k)G(k)dk

F,Gf,gsinx/xrect(k)cos(λcosx)λ|Jn(x)|n>|x|

πJ0(λ)[0,2π]

smh
fuente
Gracias, esta es una muy buena idea!
doetoe
7

πN+π2

Asintóticas

I(λ)2πλ[cos(λπ4)+c1sin(λπ4)λ+c2cos(λπ4)λ2+c3sin(λπ4)λ3+]
c1=18

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, obtienes un seno bastante agradable que coincide con el que obtuviste anteriormente.

18 años

Si desea encontrar los siguientes coeficientes, un código un poco más sofisticado si es necesario. La idea del siguiente código es tomar varios valores límite superiores altos y "promediar" sus 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]

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

Explicación

Ejemplo simple

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

seno

S(x)

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

Tu problema

Volviendo a la integral del curso de Konstantin y Yaroslav, puede ver que se comporta exactamente de la misma manera que el seno: integral en función del límite superior. Eso significa que solo necesita calcular los valores con . A continuación se muestra la gráfica de varios de estos valores con .

Ix0(λ)=20x0cos(λcos(x))sinc(x)dx
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}]

acc

Aquí puede ver el resultado de otro método de aceleración. sumas parciales de la siguiente manera y nueva secuencia que converge mucho más rápido. Ese truco también resulta útil si desea evaluar integral con alta precisión.

SN=12(SN+SN+1)
SN

David Saykin
fuente
¡Agradable! ¿Son los instructores del curso sus profesores de la vida real? Su curso es fantástico, aunque muy duro y rápido
doetoe
@doetoe sí, soy estudiante de Konstantin. Compartió conmigo un enlace a su pregunta.
David Saykin
6

El método de Ooura para las integrales seno de Fourier funciona aquí, ver:

Ooura, Takuya y Masatake Mori, una robusta fórmula doble exponencial para integrales de tipo Fourier. Revista de matemática computacional y aplicada 112.1-2 (1999): 229-241.

Escribí una implementación de este algoritmo, pero nunca me puse a trabajar para hacerlo rápido (por ejemplo, nodos de almacenamiento en caché / pesos), pero no obstante, obtengo resultados consistentes en todo más allá de la precisión de flotación:

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

Aquí está el 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";
}

Realmente no se puede ver la diferencia entre la cuadratura y el asintótico porque se encuentran uno encima del otro, excepto como : λ0ingrese la descripción de la imagen aquí

usuario14717
fuente
Gracias, esto es realmente bueno! Todavía no lo hice funcionar, mi instalación de impulso no es compatible, pero ahora estoy descargando la última versión.
doetoe
Solo para estar seguro: en 23 tienes cos (lambda * cos (x)) / x, sin el factor sin (x) del integrando. ¿Es ooura_fourier_sin el que supone que este factor sin (x) multiplica el integrando que se le pasa?
doetoe
Lo tengo funcionando. Parece que todas sus dependencias son solo encabezado, por lo que ni siquiera tuve que instalarlo ni compilarlo (excepto el ejecutable). ¡Espero que se incluya en boost!
doetoe
@doetoe: Sí, el factor está implícito e incorporado en los pesos de la cuadratura. Estoy considerando acelerarlo refactorizándolo para almacenar en caché los nodos y los pesos; ¡ya veremos! sin(x)
user14717
@doetoe: se fusionó con Boost 1.71. La API es un poco diferente a la que da esta respuesta.
user14717