Évaluation numérique de l'intégrale fortement oscillatoire

11

Dans ce cours avancé sur les applications de la théorie des fonctions complexes à un moment donné dans un exercice l'intégrale hautement oscillatoire

je(λ)=-cos(λcosX)péchéXXX

doit être approximée pour les grandes valeurs de utilisant la méthode du point de selle dans le plan complexe.λ

En raison de sa nature hautement oscillatoire, cette intégrale est très difficile à évaluer en utilisant la plupart des autres méthodes. Ce sont deux fragments du graphe de l'intégrande pour à différentes échelles:λ=dix

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

Une approximation asymptotique d'ordre supérieur est

je1(λ)=cos(λ-14π)2πλ

et un autre raffinement (beaucoup plus petit) ajoute le terme

je2(λ)=18péché(λ-14π)2πλ3

Un graphique des valeurs approximatives en fonction de se présente comme suit:λ

I (lambda) environ

Vient maintenant ma question: pour voir visuellement à quel point l'approximation est bonne, je voudrais la comparer à la "valeur réelle" de l'intégrale, ou plus précisément à une bonne approximation de la même intégrale en utilisant un algorithme indépendant. En raison de la petitesse de la correction de sous-reliure, je m'attendrais à ce que ce soit très proche.

J'ai essayé d'évaluer l'intégrale pour certains utilisant d'autres algorithmes, mais avec très peu de succès: Mathematica et Matlab utilisant l'intégrateur numérique par défaut ne parviennent pas à produire une valeur significative (et le rapportent explicitement), mpmath utilisant à la fois le doublement exponentiel substitution et la méthode de Gauss-Legendre produisent des résultats très bruyants, bien qu'elle ait une légère tendance à osciller autour des valeurs fournies par la méthode du point de selle, comme le montre ce graphique:λtanh(sinh)

mpmath environ

Enfin, j'ai tenté ma chance avec un intégrateur Monte-Carlo en utilisant un échantillon d'importance que j'ai implémenté, mais je n'ai pas réussi à obtenir de résultats stables non plus.

Quelqu'un at-il une idée de la façon dont cette intégrale pourrait être évaluée indépendamment pour une valeur fixe de ou plus?λ>1

doetoe
la source
La fonction est-elle uniforme?
nicoguaro
Oui, c'est même
doetoe
Avez-vous essayé de transformer votre intégrale en ODE?
nicoguaro
1
Non, différencier wrt puis résoudre numériquement l'équation différentielle. X
nicoguaro
1
Votre premier tracé semble montrer une fonction différente de celle de votre intégrale. À savoir, il semble que remplacé par . C'est-à-dire que l'intrigue est de la fonctionλ x x ( cos ( λ x cos x ) sinc x )λλXX(cos(λXcosX)sincX)
Ruslan

Réponses:

12

Utilisez le théorème de Plancherel pour évaluer cette intégrale.

L'idée de base est que pour deux fonctions ,F,g

je=-F(X)g(X)X=-F(k)g(k)k

où sont les transformées de Fourier de . Vos fonctions ont toutes deux un support relativement faible dans le domaine spectral. Ici, et devraient avoir une transformée de Fourier (ou série) analytique, comme l' expansion de Jacobi-Anger . Vous pouvez tronquer la série infinie à environ raison de la décroissance super-exponentielle de la fonction de Besselpour. J'espère que cela t'aides.F,gF,gpéchéX/Xrect(k)cos(λcosX)λ | J n ( x ) | n > | x |λ|Jn(X)|n>|X|

Edit : En fait, vous devez utiliser les représentations de la série de Fourier ici au lieu des transformations. Le chemin de transformation conduit à dériver la représentation asymptotique que vous avez déjà (il s’agit juste de ). Le théorème de Plancherel ci-dessus fonctionne également pour les séries de Fourier avec un domaine d'intégration de sur la dernière intégrale.πJ0(λ)[0,2π]

smh
la source
Merci, c'est une très bonne idée!
doetoe
7

πN+π2

Asymptotiques

je(λ)2πλ[cos(λ-π4)+c1péché(λ-π4)λ+c2cos(λ-π4)λ2+c3péché(λ-π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}]

En sortie, vous obtenez un joli sinus qui coïncide avec celui que vous avez dérivé ci-dessus.

18

Si vous voulez trouver les coefficients suivants, un morceau de code un peu plus sophistiqué si nécessaire. L'idée du code ci-dessous est de prendre plusieurs valeurs limites hautes et de "moyenne" leurs résultats.

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,

Explication

Exemple simple

S(X)=0Xpéché(y)yy.
S()=π2

sinus

S(X)

SN=n=1N(-1)nn.
SSN+12(-1)N+1N+1.
S(X)0πN+π2péchéXXX
max|S(X)|

Ton problème

jeX0(λ)=20X0cos(λcos(X))sinc(X)X
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

SN=12(SN+SN+1)
SN

David Saykin
la source
Agréable! Les instructeurs du cours sont-ils vos professeurs réels? Leur parcours est fantastique, bien que très difficile et rapide
doetoe
@doetoe oui, je suis un étudiant de Konstantin. Il a partagé avec moi un lien vers votre question.
David Saykin
6

La méthode d'Ooura pour les intégrales de sinus de Fourier fonctionne ici, voir:

Ooura, Takuya et Masatake Mori, Une formule exponentielle double robuste pour les intégrales de type Fourier. Journal des mathématiques computationnelles et appliquées 112.1-2 (1999): 229-241.

J'ai écrit une implémentation de cet algorithme mais je n'ai jamais travaillé pour l'obtenir rapidement (par exemple, en mettant en cache les nœuds / poids), mais néanmoins, j'obtiens des résultats cohérents pour tout ce qui dépasse la précision flottante:

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

Voici le code:

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

λ0entrez la description de l'image ici

user14717
la source
Merci, c'est vraiment sympa! Je ne l'ai pas encore fait fonctionner, mon installation boost n'est pas compatible, mais je télécharge la dernière version maintenant.
doetoe
Juste pour être sûr: en 23, vous avez cos (lambda * cos (x)) / x, sans le facteur sin (x) de l'intégrande. Est-ce ooura_fourier_sin qui suppose ce facteur sin (x) pour multiplier l'intégrande qui lui est passée?
doetoe
Je l'ai fait fonctionner. Il et ses dépendances ne semblent être que des en-têtes, donc je n'ai même pas eu à installer ou à compiler (sauf pour l'exécutable). J'espère que cela sera inclus dans le boost!
doetoe
péché(X)
@doetoe: Il a été fusionné dans Boost 1.71. L'API est un peu différente de celle donnée par cette réponse.
user14717