# Errors appear after integration of an interpolating function

I have a little question about integration of an interpolating function. There is an error in assumptions of the integral.

I have a function that was obtained with an interpolation method:

try[1] = Interpolation[simmatrix[1], InterpolationOrder -> 3];


where 'simmatrix1' is a list of values. So I can plot the previos function with the following:

Plot[try[1][x], {x, 1.5, 2}]


and it is the output:

If I write

try[1]


I obtaine (I post a figure for a clear visualization):

Now let me define the following recursive funcions:

l[0] := try[1][x[0]]
l[h_] := Integrate[l[h - 1], {x[h - 1], 1.5, x[h]}]


So if I write

l[1]


I obtaine

but if I write

l[2]


I obtaine

with the error General::ivar.

If now, for example, I need to evaluate the numerical integral over 'l2', in specific boundaries

NIntegrate[l[2], {x[1], 0, 0.1}, AccuracyGoal -> 5]


I obtaine the following errors (General::ivar, NIntegrate::inumr):

Is there a way to fix or solve the errors messages ? I tried with several assumptions but it didn't work. Thanks for any tips and helps!

• Where did you define x[i]? Commented Jan 9, 2022 at 12:32
• I don't define 'x[i]', I think it is not necessary; they are mute variables. Commented Jan 9, 2022 at 13:42
• Does l[0] := try[1][x]; l[h_Integer?Positive] := Integrate[l[h - 1], x]; work for you? Note that Integrate[InterpolatingFunction[..][x], x] in effect somputes $\int_{x_0}^x f(t)\,dt$, where $x_0$ is the beginning of the domain of the interpolating function. The variable stays x in all cases instead of x[0], x[1], etc. But that's what I'm asking, do the variables matter or can they all be x? Commented Jan 9, 2022 at 15:32
• Yes, they can be all 'x'. Commented Jan 9, 2022 at 15:58

Since no data was supplied, I made some up. See What does the construct f[x_] := f[x] = ... mean? for memoization and the mem pattern. The main idea below is to use Integrate[InterpolatingFunction[<..>][x], x] to compute the integral $$\int_0^x f(t) \; dt$$, which it does by integrating symbolically the piecewise polynomial function of the interpolating function. In this way, one avoids the numerical difficulties arising from the weak singularities at the interpolation nodes. An alternative is to use Method -> "InterpolationPointsSubdivision" with NIntegrate, but I don't recommend it for the OP's use-case; see How to integrate functions of linearly interpolated data?

SeedRandom[0];
data = Transpose@{#, Cos[#]} &@RandomReal[{Pi, 2 Pi}, 100];
try[1] = Interpolation[data,
InterpolationOrder -> 3]; (* default = 3
InterpolationOrder -> 3 is unnecessary *)

l // ClearAll;
mem : l[0] := mem = try[1][x];
mem : l[h_Integer?Positive] := mem = Integrate[l[h - 1], x];


Domain:

xdom = First[l[0] /. x -> "Domain"]
(*  {3.20169, 6.24452}  *)


Inspect l[1] $$\approx \sin x, \pi\le x \le 2\pi$$:

Another check:

Plot[{l[0], l[1], l[2]},
Prepend[xdom, x],
PlotLegends -> "Expressions"]


### Alternative

Instead of being an expression depending on x, one can also define l[h] to behave like a pure function so you could write l[2][4.1] or l[1][x]. Note that Derivative[-1][f] effectively calls Integrate to antidifferentiate f.

l // ClearAll;
mem : l[0] := mem = try[1];
mem : l[h_Integer] := mem = Derivative[-h][try[1]];
(* optional formatting *)
Format[HoldPattern@l[0]] = Superscript[l, "(0)"];
Format[l[h_Integer]] = **strong text**Derivative[-h][l]; (* use h instead of -h if desired *)

Plot[{l[0][x], l[1][x], l[2][x]}, Prepend[xdom, x],
PlotLegends -> "Expressions"]


Evaluate:

l[1][4.1]

(*  -0.758219  *)