Difficulty with numerically integrating highly-oscillatory function

Question

I'm interested in numerically solving a "highly-oscillatory" function. I've tried increasing the "max-recursion" and "precisiongoal", but it appears as though my answer is still not converging to the correct behavior.

nm = 10^-9;
um = 10^-6;
mm = 10^-3;
cm = 10^-2;
GHz = 10^9;
THz = 10^12;
MHz = 10^6;
cavityLength = 62 mm;
T = Tmax/(
   1 + (2 F/π)^2 Sin[2 n 2 π/λ L/2]^2) /. { n -> 2, 
    F -> 300, Tmax -> 1};
Tν = T /. {λ -> 10^8/ν, L -> cavityLength}

StartingFreq = 10^8/λ /. λ -> 795 nm;
FreqSpan = 10 THz;
spdcRANGE = 4 THz;

Plot[(Tν /. ν -> 10 MHz + ν)*
  PDF[NormalDistribution[StartingFreq, spdcRANGE], ν], {ν, 
  StartingFreq - 10 GHz, StartingFreq + 10 GHz }, PlotRange -> All, 
 Mesh -> All, PlotPoints -> 10000]
Plot[Tν*
  PDF[NormalDistribution[StartingFreq, spdcRANGE], ν], {ν, 
  StartingFreq - 10 THz, StartingFreq + 10 THz }, PlotRange -> All, 
 Mesh -> All, PlotPoints -> 100000]
Plot[NIntegrate[(Tν /. ν -> SHIFT + ν)*
   PDF[NormalDistribution[StartingFreq, spdcRANGE], ν], {ν, 
   StartingFreq - 10 THz, StartingFreq + 10 THz }, 
  Method -> "GlobalAdaptive", MaxRecursion -> 80, 
  PrecisionGoal -> 24], {SHIFT, 0, 500 MHz}, PlotPoints -> 200 , 
 Mesh -> All, MaxRecursion -> 0] 

Any ideas on what can be done? Should I just keep increasing these numbers until it works?

Answer 1

The NIntegrate can do smaller chunks, so you just need to break up your integral into smaller pieces. Define a function to integrate an interval, note it takes a list as input. We need this form below:

subInt[{l_, r_}] := NIntegrate[Tν PDF[NormalDistribution[StartingFreq, spdcRANGE], ν], {ν, l, r}]

Now divide your larger integration range into smaller chunks. Break it up into 100 pieces. Then use MovingMap to integrate each of the smaller intervals, and sum it all up.

spread = Subdivide[StartingFreq - 10 GHz, StartingFreq + 10 GHz, 100];
MovingMap[subInt, spread, 1] // Total 
(*   0.0000105274 *) 

You can subdivide it into a minimum of 80 pieces and still have it work. Get the same answer no matter the number of subdivisions (I tried up to 1000 intervals), as long as you have enough so the NIntegrate doesn't burp.

EDIT

The Tν basically forms a pulse train where the period of Tνis period = (12500 MHz)/31. The integral over a single period of Tν is

intPulse = Integrate[Tν, {ν, 0, period}]

(*  (12500000000 π)/(31 Sqrt[360000 + π^2]) *)

The number of pulse over the integration range is

pulses = 20 THz/period = 49600

So a reasonable approximation to your integral is smear out the pulse train into a constant multiplier

α = intPulse/period

(*    π/Sqrt[360000 + π^2]    *)

An approximation to your integral is

Integrate[α PDF[NormalDistribution[StartingFreq, spdcRANGE], ν], 
          {ν, StartingFreq - 10 THz, StartingFreq + 10 THz}]

(π Erf[5/(2 Sqrt[2])])/Sqrt[360000 + π^2]

In your plot, as you shift ν over the range {0, 500 MHz}, because Tν is periodic and essentially a pulse, the only real change in value should be from whether you include a pulse peak or not, so you have 49599, 49600, or 49601 pulses in your integration window. This is occurring at the tails of the distribution to boot, reducing the impact oven further.