Errors appear after integration of an interpolating function

Question

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!

Answer 1

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  *)