# Difficulty with numerically integrating highly-oscillatory function

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?

• Your code doesn't run, you left off THz Jun 14, 2019 at 1:07
• Thanks, I added that in. Jun 14, 2019 at 19:02
• Just curious, how accurate of an answer do you need? You have two functions, the one with the Sine, and the PDF, which vary fast and slow respectively. The Sine function can be integrated in closed form over an interval 0 to $\pi$, and then you could take that value and use it as a constant for the PDF integration. I don't expect to see the shift have much of an impact. Jun 14, 2019 at 19:13
• Thanks for the reply. I'm interested in summing across the the entire slow pdf, so I think that I will need to add up all of the "fast" contributions. I think if I only needed to look over a small range (such that the slow contributions are approximately constant) than I could use your trick. Jun 14, 2019 at 19:17
• If you took the area under the fast contribution from valley to valley (via integral), and then multiplied it by value of the slow contribution at each of the fast peaks (summation) you'd be pretty darn good, is my guess. The fast contribution is symmetric, and the slow is approximately linear for the range of the period of the fast function. Worth a try. Jun 14, 2019 at 19:27

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.

• Thanks. The way you set up this function is interesting. Unfortunately though, since the plot that I have to generate has to integrate over THz (and not GHz) the range is 1000 times larger, I will have to break it up into 100,000 chunks. (And then do this for each point in the plot). On my computer it seems like it takes 15 minutes to run per point. I haven't tried the plot yet, but this seems like it will take a long time but will work! Jun 15, 2019 at 17:22
• Ahh, I’ll have to fiddle with it a little more. It’s an interesting problem! Jun 15, 2019 at 20:30
• Give the update a look, I don't think you should have much of a plot as you vary ν. Jun 16, 2019 at 17:30