Intégration numérique par la méthode double-exponentielle

• f(x) < M1 e^{α\abs{x}^β} sur [0,\infty[ ;
• f(z) < M2 sur un cône d'ouverture tau ;

-> integration_SE(f,M1,alpha,beta,tau,M2)

gp > f(x)=sqrt(x)/gamma(x+2);

$f$ a une décroissance en $\exp(-x \log x)$ en l'infini et un comportement en $\sqrt{x}$ en 0. On effectue un changement $φ(t)=\exp(t-α \exp(-β t))$ avec $α=1.5$ et $β=1$.

gp > \p100
   realprecision = 115 significant digits (100 digits displayed)
gp > \g1
   debug = 1
gp > alpha=1.5; beta=1;
gp > M1=7 \\ ploth(x=0,10,f(x)*exp(alpha*x^beta))
% 7
gp > tau = 0.78; \\ ~Pi/4;
gp > M2 = 3 \\ plot_Ztau_exp(x->abs(f(x)),tau,alpha,beta)
% 3
gp > integration_Rplus(f,M1,alpha,beta,tau,M2)
  t=tau-3.27 e-3
  Ct=3.96 e0
  h=2.054 e-2 [pas_h_DE1]
  npoints = 248 [Rplus]
  [intnum] precision lost = 3 digits
% 1.026308986088920438152866878884727476112511603518869085278521223463073931130358019672349854444849701

on peut accélérer en prenant une plus grande valeur de tau

gp > tau_bis = 1.26;
gp > M2_bis = 1.5e7 \\ plot_Ztau_exp(x->abs(f(x)),tau_bis,alpha,beta)
% 15000000.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
gp > integration_Rplus(f,M1,alpha,beta,tau_bis,M2_bis)
  t=tau-4.97 e-3
  Ct=3.85 e0
  h=3.118 e-2 [pas_h_DE1]
  npoints = 163 [Rplus]
  [intnum] precision lost = 3 digits
% 1.026308986088920438152866878884727476112511603518869085278521223463073931130358019672349854444849701

le facteur de proportionalité est clair

gp > \p500
   realprecision = 500 significant digits
gp > alpha=1.5;tau=0.78;tau_bis=1.26;M2_bis=1.5e7; \\ redéfinir pour avoir la bonne précision
gp > I500=integration_Rplus(f,M1,alpha,beta,tau,M2)
  t=tau-6.72 e-4
  Ct=5.55 e0
  h=4.220 e-3 [pas_h_DE1]
  npoints = 1577 [Rplus]
  [intnum] precision lost = 4 digits
% 1.0263089860889204381528668788847274761125116035188690852785212234630739311303580196723498544448497007873580061101556090579147164110354872779716507759227768380041446579741063269810539725908402784320199013115395057988286866959879902621200365381274838898030238643533495561773416132700037429354067265880043838022564805724067163572295473124777309413177179979553218113399774640786976773455824091228001014980957959541024927737901404811585250721105019066538379950368147446768838942830952347398267322344287769
gp > integration_Rplus(f,M1,alpha,beta,tau_bis,M2_bis)
  t=tau-1.07 e-3
  Ct=5.38 e0
  h=6.729 e-3 [pas_h_DE1]
  npoints = 989 [Rplus]
  [intnum] precision lost = 3 digits
% 1.0263089860889204381528668788847274761125116035188690852785212234630739311303580196723498544448497007873580061101556090579147164110354872779716507759227768380041446579741063269810539725908402784320199013115395057988286866959879902621200365381274838898030238643533495561773416132700037429354067265880043838022564805724067163572295473124777309413177179979553218113399774640786976773455824091228001014980957959541024927737901404811585250721105019066538379950368147446768838942830952347398267322344287768

la borne de troncature est assez bonne

gp > integration_exp_until(f,1.5,1,6.729e-3)
  npoints = 866 [exp]
  0
  [intnum] precision lost = 3 digits
% 1.0263089860889204381528668788847274761125116035188690852785212234630739311303580196723498544448497007873580061101556090579147164110354872779716507759227768380041446579741063269810539725908402784320199013115395057988286866959879902621200365381274838898030238643533495561773416132700037429354067265880043838022564805724067163572295473124777309413177179979553218113399774640786976773455824091228001014980957959541024927737901404811585250721105019066538379950368147446768838942830952347398267322344287776

et le pas d'intégration également assez précis

gp > I500-integration_exp_until(f,1.5,1,7e-3)
  npoints = 833 [exp]
  0
  [intnum] precision lost = 3 digits
% 1.0959773805847659286 E-494

• f(x) < M1 e^{α\abs{x}^β} sur \R ;
• f(z) < M2 sur un cône d'ouverture tau ;

-> integration_SE(f,M1,alpha,beta,tau,M2)

Gaussienne et fonctions de Hermite

gp > \g0
   debug = 0
gp > \p 1000
   realprecision = 1001 significant digits (1000 digits displayed)
gp > f(z)=1/(1+z^2);

paramètres d'intégration pour 1000 chiffres (h=0.004;n=2200)

gp > D=precision(1.)*log(10);
gp > tau=Pi/2.1;
gp > M1=1; u=1; M2=2; \\ plot_Ztau_shsh(f,tau,u)
gp > h=pas_h_shsh(D,tau,M2,u);show(h)
% 0.004071096256146010121
gp > npoints=ceil(tmax_shsh(D,M1,/*alpha=*/2)/h)
% 2073

calcul de référence

gp > Pi-integration_shsh(f,h,npoints,/*lambda=*/1)
% -5.946130296599909037 E-1002

décalage horizontal : le résultat devient très mauvais du fait des pôles

gp > diff=Pi-integration_shsh(z->f(z-15),h,npoints,1);show(diff)
% 1.1644390926442540736 E-12

on retrouve un résultat correct en choisissant tau assez petit pour éviter les pôles.

gp > show(asinh(asinh(10+I)))
% 1.8198999017302092895 + 0.03132601926270185771*I
gp > tau_dec=0.03;
gp > M2_dec = 2; \\ plot_Ztau_shsh(x->f(x-15),tau,u)
gp > h_dec=pas_h_shsh(D,tau_dec,M2_dec,u);show(h)
% 0.004071096256146010121

Le nombre de points est alors considérable :

gp > npoints_dec=ceil(asinh(asinh(2*10^1000+10))/h_dec)
% 103211
gp > \\Pi-integration_shsh(z->f(z-10),h_dec,npoints_dec,1)

mais la meilleure solution reste de calculer les termes dûs aux pôles, qui corrigent exactement l'erreur

gp > Rho=[15+I,15-I]; Res=[-I/2,I/2];
gp > diff+poles_shsh(Rho,Res,h)
% -6.673167885609415581 E-1002 + 0.E-1014*I

La fonction $f(x)=e^{-\ch(x)+i\ch(2x)}$ a une décroissance double exponentielle, mais oscille plus qu'elle ne décroît.

gp > f(z)=exp(-cosh(z)+I*cosh(2*z));
gp > \p50
   realprecision = 57 significant digits (50 digits displayed)
gp > f(5.6) \\ nh=5.6 convient
% -1.8444123454949663188659922688265233544039431479976 E-59 - 4.2025509810487181931181431963759407647201701526148 E-60*I
gp > my(h=0.01);integration(f,h,5.6/h)
% -0.0078972187642165579381310641669378229064813703996416 + 0.42873574152985099568629064252831104846577365576273*I

pour obtenir 50 chiffres, on a besoin d'un pas très petit

gp > #
   timer = 1 (on)
gp > my(h=0.0001);integration(f,h,5.6/h)
time = 3,840 ms.
% -0.0078977140259439402764638281328125967392134231974354 + 0.42873579344423196289801335397613892988520816667985*I

Les oscillations sont due au fait que \R est la limite entre la zone où f est très fortement croissante et celle où elle est fortement décroissante. Par un changement de variable, on peut déplacer le chemin pour se situer toujours dans les zones de forte décroissance : en +iπ/4 pour x>0, et en -iπ/4 pour x<0. On le fait par exemple avec la fonction tanh.

gp > g(t)=f(t+I*Pi/4*tanh(t))*(1+I*Pi/(4*cosh(t)^2))
time = 0 ms.
% (t)->f(t+I*Pi/4*tanh(t))*(1+I*Pi/(4*cosh(t)^2))
gp > g(2.8) \\ la décroissance est encore plus rapide
time = 0 ms.
% -2.7363507743417357988597804623217568122612654463066 E-62 + 4.6704776781954731984268092639742316768800602057873 E-62*I

l'intégration est considérablement plus efficace

gp > my(h=0.04);integration(g,h,2.8/h)
time = 0 ms.
% -0.0078977140259439402764638281328125967392134231974354 + 0.42873579344423196289801335397613892988520816667985*I

$f(x)=\Gamma(1+ix)$ n'est pas bornée sur un cône à cause de sa croissance exponentielle à droite.

gp > f(x)=gamma(1+I*x);
gp > ref = 2*Pi*exp(-1); \\ = \int_\R f
gp > \\ voir avec plot_Ztau_sh(f,Pi/10)
gp > D = precision(1.)*log(10);

l'intégration double exponentielle ne fonctionne pas

gp > integration_sh_until(f,pas_h_sh(D,Pi/5,5))/ref
time = 60 ms.
% 1.0000000000000000000000000000000000320635836429609 + 0.E-79*I

une étude précise donnerait ce pas d'intégration

gp > h = Pi^2/(D*log(D+1))
time = 0 ms.
% 0.015394784482143166516588760072417670356156469470832
gp > integration_sh_until(f,h)/ref
time = 110 ms.
% 1.0000000000000000000000000000000000000000000000000 + 0.E-80*I

un décalage vers la gauche permet de résoudre le problème

gp > alpha=Pi/2; beta=1; tau=Pi/4;
gp > M1 = 2;
gp > h = pas_h_sh(,Pi/4,M1)
time = 0 ms.
% 0.040480070124782407726623398265139794428475517767974
gp > integration_shteta_until(f,Pi/4,h)/ref
time = 30 ms.
% 1.0000000000000000000000000000000000000000000000000 + 2.212207383793983914 E-59*I