# Weighted mean of complex exponential function using NIntegrate

I have defined the following functions:

γ[r_, v_, rDet_] :=
Which[
Abs[r - v tDet] >= rDet, 0,
r + v tDet <= rDet, π,
True, ArcCos[((v tDet)^2 + r^2 - rDet^2)/(2 v tDet r)]];
ρ[r_, v_, v0_] := Exp[-(1/2) (r/rCloud)^2] r Exp[-(1/2) (v/v0)^2]  v;
ΔΦ[v_] := (2 π )/λ c (4 (v*t)^2)/(rBeam)^2;


I want to weight ΔΦ by ρ and γ. I do so by defining a function Awfc that numerically integrates the product of the three functions over v and r and divides by the numerical integral over the two weighting functions:

Awfc[rDet_?NumericQ, v0_?NumericQ] :=
NIntegrate[ γ[r,v,rDet]*ρ[r,v,v0]* Exp[I ΔΦ[v]], {r,0,∞}, {v,0,∞ }}]/
NIntegrate[ γ[r,v,rDet]*ρ[r,v,v0], {r,0,∞}, {v,0,∞}}];


I calculate Awfc for some parameters rDet and v0

AwfcTable =
ParallelTable[
{rDet, Awfc[rDet, 0.5 v0], Awfc[rDet, v0], Awfc[rDet,2 v0], Awfc[rDet,3.7 v0]},
{rDet,0.0005,0.010,0.0005}];


using these values for the other parameters:

tDet = 0.7;
t = 0.230;
rCloud = 0.0025;
λ =  780 10 ^-9;
c = λ/20;
v0 = 0.00588;
rDet = 0.008/2;
rBeam = 0.015;


All kernels throw warnings:

NIntegrate::slwcon: Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small.

I think that ΔΦ is a hightly osciallating function so I don't see a way to fix this problem.

I use TableForm[Abs[AwfcTable]] to get the amplitude of the complex numbers in AwfcTable. The problem with the output is, that this amplitude is significantly larger than 1, especially for larger v0 and small rDet. The function ΔΦ I am calculating the weighted mean of, has a maximum amplitude of 1 so I think the weighted mean should also have a maximum of 1. (The complex phase is consistent with what I expect from a different approach to the problem I did using MATLAB)

Is there something wrong with my reasoning (meaning that the result Mathematica gives is actually correct), or is does the fault lie with my implementation of the problem?

-

Observation:

First of all I think it is always a useful trick to plot your problematic integrand if possible. It gives us often the clue in case NIntegrate complained about the particular integrand. If we can track down the issue we can often come up with a remedy. Given the following input if one sweeps over the rDet we get the following plots. Given your the quantity of interest which is ratio of two integrals.

Awfc[rDet_?NumericQ, v0_?NumericQ] := NIntegrate[ γ[r,v,rDet]*ρ[r,v,v0]* Exp[I ΔΦ[v]],
{r,0,∞}, {v,0,∞ }}]/NIntegrate[ γ[r,v,rDet]*ρ[r,v,v0], {r,0,∞}, {v,0,∞}}];


Conclusion:

The numerator integral involves a complex function that rapidly becomes zero. The denominator is a real function that is also fast asymptotic. Hence for faster evaluation of your ParallelTable we can NIntegrate on a shorter interval. But in this case we can also do the full integration by tuning up some options. One can also scale the integrals beforehand as that will avoid numerical error to some extent.

AwfcInf[rDet_?NumericQ, v0_?NumericQ] :=NIntegrate[γ[r, v, rDet]*ρ[r, v, v0]*Exp[I ΔΦ[v]],
{r, 0, Infinity}, {v, 0,Infinity}, AccuracyGoal -> 12,PrecisionGoal -> 30]/
(Re@NIntegrate[γ[r, v, rDet]*ρ[r, v, v0], {r, 0,Infinity}, {v, 0, Infinity},
AccuracyGoal -> 12,PrecisionGoal -> 30])


Now you can call your table safely.

AwfcTable = ParallelTable[{rDet, AwfcInf[rDet, 0.5 v0], AwfcInf[rDet, v0],
AwfcInf[rDet, 2 v0],AwfcInf[rDet, 3.7 v0]},
{rDet, 0.0005,0.010, 0.0005}]; // AbsoluteTiming


{11.972305, Null}

Ans now this is how the last four entries vary in AwfcTable as you sweep over rDet.

Going faster: Following the visual observation and the conclusion we can try to decrease the integration interval considerably! We first check for the numerator and denominator integral for convergence.

NIntegrate[γ[r, v, rDet]*ρ[r, v, v0]*
Exp[I ΔΦ[v]], {r, 0, Evaluate@#}, {v, 0,Evaluate@#}, AccuracyGoal -> 20,
PrecisionGoal -> 30] & /@ {.1, .25, 1, 10, 100, Infinity}


{1.98046*10^-10 + 1.57094*10^-12 I, 1.98046*10^-10 + 1.57094*10^-12 I, 1.98046*10^-10 + 1.57094*10^-12 I, 1.98046*10^-10 + 1.57094*10^-12 I, 1.98046*10^-10 + 1.57094*10^-12 I, 1.98046*10^-10 + 1.57094*10^-12 I}

The numerator converges to even if we replace Infinity with $0.25$! Now do the same for the denominator integral.

Re@NIntegrate[γ[r, v, rDet]*ρ[r, v, v0], {r, 0,Evaluate@#},
{v, 0, Evaluate@#}, AccuracyGoal -> 20,
PrecisionGoal -> 30] & /@ {.1, .25, 1, 10, 100, Infinity}


{1.98058*10^-10, 1.98058*10^-10, 1.60127*10^-10, 1.60127*10^-10, 1.60127*10^-10, 1.60127*10^-10}

The numerator converges even if we assume an integration limit of $1$ in both dimensions in place of Infinity. Following this we can have a new Awfc that will be much faster than your infinite integrals. Re is used in the denominator to get rid of the the imaginary part that arises due to numerical error.

Clear[Awfc];
Awfc[rDet_?NumericQ, v0_?NumericQ] :=
Quiet@(NIntegrate[γ[r, v, rDet]*ρ[r, v, v0]*
Exp[I ΔΦ[v]], {r, 0, .25}, {v,0, .25},
AccuracyGoal -> 12,PrecisionGoal -> 30]/
(Re@NIntegrate[γ[r, v, rDet]*ρ[r, v, v0], {r, 0,1}, {v, 0, 1},
AccuracyGoal -> 12, PrecisionGoal -> 30]));


Testing how fast it is.

AwfcTableLess =
ParallelTable[{rDet, Awfc[rDet, 0.5 v0], Awfc[rDet, v0],
Awfc[rDet, 2 v0], Awfc[rDet, 3.7 v0]}, {rDet, 0.0005, 0.010,
0.0005}]; // AbsoluteTiming


{3.879663, Null}

Compare the results with above plot.

Hope this helps you!

-
Your answer does indeed speed things up. But speed isn't really the problem here: The settings for AccuracyGoal and PrecisionGoal do alter the Arg[] of the values in AwfcTable and this value is also really important. Also, if you look at TableForm[Abs[AwfcInfTable]] you'll see that these values are greater than 1 but the weighted function Exp[I ΔΦ[v] is smaller or equal to 1. – frankundfrei Nov 26 '13 at 19:47
@frankundfrei good to hear few feed backs from you! Will give a deeper look at the issue here. Give me some time ;) – PlatoManiac Nov 26 '13 at 19:54
I think I found the mistake I've made. I'm looking into it a bit further and will post an answer if it turns out it is actually working. Thanks anyway so far, I did learn from your answer! – frankundfrei Nov 26 '13 at 21:43
By the way I noticed just now that Exp[I \[CapitalDelta]\[CapitalPhi][v]] is a function of single variable and probably not absolutely integrable with respect to your 2D weight function γ[r, v, rDet]*ρ[r, v, v0]. – PlatoManiac Nov 26 '13 at 21:59

The problems were apparently caused by how the parameters were chosen.

In the problem at hand they represent physical quantities: parameters starting with r have dimension meter, those with t have dimension second and those starting with v have meter/second.

A much more resonable choice would be to use milimeters and miliseconds, i.e. multiplying all variables starting with r or t by a factor $1000$:

tDet = 700;
t = 230;
rCloud = 2.5;
c = λ/20;
v0 = 0.00588;
rDet = 8.0/2;
rBeam = 15;


This way in AwfcTable instead of the intervall one can use a different intervall, thus both avoiding the NIntegrate::slwcon warnings and yielding reasonable results for Abs[AwfcTable]$\leq 1$:

AwfcTable = ParallelTable[
{rDet, Awfc[rDet, 0.5 v0], Awfc[rDet, v0], Awfc[rDet, 2 v0],  Awfc[rDet, 3.7 v0]},
{rDet, 0.5, 10, 0.5}];

-
Changing the units results into what I meant by scaling in my answer! Physics and math both prefers to be in agreement most of the time. – PlatoManiac Nov 27 '13 at 3:55