# Numerical integration with NIntegrate - faster or better method?

Long story short: I need to get a function ($f(r,z,\alpha)$) via numerical integration, and a number from another numerical integration for this function. The method I have at the moment does converse & gives an answer, it is based on NIntegrate but takes a long time.
I am looking for suggestions for faster method. Even less precise ones, as I just a rough number.

I found this very tempting method but could not generalise is for multivariables and tables, like I will need to do for the three arguments of the function.

I am calculating the intensity of a laser beam going through a particular optical component, an $axicon$. I found the data from this paper.

The electric field is given by:

$$E(r,z) = \frac{2\pi}{\lambda z}w_0 w(z) \int_0^{\infty} e^{-x^2}e^{i\pi(ax^2 - \frac{x}{\Delta})}J_0(\beta x)x\, \mathrm{d}x,$$ where $\beta= \beta(r,z), \, a=a(z), \,\Delta = \Delta(\alpha,z), \,$.
I am writing this as:

Efield[r_?NumericQ, z_?NumericQ, α_?NumericQ] := (
2 π)/(λlaser z)*w0*w[z]*
NIntegrate[
Exp[-x^2]*
Exp[I*π*(a [z] x^2 - x/Δ[α, z])]*
BesselJ[0, β [r, z] x]*x, {x, 0, ∞}]


The intensity of the beam is given by: $$I(r,z) = |E(r,z)|^2,$$

and I am writing as:

Intensity[r_?NumericQ, z_?NumericQ, α_?NumericQ] := Abs[Efield[r, z, α]]^2


Now I need the power $P$:

$$P = \int^{\infty}_0 2 \pi I(r,\mathrm{distance})r\,\mathrm{d}dr,$$ where the $z$ coordinate does not matter as it should give the name answer $\forall z$.

I am writing this as:

PowerLens = NIntegrate[(Intensity[r, distance, conv*20]) * r*2 π, {r, 0, 1}]

• I need to look at the intensity (vs. $r$) for different values of $\alpha$ and $z$. I can do this with Table but it would take an NIntegrate each time, hence the long time.

• The power needs two NIntegrate, so takes a long time too.

Full code:

nm = 10^-6 mm;
mm = 10^-3;
μm = 10^-3 mm;
λlaser = 532 nm;
klaser = (2 π)/λlaser;
w0 = 5 mm;
flens = 50 mm;
ng = 1.5;
conv = (2 π)/360;
Gaussian beam parameters

zR = π w0^2/λlaser // N
w[z_] = w0*Sqrt[1 + z^2/zR^2]
Rc[z_] = z*(1 + zR^2/z^2)
Guoy[z_] = ArcTan[z/zR]

GaussianBeamField[ρ_, z_] = (w0/w[z])*Exp[-ρ^2/w[z]^2]*
Exp[I*(klaser z - Guoy[z] + (klaser ρ^2)/(2 Rc[z]))]

GaussianBeamIntensity[ρ_, z_] =
Abs[GaussianBeamField[ρ, z] ]^2 // ComplexExpand

Normalisation[z_] :=
Integrate[2 π*ρ*GaussianBeamIntensity[ρ, z], {ρ,
0, ∞}, Assumptions -> {Element[{ρ, z}, Reals]}]

GaussianBeamIntensityNormalised[ρ_, z_] = (1/Normalisation[0])*
GaussianBeamIntensity[ρ, z]

b[α_] = (2 π (ng - 1) )/λlaser Tan[α]
Δ[α_, z_] = π/(b[α]*w[z])
β[r_, z_] = (2 π w[z] r)/(λlaser z)
a[z_] = w[z]^2/λlaser*(1/z + (1/Rc[z] - 1/flens))

Efield[r_?NumericQ, z_?NumericQ, α_?NumericQ] := (
2 π)/(λlaser z)*w0*w[z]*
NIntegrate[ Exp[-x^2]*Exp[I*π*(a [z] x^2 - x/Δ[α, z])]*
BesselJ[0, β [r, z] x]*x, {x, 0, ∞}]

Intensity[r_?NumericQ, z_?NumericQ, α_?NumericQ] := Abs[Efield[r, z, α]]^2

p1 = Plot[GaussianBeamIntensityNormalised[r, flens], {r, 0, 15 mm}]

p2 = Plot[IntensityNormalised[r, flens, 20*conv], {r, 0, 15 mm}]

Show[p1, p2]

PowerLens = NIntegrate[Intensity[r, flens, conv*20]*r*2 π, {r, 0, ∞}]

• your intensity appears to have a single sharp peak over a very narrow range of r and is zero everywhere else. You'll get better/faster results of you can integrate over a smaller relevant domain instead of 0,Infinity – george2079 Jan 23 '18 at 22:40
• In a hurry, I suggest you use the Method option in NIntegrate. I remember that a fast method exists for integrands involving Bessel functions. Another posibility could be to implement the numerical Bessel Transform... – José Antonio Díaz Navas Jan 23 '18 at 22:42
• @george2079 yes that would be ideal, the thing is I don’t always know where that peak is located – SuperCiocia Jan 23 '18 at 23:11