0
$\begingroup$

I am trying to numerically integrate a function using nested NIntegrate:

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

where the kernel of the integration, $K(N,z'x,x')$, is a messy expression defined in the mathematica code below, and $g_{x',s'}$ is a bi-variate gaussian defined by:

$$ 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).$$

The tricky part(s) is that:

  1. $z'$ in the $K(N,z',x,x')$ needs to be solved for numerically using FindRoot and will have a $s'$ dependence.
  2. The integration upper limit over $ds'$ is a variable $s$.
  3. I suspect the kernel is oscillatory with $N$ (denoted "Kernel" in the code below) so maybe an averaging of the kernel over $N$ can be done to simplify the kernel and eliminate $N$ if the integrations prove to be too time consuming.

At the end, I would like a function, F(N,x,s), that would be able to plot across $s$ for a given $(N,x)$ values i.e. 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]
$\endgroup$
0
$\begingroup$

This case is not so different from what we discussed here. We have used same approach and get answer

(*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

Now we can calculate and plot lists

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]

Figure 1

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.