Performing iterated integrals

Question

I would like to define a function $f_n(x)$ and plot it with Manipulate. The function is defined by $$ f_{n+1}(x)=\frac{1}{\int_0^1(f_n(x'))^2dx'}\left(\int_0^xf_n(x')dx'-\int_0^1f_n(x')dx'\right) $$ The key part is $\int_0^xf_n(x')dx'$, the iterated integral of the original function $f_1(x)$ and the rest is details - the other term in parentheses subtracts the average in each iteration, and the fraction normalizes in each iteration.

Before doing this more complex mapping, the first step is to be able to do iterated integrals, which I test with sine:

n = 2;
f[1][x_] := Sin[x];
For[j = 1, j < n, j++, f[j + 1][x_] := Integrate[f[j][x], x]]
Plot[f[1][x], {x, 0, 10}]

Disappointingly, I get

SetDelayed::write: Tag Sin in Sin[1][x_] is Protected.
SetDelayed::write: Tag Sin in Sin[2][x_] is Protected.

and an empty plot.

If I delete the index j

n = 2;
f[x_] := Sin[x];
For[j = 1, j < n, j++, f[x_] := Integrate[f[x], x]]
Plot[f[x], {x, 0, 10}]

gives recursion errors as the plot is attempted:

RecursionLimit::reclim2: Recursion depth of 1024 exceeded during evaluation of f[0.000204286].

What is the proper way to do an iterated integral like this?

Ideally, I'd be able to a sequence of already-defined functions $f_n(x)$ so they could be displayed by Manipulate later without needing to be recalculated.

Answer 1

Why not just nest the integration?

integrals = NestList[Integrate[#, {x, 0, x}] &, Sin[x], 5]
Plot[integrals, {x, -1, 1}, PlotLegends -> "Expressions"]

You can then get the function $f_n(x)$ for a given $n$ with integrals[[n]]

Since you mentioned Manipulate, you can do this and control parameters like a and the number of integration levels. It would not take much more to adapt this to the particular formula in your question:

Manipulate[
 With[{fn = Nest[Integrate[#, {x, 0, x}] &, Sin[2 π a*x], n]},
  Plot[fn, {x, -4 π, 4 π}, PlotRange -> {-4, 4}]
  ], {a, 0, 3}, {n, 0, 5, 1}]

Answer 2

Here is an example of such a definition. I will use f[1] = 1 rather than Sin[x], so all results are distinct:

ClearAll[f]
f[1][x_] = 1;
f[i_Integer][x_] := Integrate[f[i - 1][x], x]

f[5][t]
(* Out: t^4/24 *) 

Plot[
  Evaluate[ f[#][x] & /@ Range[5] ], {x, -2, 2},
  PlotLegends -> "Expressions"
]

---

Here is the same idea applied to Sin[x] as a starting point:

ClearAll[g]
g[1][x_] = Sin[x];
g[i_Integer][x_] := Integrate[g[i - 1][x], x]

Plot[
  Evaluate[g[#][x] & /@ Range[4] ], {x, -2, 2},
  PlotLegends -> "Expressions"
]

Answer 3

On the off chance that an explanation would be appreciated, one that explains how to think about definitions and would help one avoid making such mistakes, instead of just some working code, I offer the following.

First, gratuitous advice: Why should I avoid the For loop in Mathematica? Using Do[body, {j, n-1}] instead of For[j = 1, j < n, j++, body] won't solve any of the problems in the OP (hence, "gratuitous"), but the advantages are discussed in the linked Q&A. When body effectively has the form y = G[y], then Nest[] is an appropriate tool to use instead of For[] or Do[] (see @flinty's answer).

1A. The first error is explained by this:

• Lingering Definitions: when calculations go bad

If I execute the OP's first-example code, I don't get any errors on a clean start. The errors are due to a lingering definition, probably f[x_] := Sin[x], but possibly f = Sin. The following probably should be the advice in the link above (but it's not):

• Start your definition of f with ClearAll[f].

Add the line ClearAll[f]; to the first example and the errors disappear.

1B. Even with the ClearAll fix, the first example contains another bug not exhibited in the example code:

f[2][x]

$RecursionLimit::reclim2: Recursion depth of 1024 exceeded during evaluation of f[j][x].

(*  Hold[\[Integral]f[j][x] \[DifferentialD]x]  *)

There are various ways to look at the issues. They concerned the order of evaluation, including held arguments, and the following:

• Understand the difference between Set and SetDelayed

• Understand what Set (=) really does

• What is the distinction between DownValues, UpValues, SubValues, and OwnValues?

You can see the effect of Set or SetDelayed with ? f, which shows the effect is to define SubValues for f. (in the first example you would need SubValues[f], which you can read about in the third link in the list.

ClearAll[f, j];
n = 2;
f[1][x_] := Sin[x];
For[j = 1, j < n, j++, f[j + 1][x_] := Integrate[f[j][x], x]];

SubValues[f]

(*
  {HoldPattern[f[1][x_]] :> Sin[x], 
   HoldPattern[f[2][x_]] :> Integrate[f[j][x], x]}
*)

Note the effect of the definition in the For loop. Because the right-hand side is "held" (not evaluated), the literal f[j] appears in replace rule in SubValues. At this point, j has a value 2 after the For loop. So f[2][x] evaluates as follows:

f[2][x] -->  apply the SubValue for f[2][x_]
Integrate[f[j][x], x] -->  Integrate now evaluates f[j][x]
  f[j][x] -->  j=2, so...
  f[2][x] -->  apply the SubValue for f[2][x_]
  Integrate[f[j][x], x], x] -->  Integrate now evaluates f[j][x]
    f[j][x] -->  j=2, so...
    f[2][x] --> ... ad infinitum or until $RecursionLimit is reached

You could clear j or set j = 3 or to another value to see how it affects the result.

i) One workaround is to use Set instead of SetDelayed in the For loop:

f[j + 1][x_] = Integrate[f[j][x], x]

With Set, the RHS is evaluated before the definition is made. Then the SubValues become the following and contain no recursive reference to f:

SubValues[f]

(*
  {HoldPattern[f[1][x_]] :> Sin[x],
   HoldPattern[f[2][x_]] :> -Cos[x]}
*)

ii) Another workaround is to inject the value of j into the definition using With:

With[{j = j}, f[j + 1][x_] := Integrate[f[j][x], x]]

Now the SubValues are the following:

SubValues[f]

(*
  {HoldPattern[f[1][x_]] :> Sin[x], 
   HoldPattern[f[2][x$_]] :> Integrate[f[1][x$], x$]}
*)

Note that the reference to f in the definition of f[2] is to f[1], so we won't get an infinite recursion.

2. The issue here is similar to 1B. In this case you can see the effect of SetDelayed with ? f or with DownValues[f]):

ClearAll[f, j];
n = 2;
f[x_] := Sin[x];
For[j = 1, j < n, j++, f[x_] := Integrate[f[x], x]]

DownValues[f]

(*
   HoldPattern[f[x_]] :> Integrate[f[x], x]}
*)

Because the right-hand side of SetDelayed is "held" (not evaluated), the expression f[x] appears in replacement rule in DownValues. Thus we get infinite recursion. The first workaround in 1B, using Set, can be applied here:

ClearAll[f, j, x];
n = 2;
f[x_] := Sin[x];
For[j = 1, j < n, j++, f[x_] = Integrate[f[x], x]]

DownValues[f]

(*  {HoldPattern[f[x_]] :> -Cos[x]}  *)

Of course, now the body has the form y = G[y], and thus Nest[] is an appropriate tool to use instead of For[].