Método para la integración numérica de la integral oscilatoria difícil.

25

Necesito evaluar numéricamente la integral a continuación:

0sinc(xr)rE(r)dr

donde E(r)=r4(λκ2+r2)ν5/2Kν5/2(λκ2+r2),xR+yλ,κ,ν>0. AquíKes la función Bessel modificada del segundo tipo. En mi caso particular tengoλ=0.00313,κ=0.00825yν=0,33.

Estoy usando MATLAB, y he probado las funciones integradas integraly quadgk, lo que me da muchos errores (ver más abajo). Naturalmente, también he intentado muchas otras cosas, como la integración por partes y la suma de integrales de a ( k + 1 ) x π .kXπ(k+1)Xπ

Entonces, ¿tiene alguna sugerencia sobre qué método debería probar a continuación?

ACTUALIZACIÓN (preguntas adicionales)
Leí el documento al que se vinculó @Pedro, y no creo que fuera demasiado difícil de entender. Sin embargo, tengo algunas preguntas:

  • ¿Estaría bien usar como elementos básicos ψ k , en el método univariado de Levin descrito?Xkψk
  • ¿Podría usar un método Filon, ya que la frecuencia de las oscilaciones es fija?

Código de ejemplo
>> 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
fuente
¿Qué es en tu integral? X
Pedro
Cualquier número positivo, real. Acabo de actualizar mi publicación.
torbonde
Si pudieras mostrar algunos códigos y errores, probablemente no sea demasiado difícil resolver la mayoría de ellos. Por supuesto, primero intente leer el error detenidamente y vea si puede hacerlo desaparecer por su cuenta.
Dennis Jaheruddin el
Haré un comentario más tarde hoy con algunos códigos y errores. O mañana.
torbonde
Está bien, así que lo olvidé. Pero ahora actualicé mi publicación con un ejemplo (he dividido la integral en dos calculando explícitamente). syonortedo
torbonde

Respuestas:

12

He escrito mi propio integrador, quadccque maneja sustancialmente mejor que los integradores de Matlab con singularidades, y proporciona una estimación de error más confiable.

Para usarlo para su problema, hice lo siguiente:

>> 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) );

La función fahora es tu integrando. Tenga en cuenta que acabo de asignar cualquier valor anterior a x.

Para integrar en un dominio infinito, aplico una sustitución de variables:

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

es decir, integrar gde 0 a 1 debería ser lo mismo que integrar fde 0 a . Diferentes transformaciones pueden producir resultados de diferente calidad: matemáticamente todas las transformaciones deberían dar el mismo resultado, pero diferentes transformaciones pueden producir s más suaves o más fácilmente integrables .g

Luego llamo a mi propio integrador, quadccque puede manejar el NaNs en ambos extremos:

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

Tenga en cuenta que la estimación del error es enorme, es decir quadcc, no tiene mucha confianza en los resultados. Sin embargo, al observar la función, esto no es sorprendente, ya que oscila en valores de tres órdenes de magnitud por encima de la integral real. Nuevamente, usar una transformación de intervalo diferente puede producir mejores resultados.

También es posible que desee ver métodos más específicos como este . Es un poco más complicado, pero definitivamente es el método correcto para este tipo de problema.

Pedro
fuente
Muchas gracias. Echaré un vistazo a los diferentes métodos. Para mis propósitos, el error no necesita ser tan pequeño como es estándar en la ecuación integral(1e-10, creo), pero 1.7e + 07 sigue siendo muy, muy grande. Quizás otra transformación sirva, como mencionas.
torbonde
@ cimrg.joe: Tenga en cuenta que la estimación del error es una estimación del error absoluto basada, entre otros, en los valores absolutos máximos del integrando. En algunos casos extremos, el valor devuelto puede estar bastante bien. Si está buscando diez dígitos de precisión, le recomiendo usar los métodos de tipo Levin que mencioné al final de mi publicación.
Pedro
Tal vez no necesito diez dígitos de precisión, pero creo que necesito al menos cinco. ¿Puede tu método producir eso?
torbonde
El método no puede garantizar ese tipo de precisión para su integral, ya que los valores en el extremo derecho del intervalo son varios órdenes de magnitud más grandes que la integral en sí.
Pedro
11

Como señala Pedro, los métodos de tipo Levin son los mejores métodos establecidos para este tipo de problemas.

¿Tienes acceso a Mathematica? Para este problema, Mathematica los detectará y usará por defecto:

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.

Aquí hay una gráfica sobre un rango de valores de x:

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

Trazar desde x = 0.5 a x = 10

También puede especificar manualmente el método particular de tipo Levin para aplicar, que en este caso puede producir una ligera mejora en el rendimiento:

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.

Consulte la documentación para obtener detalles sobre los métodos de tipo Levin en Mathematica .

Andrew Moylan
fuente
Lamentablemente no tengo acceso a Mathematica, solo a MATLAB. Solo actualizaré mi pregunta con algunas preguntas adicionales, en relación con el documento @Pedro vinculado.
torbonde
OK, como dices, tendrás que conformarte con Matlab. Agregaré otra respuesta sobre eso.
Andrew Moylan
5

Si no tiene acceso a Mathematica, podría escribir un método de tipo Levin (u otro oscilatorio especializado) en Matlab como sugiere Pedro.

¿Utiliza la biblioteca chebfun para Matlab? Me acabo de enterar que contiene una implementación de un método básico de tipo Levin aquí . La implementación está escrita por Olver (uno de los expertos en el campo de la cuadratura oscilatoria). No trata con singularidades, subdivisión adaptativa, etc., pero puede ser justo lo que necesita para comenzar.

Andrew Moylan
fuente
He pensado en implementar un método de Levin, pero no estoy seguro de si estoy preparado para el desafío todavía. Creo que necesito entender el método un poco mejor. Tal vez podría hablar con mi asesor al respecto. De todos modos, la razón por la que pregunté sobre los métodos de Filon, es que parecen más fáciles de implementar. Y como no necesito una precisión extremadamente alta, pero esto es parte de mi tesis de maestría, la dificultad pesa.
torbonde
He echado un vistazo a la biblioteca chebfun (que es impresionante) y al ejemplo de integración de Levin. Pero no puedo hacerlo funcionar. De hecho, he publicado una pregunta al respecto aquí .
torbonde
0

La transformación recomendada por Pedro es una gran idea. ¿Has intentado jugar con los parámetros en la función "quadgk" de Matlab? Por ejemplo, usando la transformación de Pedro, puede hacer lo siguiente:
quadgk(f, 0.0+eps, 1.0-eps, 'AbsTol', eps, 'MaxIntervalCount', 100000)
Usar esto me da una solución de:
-2184689.50220729
y solo toma 0.8 segundos (usando los valores mencionados anteriormente: x = 10)
Walter Gander y Walter Gautschi tienen un documento sobre la cuadratura adaptativa con Matlab código que puede usar también (enlace aquí )

Joel Vroom
fuente