# Slow integration for simple functions involving UnitBox

I would like to integrate periodic functions involving step functions. For example,

f[x_] = UnitBox[1 - Mod[x, 3]]
int[x_] := Integrate[f[u], {u, -x, x}]
intn[x_] := NIntegrate[f[u], {u, -x, x}]


Computing int[20] or intn[20] takes 0.06 seconds, plotting Plot[int[x], {x, 2, 10}] or Plot[intn[x],{x,2,10}] about 6 seconds. Changing UnitBox to HeavisidePi made things much worse, probably because the latter is not a monovalued function.

When combining such functions, it quicky takes too much time to be usable. However, from basic geometry, the integral can be computed quite easily. So my question is whether it is possible to reduce the computation time. Possibilities could i) deal with the singular points, ii) take advantage of the piecewise-constant property of f, iii) use periodicity. For iii), one workaround would be to modify the boundaries of the integral to an interval of span at most the period, and add the number of periods times the integral over one period, but I am hoping for a more straightfoward solution, if possible.

Another slightly less trivial example:

f[x_?NumericQ] =
UnitBox[0.5 - 1.25 (3. - Mod[x, 3.])] +
UnitBox[0.5 - 1.25 Mod[x, 3.]];
u[x_, t_] := NIntegrate[psitilde[xi], {xi, x - t, x + t}]
Plot[{u[x, 1.3]}, {x, 0, 10}, Exclusions -> None] // AbsoluteTiming


The Plot takes about 20s, plus it throws convergence errors (that's understandable) and the graph is hence "noisy". +6

Edit As pointed out by yarchik,

i[x_] = Integrate[f[u], u] /. {u -> x}
Plot[{i[x] - i[-x]}, {x, 2, 10}, PlotRange -> Full]


takes 0.02s and

Plot[{(Integrate[f[u], u] /. {u -> x}) - (Integrate[f[u],
u] /. {u -> -x})}, {x, 2, 10}]


takes 1.3s.

However, both curves are discontinuous.

• Evaluating int[20] takes far longer (1.4s) on my machine (macOS Mathematica11) while intn[20] gives a warning: 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. Sep 4, 2016 at 17:06
• @Sascha Yes, because I had used HeavisidePi instead of UnitBox. Edited. Sep 4, 2016 at 17:27
• Plot[Integrate[f[u], u] /. {u -> x}, {x, 2, 10}, PlotRange -> All] // Timing takes only half a second... Sep 4, 2016 at 17:43
• @yarchik Interesting, added a few words in the question about that. Sep 4, 2016 at 18:40

One idea: Add singularities to integration iterator of NIntegrate:

i2[x_?NumericQ] := NIntegrate[
f[u], {u, -x, Sequence @@ Range[1/2 + Ceiling[-x - 1/2], x], x},
Method -> {"GaussKronrodRule", "Points" -> 3}, "SymbolicProcessing" -> 0},
AccuracyGoal -> If[Abs[x] < 1, 15, Infinity]];


The Method option settings assume f[x] is piecewise linear (or a low degree polynomial). Remove "Points" and possibly "SymbolicProcessing" if needed.

Plot[i2[x], {x, 2, 10}] // AbsoluteTiming


Second idea: It's like @yarchik's idea, but using DSolve. You have to know what range of inputs you would like to have before running DSolve:

if = y /. First@DSolve[{y'[x] == f[x], y[0] == 0}, y, {x, -10, 10}]; // AbsoluteTiming
(*  {0.161635, Null}  *)

Plot[if[x] - if[-x], {x, 0, 5}] // AbsoluteTiming


DSolve returns a Piecewise function, which leads Plot to poke holes in the graph. Use Exclusions -> None to remove them. (It plots faster, too.)

• I get errors when I use i2[x] whether I use "Points" -> 3 and "SymbolicProcessing" -> 0, although it generates a plot that appears correct. The DSolve solution works great! Sep 5, 2016 at 1:03
• @JackLaVigne What version are you using? I get none in V10.4.1, on rechecking with a fresh kernel. (Thanks for letting me know, in any case.) Sep 5, 2016 at 1:05
• I was using version 11.0. I re-ran using version 10.4.1.0 and got no errors. There is progress for you :) Sep 5, 2016 at 1:08
• @JackLaVigne I added ?NumericQ` protection. Maybe that'll help. Normally I would've done that, but I was following the OP's lead. Sep 5, 2016 at 1:14
• Yup, that did the trick. Sep 5, 2016 at 1:18