# NIntegrate extremely piecewised functions

I often need to integrate extremely piecewised functions, like the following one (not extreme, but gives an idea):

q[time_?NumericQ] :=
Module[{q1 = 0.7, d1 = 600, p1 = 2700, ph1 = 0, q2 = 0.37, d2 = 870.,
p2 = 5000, ph2 = 0},
If[Mod[time, p1] > ph1 && Mod[time, p1] < ph1 + d1, q1, 0] +
If[Mod[time, p2] > ph2 && Mod[time, p2] < ph2 + d2, q2, 0]]


I know that I can list all the mandatory points in the NIntegrate interval List. But most of the time I don't know them.

Running NIntegrategenerally returns singularity, accuracy, etc, type of errors.

NIntegrate[q[t], {t, 0, 7*24*3600}]


I can play around with WorkingPrecision, MaximumRecursion, etc, but it is not that easy to arrive to something that returns no errors, or that runs in a reasonable time.

I generally end up running something like:

NIntegrateBis[funct_, list_] :=
Module[{step = 0.1, time = list[[1]], int = 0},
While[time < list[[2]], int = int + funct[time]*step;
time = time + step]; int]

NIntegrateBis[q, {0, 10000}]


Simple and with no messages... but no error control, no variable step, etc (it can even overpass the boundaries, although it is easy to correct that part...)

What is the best strategy to integrate my functions? Use NIntegrate with a specific set of options (that I don't know of)? Built a personal NIntegrate?

PS - I'm running optimization simulations with those functions, and so, I need it to run pretty fast.

-
Are all functions have the same form as in this example? – Andrew Jul 20 '13 at 15:43
@Andrew No. Some of them are based on order 0 or 1 Interpolation. And they are all joined up at the end. That's why I would prefer a solution on the integration end, not looking at the function itself. – P. Fonseca Jul 20 '13 at 17:05
If they are piecewise constant or piecewise linear, then they can be computed exactly (well, with MachinePrecision) if the points can be found where the pieces change (e.g. use the interpolation grid). Is that possible for your functions? – Michael E2 Jul 20 '13 at 17:43
@MichaelE2 I have done that already on everything that was possible. Please see this example as something very close to what I have. But I have a huge variation of them. I need to find a numerical solution that can master these abrupt variations. It can be very simple as the one on my example. But I would have imagined that NIntegrate could master it with a specific set of options. The simplest numerical method of integration that we first learn is the one I've recreated with my basic functions. Can't NIntegrate be degraded to that level? If not, how can I make a real fast basic one? – P. Fonseca Jul 20 '13 at 19:32
Interesting question. Unfortunately I don't have an answer for you, but thought I should offer some moral support anyway. – Oleksandr R. Jul 26 '13 at 15:57

Assuming the discontinuities in your function cannot be studied symbolically (which is the case in your example due to the inclusion of _?NumericQ), you will indeed need to roll something by hand to solve this efficiently. The built-in strategies for NIntegrate are suitable for piecewise smooth functions where the piecewise behavior can be elicited from the symbolic structure of the integrand.

Essentially you need to create a FindAllDiscontinuities function that attempts to find all function discontinuities numerically. This is only possible if you can make some assumptions about the scale of the discontinuities, which perhaps you can in your real scenario.

A rough outline of the function would be to first evenly, densely* sample the function into small sub-regions, then recursively bisect these until you are convinced* they look smooth: As you bisect, if the apparent function slope increases instead of staying roughly constant, you probably have a discontinuity.

*The asterisks are parameters that depend on the expected behavior of your real function.

-
That's it. I'm able to say that there's no peack that last shorter than 10 seconds (or that they are very rare, and it is not a problem to miss them), and that I can live with a precision of 0,1 seconds. Implement what you have described with do loops is not hard. Do you have an idea for a functional implementation, or a practical compiled implementation (where the q function is not included in the compiled function, or if it is included, it is a result of some automatic meta programming with memorisation. – P. Fonseca Jul 27 '13 at 7:34

Neither fast nor slow, but it gives no errors. Originally I used PiecewiseExpand to convert the function before passing it to NIntegrate. But it turns out NIntegrate will do it for me, if I remove the ?NumericQ from the definition of q and raise the $MaxPiecewiseCases limit to let NIntegrate do some symbolic processing. q[time_] := Module[{q1 = 0.7, d1 = 600, p1 = 2700, ph1 = 0, q2 = 0.37, d2 = 870., p2 = 5000, ph2 = 0}, If[Mod[time, p1] > ph1 && Mod[time, p1] < ph1 + d1, q1, 0] + If[Mod[time, p2] > ph2 && Mod[time, p2] < ph2 + d2, q2, 0]]; Block[{$MaxPiecewiseCases = 500},
NIntegrate[q[t], {t, 0, 7*24*3600}]
] // AbsoluteTiming
(*
{4.384593, 133030.}
*)

-

One can try to average the function say, on $[0,50000]$ with step 100 (and compile to get better performance): $$f(t)=\sum_{k=0}^{5\cdot10^2-1} q(t+100k)$$

f = Compile[{{time, _Real}},
Sum[
Module[{q1 = 0.7, d1 = 600, p1 = 2700, ph1 = 0, q2 = 0.37,
d2 = 870., p2 = 5000, ph2 = 0},
If[Mod[tt, p1] > ph1 && Mod[tt, p1] < ph1 + d1, q1, 0] +
If[Mod[tt, p2] > ph2 && Mod[tt, p2] < ph2 + d2, q2, 0]], {tt,
time, time + 5 10^4 - 100, 100}]
];


hoping the result will be "less piecewise".

Plot[f[x], {x, 0, 100}]


Trying to integrate f gets the result, still producing error messages though. But since the integrand is not smooth, the precision to be had directly depends on number of points at which the function is evaluated. One can play here with methods/parameters/precision, of course, to obtain better results.

-

In my case, just using Compile speeds up everything pretty much. The speed up gain is more than 100 times.

q2 = Compile[{{time, _Real}}, With[{q1 = 0.7, d1 = 600, p1 = 2700, ph1 = 0, q2 = 0.37, d2 = 870., p2 = 5000, ph2 = 0}, If[Mod[time, p1] > ph1 && Mod[time, p1] < ph1 + d1, q1, 0] + If[Mod[time, p2] > ph2 && Mod[time, p2] < ph2 + d2, q2, 0]], CompilationTarget -> "C"];

NIntegrateBis2 = Compile[{{list1, _Real}, {list2, _Real}}, Block[{step = 0.1, time = list1, int = 0.0}, While[time < list2, int = int + q2[time]*step; time = time + step]; int], CompilationTarget -> "C"];

In[14]:= NIntegrateBis2[0, 10000] // Timing
Out[14]= {0.028000, 2323.69}

-
Thank you for the analysis. But that obliges me to write an integration function for each different q function. Is there a way of passing the function into the compiled integration? Or to meta-program the compilation? – P. Fonseca Jul 20 '13 at 12:21
@P.Fonseca I'm not sure, I don't think you can give a function name as an input in the Compile command. Anyway, although probably this is not the most elegant solution to your problem(it might be the fastest though); it is generally good practice to consider using Compile . – Ali Jul 20 '13 at 12:33