7
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – Sascha Sep 4 '16 at 17:06
  • $\begingroup$ @Sascha Yes, because I had used HeavisidePi instead of UnitBox. Edited. $\endgroup$ – anderstood Sep 4 '16 at 17:27
  • $\begingroup$ Plot[Integrate[f[u], u] /. {u -> x}, {x, 2, 10}, PlotRange -> All] // Timing takes only half a second... $\endgroup$ – yarchik Sep 4 '16 at 17:43
  • $\begingroup$ @yarchik Interesting, added a few words in the question about that. $\endgroup$ – anderstood Sep 4 '16 at 18:40
5
$\begingroup$

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 -> {"GlobalAdaptive", 
    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

Mathematica graphics

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

Mathematica graphics

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

$\endgroup$
  • $\begingroup$ 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! $\endgroup$ – Jack LaVigne Sep 5 '16 at 1:03
  • $\begingroup$ @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.) $\endgroup$ – Michael E2 Sep 5 '16 at 1:05
  • $\begingroup$ I was using version 11.0. I re-ran using version 10.4.1.0 and got no errors. There is progress for you :) $\endgroup$ – Jack LaVigne Sep 5 '16 at 1:08
  • $\begingroup$ @JackLaVigne I added ?NumericQ protection. Maybe that'll help. Normally I would've done that, but I was following the OP's lead. $\endgroup$ – Michael E2 Sep 5 '16 at 1:14
  • $\begingroup$ Yup, that did the trick. $\endgroup$ – Jack LaVigne Sep 5 '16 at 1:18

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.