Méthode d'intégration numérique d'une intégrale oscillatoire difficile

25

J'ai besoin d'évaluer numériquement l'intégrale ci-dessous:

0sinc(xr)rE(r)dr

où , et . Ici est la fonction de Bessel modifiée du second type. Dans mon cas particulier, j'ai , etxR+λ,κ,ν>0Kλ=0,00313κ=0,00825ν=0,33E(r)=r4(λκ2+r2)ν5/2Kν5/2(λκ2+r2)xR+λ,κ,ν>0Kλ=0.00313κ=0.00825ν=0.33 .

J'utilise MATLAB, et j'ai essayé les fonctions intégrées integralet quadgk, ce qui me donne beaucoup d'erreurs (voir ci-dessous). J'ai naturellement essayé de nombreuses autres choses aussi, comme l'intégration par parties et la sommation d'intégrales de à( k + 1 ) x πkxπ(k+1)xπ .

Alors, avez-vous des suggestions sur la méthode que je devrais essayer ensuite?

MISE À JOUR (questions ajoutées)
J'ai lu l'article @Pedro lié à, et je ne pense pas que ce soit trop difficile à comprendre. Cependant, j'ai quelques questions:

  • Serait-il correct d'utiliser comme éléments de baseψ kxkψk , dans la méthode univariée de Levin décrite?
  • Pourrais-je à la place simplement utiliser une méthode Filon, car la fréquence des oscillations est fixe?

Exemple de code
>> 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
la source
Qu'est-ce que dans votre intégrale? x
Pedro
Tout nombre positif et réel. Je viens de mettre à jour mon message.
torbonde
Si vous pouviez afficher du code et des erreurs, il n'est probablement pas trop difficile de résoudre la plupart d'entre eux. Bien sûr, veuillez d'abord lire attentivement l'erreur et voir si vous pouvez la faire disparaître par vous-même.
Dennis Jaheruddin
Je ferai un commentaire plus tard aujourd'hui avec du code et des erreurs. Ou demain.
torbonde
D'accord, j'ai donc oublié. Mais maintenant, j'ai mis à jour mon article avec un exemple (j'ai divisé l'intégrale en deux en calculant explicitement). sinc
torbonde

Réponses:

12

J'ai écrit mon propre intégrateur, quadccqui s'adapte nettement mieux que les intégrateurs Matlab avec des singularités, et fournit une estimation d'erreur plus fiable.

Pour l'utiliser pour votre problème, j'ai fait ce qui suit:

>> 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 fonction fest maintenant votre intégrale. Notez que je viens d'attribuer une ancienne valeur à x.

Afin de m'intégrer sur un domaine infini, j'applique une substitution de variables:

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

c'est-à-dire que l'intégration gde 0 à 1 devrait être la même chose que l'intégration fde 0 à . Différentes transformations peuvent produire des résultats de qualité différente: mathématiquement, toutes les transformations doivent donner le même résultat, mais des transformations différentes peuvent produire des s plus lisses ou plus facilement intégrables .g

J'appelle ensuite mon propre intégrateur, quadccqui peut gérer les NaNs aux deux extrémités:

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

Notez que l'estimation d'erreur est énorme, c'est quadcc-à- dire qu'elle n'a pas beaucoup confiance dans les résultats. En regardant la fonction, cependant, cela n'est pas surprenant car elle oscille à des valeurs de trois ordres de grandeur au-dessus de l'intégrale réelle. Encore une fois, l'utilisation d'une transformation d'intervalle différente peut produire de meilleurs résultats.

Vous pouvez également envisager des méthodes plus spécifiques telles que celle-ci . C'est un peu plus compliqué, mais certainement la bonne méthode pour ce type de problème.

Pedro
la source
Merci beaucoup. J'examinerai les différentes méthodes. Pour mes besoins, l'erreur n'a pas besoin d'être aussi petite que la norme en eq integral(1e-10 je pense), mais 1.7e + 07 est toujours vraiment, vraiment grande. Peut-être qu'une autre transformation fera du bien, comme vous le mentionnez.
torbonde
@ cimrg.joe: Notez que l'estimation d'erreur est une estimation de l'erreur absolue basée, entre autres, sur les valeurs absolues maximales de l'intégrande. Dans certains cas extrêmes, la valeur renvoyée peut en fait être tout à fait correcte. Si vous recherchez dix chiffres de précision, je vous recommande fortement d'utiliser les méthodes de type Levin que j'ai mentionnées à la fin de mon article.
Pedro
Je n'ai peut-être pas besoin de dix chiffres de précision, mais je pense que j'en ai besoin d'au moins cinq. Votre méthode peut-elle produire cela?
torbonde
La méthode ne peut pas garantir ce type de précision pour votre intégrale, car les valeurs à l'extrémité droite de l'intervalle sont supérieures de plusieurs ordres de grandeur à l'intégrale elle-même.
Pedro
11

Comme le souligne Pedro, les méthodes de type Levin sont les méthodes les mieux établies pour ce type de problèmes.

Avez-vous accès à Mathematica? Pour ce problème, Mathematica les détectera et les utilisera par défaut:

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.

Voici un tracé sur une plage de valeurs de x:

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

Tracé de x = 0,5 à x = 10

Vous pouvez également spécifier manuellement la méthode de type Levin particulière à appliquer, ce qui dans ce cas peut entraîner une légère amélioration des performances:

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.

Voir la documentation pour plus de détails sur les méthodes de type Levin dans Mathematica .

Andrew Moylan
la source
Malheureusement, je n'ai pas accès à Mathematica - seulement MATLAB. Je vais juste mettre à jour ma question avec quelques questions supplémentaires, concernant le papier @Pedro lié à.
torbonde
D'accord, comme vous dites, vous devrez vous contenter de Matlab. J'ajouterai une autre réponse à ce sujet.
Andrew Moylan
5

Si vous n'avez pas accès à Mathematica, vous pouvez écrire une méthode de type Levin (ou autre méthode oscillatoire spécialisée) dans Matlab comme le suggère Pedro.

Utilisez-vous la bibliothèque chebfun pour Matlab? Je viens d'apprendre qu'il contient une implémentation d'une méthode de base de type Levin ici . L'implémentation est écrite par Olver (l'un des experts du domaine de la quadrature oscillatoire). Il ne traite pas des singularités, de la subdivision adaptative, etc., mais il peut être juste ce dont vous avez besoin pour commencer.

Andrew Moylan
la source
J'ai moi-même pensé à implémenter une méthode Levin, mais je ne sais pas encore si je suis prêt à relever le défi. Je pense que j'ai besoin de mieux comprendre la méthode. Je pourrais peut-être en parler à mon conseiller. Quoi qu'il en soit, la raison pour laquelle j'ai posé des questions sur les méthodes Filon, c'est qu'elles semblent plus faciles à mettre en œuvre. Et comme je n'ai pas besoin d'une précision extrêmement élevée, mais cela fait partie de ma thèse, la difficulté pèse.
torbonde
J'ai jeté un œil à la bibliothèque chebfun (qui est impressionnante) et à l'exemple d'intégration de Levin. Mais je ne peux pas le faire fonctionner. J'ai en fait posté une question à ce sujet ici .
torbonde
0

La transformation recommandée par Pedro est une excellente idée. Avez-vous essayé de jouer avec les paramètres de la fonction "quadgk" de Matlab? Par exemple, en utilisant la transformation de Pedro, vous pouvez faire ce qui suit:
quadgk(f, 0.0+eps, 1.0-eps, 'AbsTol', eps, 'MaxIntervalCount', 100000)
Utiliser cela me donne une solution de:
-2184689.50220729
et ne prend que 0,8 seconde (en utilisant les valeurs mentionnées ci-dessus: x = 10)
Walter Gander et Walter Gautschi ont un document sur la quadrature adaptative avec Matlab code que vous pouvez également utiliser (lien ici )

Joel Vroom
la source