Nested NIntegrate avec FindRoot

Question

J'essaie d'intégrer numériquement une fonction à l'aide de NIntegrate imbriqué:

$$F(N,x,s)=\int_{-\infty}^s \int_{-\infty}^{+\infty} K(N,z',x,x') g_{x',s'} dx'ds' $$

où le noyau de l'intégration, $K(N,z'x,x')$, est une expression désordonnée définie dans le code mathématique ci-dessous, et $g_{x',s'}$ est un gaussien bi-varié défini par:

$$ g_{x,s'}=\frac{n}{2\pi\sigma_{x'}\sigma_{s'}}\exp\left({ -\frac{x'^2}{2\sigma_{x'}^2} }\right)\exp\left({ -\frac{s'^2}{2\sigma_{s'}^2} }\right).$$

La (les) partie (s) délicate (s) est que:

1. $z'$ dans le $K(N,z',x,x')$ doit être résolu pour l'utilisation numérique de FindRoot et aura un $s'$ dépendance.
2. La limite supérieure d'intégration au-dessus $ds'$ est une variable $s$.
3. Je soupçonne que le noyau est oscillatoire avec $N$ (noté "Kernel" dans le code ci-dessous) donc peut-être une moyenne du noyau sur $N$ peut être fait pour simplifier le noyau et éliminer $N$ si les intégrations s'avèrent trop longues.

À la fin, je voudrais une fonction, F (N, x, s), qui pourrait tracer à travers $s$ pour un donné $(N,x)$ valeurs ie Plot [F [a, b, s, {s, -1e-5,1e-5}].

(*Constants*)
e = -1.60217733*10^-19;
m = 9.109389699999999*10^-31;
epsilon = 8.854187817620391*10^-12;
re = 2.81794092*10^-15;
c = 2.99792458*10^8;
n = -10^-10/e;
KK = 1;
lw = 0.026;
kw = (2 Pi)/lw;
gamma = 4000/0.511;
beta = Sqrt[1 - 1/gamma^2];
sigmaS = 10^-5;
sigmaX = 30*10^-6;
coeff = n/(2 Pi*sigmaS*sigmaX) Exp[-(xprime^2/(2 sigmaX^2))]*
Exp[-(sprime^2/(2 sigmaS^2))];

(*Preliminary Equations*)
rs2 = {zprime, xprime + KK/(gamma*kw) Sin[kw*zprime], 0};
ro2 = {(NN + 10000)*lw, x + KK/(gamma*kw) Sin[kw*(NN + 10000)*lw], 0};

betas = {beta - KK^2/(2 gamma^2) Cos[kw*zprime]^2,KK/gamma Sin[kw*zprime], 0};
betao = {beta - KK^2/(2 gamma^2) Cos[kw*(NN + 10000)*lw]^2,KK/gamma Sin[kw*(NN + 10000)*lw], 0};

betaDot = {(c*KK^2*kw)/(2 gamma^2)Sin[2 kw*zprime], -((KK*c*kw)/gamma) Sin[kw*zprime], 0};

deltar2 = ro2 - rs2;
Rgam2 = Sqrt[deltar2[[1]]^2 + deltar2[[2]]^2];

Ec2 = (e/(4 Pi*epsilon)) (deltar2/Rgam2 - betas)/(gamma^2 Rgam2^2 (1 - (deltar2/Rgam2).betas)^3);
Erad2 = (e/(4 Pi*epsilon)) Cross[deltar2/Rgam2,Cross[deltar2/Rgam2 - betas, betaDot]]/(c*Rgam2*(1 - (deltar2/Rgam2).betas)^3);

sumElong = (Ec2[[1]] + Erad2[[1]]);
sumEtran = (Ec2[[2]] + Erad2[[2]]);

(*Numerical Functions*)

ZPRIME[NN_?NumericQ, x_?NumericQ, xprime_?NumericQ, s_?NumericQ, sprime_?NumericQ] := zprime /.FindRoot[s - sprime == (Sqrt[gamma^2 + KK^2] (EllipticE[kw*(NN + 10000)*lw,KK^2/(gamma^2 + KK^2)] - EllipticE[kw zprime, KK^2/(gamma^2 + KK^2)]))/(gamma kw) -beta Sqrt[((NN + 10000)*lw - zprime)^2 + (x - xprime + (KK Sin[kw *(NN + 10000)*lw])/(gamma kw) - (KK Sin[kw zprime])/(gamma kw))^2], {zprime, 0}]


Kernel = coeff re/gamma (sumElong*betao[[1]] + sumEtran*betao[[2]])/.{zprime -> ZPRIME[NN, x, xprime, s, sprime]};

FNxprimesprime[NN_?NumericQ, x_?NumericQ, xprime_?NumericQ, s_?NumericQ, sprime_?NumericQ]:= Kernel

FNsprime[NN_?NumericQ, x_?NumericQ, s_?NumericQ, sprime_?NumericQ] :=NIntegrate[FNxprimesprime[NN, x, xprime, s, sprime], {xprime, -300/10^6, 300/10^6}]

FN[NN_?NumericQ,x_?NumericQ, s_?NumericQ] := NIntegrate[FNsprime[NN,x, s, sprime], {sprime,-10^-4, s}]

lst1 = Table[{ss, FN[0,0, ss], PrecisionGoal -> 5] // Quiet}, {ss, -10^-5, 10^-5, 10^-6}]
ListPlot[lst1]

Answer 1

Ce cas n'est pas si différent de ce dont nous avons discuté ici . Nous avons utilisé la même approche et obtenir une réponse

(*Constants*)e = -1.60217733*10^-19;
m = 9.109389699999999*10^-31;
epsilon = 8.854187817620391*10^-12;
re = 2.81794092*10^-15;
c = 2.99792458*10^8;
n = -10^-10/e;
KK = 1;
lw = 0.026;
kw = (2 Pi)/lw;
gamma = 4000/0.511;
beta = Sqrt[1 - 1/gamma^2];
sigmaS = 10^-5;
sigmaX = 30*10^-6;
coeff = n/(2 Pi*sigmaS*sigmaX) Exp[-(xprime^2/(2 sigmaX^2))]*
   Exp[-(sprime^2/(2 sigmaS^2))];

(*Preliminary Equations*)
rs2 = {zprime, xprime + KK/(gamma*kw) Sin[kw*zprime], 0};
ro2 = {(NN + 10000)*lw, x + KK/(gamma*kw) Sin[kw*(NN + 10000)*lw], 0};

betas = {beta - KK^2/(2 gamma^2) Cos[kw*zprime]^2, 
   KK/gamma Sin[kw*zprime], 0};
betao = {beta - KK^2/(2 gamma^2) Cos[kw*(NN + 10000)*lw]^2, 
   KK/gamma Sin[kw*(NN + 10000)*lw], 0};

betaDot = {(c*KK^2*kw)/(2 gamma^2) Sin[
     2 kw*zprime], -((KK*c*kw)/gamma) Sin[kw*zprime], 0};

deltar2 = ro2 - rs2;
Rgam2 = Sqrt[deltar2[[1]]^2 + deltar2[[2]]^2];

Ec2 = (e/(4 Pi*epsilon)) (deltar2/Rgam2 - 
      betas)/(gamma^2 Rgam2^2 (1 - (deltar2/Rgam2).betas)^3);
Erad2 = (e/(4 Pi*epsilon)) Cross[deltar2/Rgam2, 
     Cross[deltar2/Rgam2 - betas, betaDot]]/(c*
      Rgam2*(1 - (deltar2/Rgam2).betas)^3);

sumElong = (Ec2[[1]] + Erad2[[1]]);
sumEtran = (Ec2[[2]] + Erad2[[2]]);

(*Numerical Functions*)

ZPRIME[NN_?NumericQ, x_?NumericQ, xprime_?NumericQ, s_?NumericQ, 
  sprime_?NumericQ] := 
 zprime /. 
  FindRoot[s - 
     sprime == (Sqrt[
         gamma^2 + 
          KK^2] (EllipticE[kw*(NN + 10000)*lw, 
           KK^2/(gamma^2 + KK^2)] - 
          EllipticE[kw zprime, KK^2/(gamma^2 + KK^2)]))/(gamma kw) - 
     beta Sqrt[((NN + 10000)*lw - zprime)^2 + (x - 
           xprime + (KK Sin[kw*(NN + 10000)*lw])/(gamma kw) - (KK Sin[
               kw zprime])/(gamma kw))^2], {zprime, 0}]

kernel = coeff re/gamma (sumElong*betao[[1]] + sumEtran*betao[[2]]);
FN[nn_?NumericQ, x0_?NumericQ, ss_?NumericQ] := 
 NIntegrate[
   Evaluate[(kernel /. 
       zprime -> ZPRIME[NN, x, xprime, s, sprime]) /. {NN -> nn, 
      x -> x0, xprime -> xp, s -> ss, sprime -> sp}], {xp, -300/10^6, 
    300/10^6}, {sp, -10^-4, ss}, PrecisionGoal -> 6] // Quiet

Maintenant, nous pouvons calculer et tracer des listes

lst1 = Table[{ss, FN[0., 0, ss]}, {ss, -10^-5, 10^-5, .5 10^-6}];
    lst2 = Table[{ss, FN[1., 0, ss]}, {ss, -10^-5, 10^-5, .5 10^-6}];
    ListLinePlot[{lst1, lst2}, PlotRange -> All,Frame -> True]