# NIntegrate not working after interpolation

I used linear interpolation between points:

T = 1;
w = 0.05;
num = 1000;
A = 1;

pulse[x_] := A*(UnitStep[x + w*T/2] - UnitStep[x - w*T/2])

fun =
Table[pulse[x] + 0.2*(RandomReal[] - 0.5), {x, -T/2, T/2,
T/(num - 1)}];
funX = Table[i, {i, -T/2, T/2, T/(num - 1)}];

funINT =
Interpolation[Transpose[{funX[[All]], fun[[All]]}],
InterpolationOrder -> 1];

ListPlot[Transpose[{funX[[All]], fun[[All]]}], PlotRange -> All,
Filling -> Axis, Frame -> True,
FrameLabel -> {"Time [s]", "Amplitude [V]"},
PlotLegends -> {"Pulse"},
ImageSize ->
Large]


This produces exactly what I want:

But in order to calculate the coefficients of complex Fourier series, I have to calculate the following integral:

    cn[k_] = NIntegrate[funINT[t]*Exp[-I*2*Pi*k*t/T]/T, {t, -T/2, T/2},
Method -> "Trapezoidal"];


which is not working. I get an error saying that the integrand has evaluated to non-numerical values. Any ideas on what the problem is?

EDIT: If there is a method to calculate the integral between the discrete points, that might be even better in my case. However, I couldn't find one.

• Try cn[k_] := , i.e. use SetDelayed. As it stands, you are telling Mathematica to evaluate the RHS immediately, but NIntegrate cannot do that for non-numerical k. – Marius Ladegård Meyer Dec 3 '16 at 22:09
• @MariusLadegårdMeyer: I am familiar with that solution, yes. However it has a huge problem: Once I try to evaluate the Fourier series for 1000 terms I takes for ages to evaluate rekon[t, 1000] if rekon[t_, st_] := Sum[cn[n]*Exp[I*2*Pi*n*t/T], {n, -st, st}]. That is why I would like to somehow do the integration with k as parameter. – skrat Dec 3 '16 at 22:16

You could find a general symbolic Fourier coefficient for a linear polynomial and use the formula to integrate the interpolating function piecewise. If you're content with machine precision (double precision), then you can Compile it for really great speed.

(* Basic integral formulas *)
ClearAll[cn0];
cn0[0][{t0_, t1_}, {x0_, x1_}] = (* k == 0 is a special case *)
Integrate[(x0 + (x1 - x0)/(t1 - t0) (t - t0)) * Exp[-I*2*Pi*k*t/T]/T /. k -> 0,
{t, t0, t1}];
cn0[k_][{t0_, t1_}, {x0_, x1_}] =
Integrate[(x0 + (x1 - x0)/(t1 - t0) (t - t0)) * Exp[-I*2*Pi*k*t/T]/T,
{t, t0, t1}];

(* Coefficient function *)
Clear[cn];
cn[0] = Total@
MapThread[ (* map over interpolation segments *)
cn0[0],
{Partition[funX, 2, 1], Partition[fun, 2, 1]}];
cn[k_] = Total@
cn0[k],
{Partition[funX, 2, 1], Partition[fun, 2, 1]}];

(* Compiled version *)
cnC = With[ (* basic integrals *)
{i0 = Function[{t0, t1, x0, x1}, (* k == 0 is a special case *)
Evaluate@ Integrate[
(x0 + (x1 - x0)/(t1 - t0) (t - t0)) * Exp[-I*2*Pi*k*t/T]/T /. k -> 0,
{t, t0, t1}]],
i = Function[{t0, t1, x0, x1},
Evaluate@ Integrate[
(x0 + (x1 - x0)/(t1 - t0) (t - t0)) * Exp[-I*2*Pi*k*t/T]/T,
{t, t0, t1}]]},
Compile[{{k, _Integer}, {t, _Real, 1}, {x, _Real, 1}},
Total@If[k == 0,
i0[Most[t], Rest[t], Most[x], Rest[x]], (* vectorized for speed *)
i[Most[t], Rest[t], Most[x], Rest[x]]]
]];


Checks and comparison of speeds:

(* OP's method for comparison *)
cn1[k_] :=
NIntegrate[funINT[t]*Exp[-I*2*Pi*k*t/T]/T, {t, -T/2, T/2},
Method -> "Trapezoidal"];

res1 = Table[cn1[k], {k, 0, 5}] // AbsoluteTiming
res2 = Table[cn[k], {k, 0, 5}] // AbsoluteTiming
res3 = Table[cnC[k, funX, fun], {k, 0, 5}] // AbsoluteTiming
(*
{6.41549, {0.0509924, 0.0485667 + 0.000384561 I,
0.0489479 + 0.000373475 I, 0.0489737 + 0.000976852 I,
0.0459565 + 0.00132399 I, 0.0452069 + 0.000833868 I}}

{0.154046, {0.0509924, 0.0485667 + 0.000384561 I,
0.0489479 + 0.000373475 I, 0.0489737 + 0.000976852 I,
0.0459565 + 0.00132399 I, 0.0452069 + 0.000833868 I}}

{0.001207, {0.0509924 + 0. I, 0.0485667 + 0.000384561 I,
0.0489479 + 0.000373475 I, 0.0489737 + 0.000976852 I,
0.0459565 + 0.00132399 I, 0.0452069 + 0.000833868 I}}
*)

res1 - res2
res2 - res3
(*
{6.26145, {-6.245*10^-17, -1.27026*10^-10 - 3.56074*10^-10 I,
5.71595*10^-12 - 8.24057*10^-11 I, -5.08276*10^-10 - 1.50366*10^-11 I,
8.25427*10^-11 - 4.73669*10^-10 I, -3.07932*10^-10 + 1.72791*10^-10 I}}

{0.152839, {-1.38778*10^-17 + 0. I, -3.7817*10^-15 - 1.49451*10^-14 I,
1.1019*10^-14 - 3.11627*10^-15 I, -1.10328*10^-15 - 2.62073*10^-15 I,
-3.42781*10^-15 - 8.26162*10^-17 I, 1.88738*10^-15 - 1.13711*10^-15 I}}
*)


So cn is almost 50 times faster than NIntegrate and cnC is over 100 times faster than cn.

• I was wandering should I answer this question or first wait and see would you answer it... :) – Anton Antonov Dec 4 '16 at 0:12
• @AntonAntonov Ha-ha. :) – Michael E2 Dec 4 '16 at 0:13
• @skrat The (x0 + (x1 - x0)/(t1 - t0) (t - t0)) is the point-slope formula for the interpolating line between {t0, x0} and {t1, x1} as a function of t, which is the formula used by the order-1 interpolation. So i (and i0 for k == 0) is the value of the integral over a general piece of the interpolation in terms of k and two interpolation points (t0, t1, x0, x1). – Michael E2 Dec 4 '16 at 13:31

There are several ways to approach answering this question.

## Not using NIntegrate

One approach is to use direct Trpezoidal formula integration as shown in this answer of "Is it possible to compute with the trapezoidal rule by numerical integration?".

## Using NIntegrate's range specification

If there is a method to calculate the integral between the discrete points, that might be even better in my case. However, I couldn't find one.

NIntegrate's range specification can take range splitting points. In this case we put the points funX in the range specification with Evaluate[{t, Sequence @@ funX}].

Clear[cn]
cn[k_] :=
NIntegrate[funINT[t]*Exp[-I*2*Pi*k*t/T]/T,
Evaluate[{t, Sequence @@ funX}], PrecisionGoal -> 4,
AccuracyGoal -> 6, MaxRecursion -> 1,
Method -> {"GlobalAdaptive", "SymbolicProcessing" -> 0}];

Clear[rekon]
rekon[t_, st_] := Sum[cn[n]*Exp[I*2*Pi*n*t/T], {n, -st, st}]

In[74]:= AbsoluteTiming[res = rekon[t, 100];]
Out[74]= {33.411, Null}

In[75]:= AbsoluteTiming[res = rekon[t, 1000];]
Out[75]= {360.723, Null}


## Using a specially made rule for a list of functions

We can apply the approach given in this answer of the question "How to calculate the numerical integral more efficiently?".

We make a list of all the integrands formed in rekon and we integrate them at once with ArrayOfFunctionsRule. The computations are 8-10 times faster.

Import["https://raw.githubusercontent.com/antononcube/\
MathematicaForPrediction/master/Misc/ArrayOfFunctionsRule.m"]

Clear[n, st, t]
With[{st = 1000},
factors = Table[Exp[I*2*Pi*n*t/T], {n, -st, st}];
funcs = Table[funINT[t]*Exp[-I*2*Pi*k*t/T]/T, {k, -st, st}];
]

AbsoluteTiming[
res = NIntegrate[1, Evaluate[{t, Sequence @@ funX}],
PrecisionGoal -> 4, AccuracyGoal -> 6, MaxRecursion -> 20,
Method -> {"UnitCubeRescaling", "FunctionalRangesOnly" -> True,
Method -> {"GlobalAdaptive", "SingularityHandler" -> None,
Method -> {ArrayOfFunctionsRule, "Functions" -> funcs}}}];
]

(* {39.5742, Null} *)

res.factors


Michael's and Anton's answers involve the manual splitting of the piecewise linear functions involved in the computation of the Fourier coefficients at the interpolation points. Yet another way to tell Mathematica to automatically split the integrand before integrating is to use the option setting Method -> "InterpolationPointsSubdivision", like so:

With[{T = 1, w = 0.05, num = 1000, A = 1},
BlockRandom[SeedRandom["many pulses"]; (* for reproducibility *)
funINT = Interpolation[Table[{x, A (UnitStep[x + w T/2] - UnitStep[x - w T/2]) +
RandomReal[{-1, 1} 0.1]},
{x, -T/2, T/2, T/(num - 1)}], InterpolationOrder -> 1]]];

With[{T = 1},
Table[NIntegrate[funINT[t] Exp[-2 π I k t/T]/T, {t, -T/2, T/2},
Method -> "InterpolationPointsSubdivision"], {k, 0, 5}]]
{0.0511259271805279, 0.05148235941558258 - 0.001596805628001697*I,
0.05010102394935553 - 0.0007387015194085259*I,
0.05130250748609571 + 0.0018324339082365982*I,
0.04620591693520239 - 0.00029540402973357425*I,
0.04467904044269826 + 0.0016773681935466633*I}

• Yeah of course I should have mentioned "InterpolationPointsSubdivision" in my answer! – Anton Antonov Dec 10 '16 at 17:24