# recursive integration

I am trying to do multiple integrations recursively. For instance, I would like to do the following equation for arbitrary integer $n$:

$\displaystyle R_n(t) = \int_0^t \mathrm dt' R_0(t-t') R_{n-1}(t')$

where $R_0(t) = e^{-k t}\cos\omega_0 t$. I can do it by brute force, but I would like to be able to do it automatically for any $n$.

I think this should work:

ClearAll[r];
r[0, t_] := Exp[-k*t]*Cos[t];
r[n_, t_] := Integrate[r[0, t - td]*r[n - 1, td], {td, 0, t}]


eg

r[2,t]

(*
(\[ExponentialE]^(-k t) (2 \[ExponentialE]^(k t) k^2 -
k (2 k + t + k^2 t) Cos[t] + (k - k^3 + t + k^2 t) Sin[t]))/(2 (1 +
k^2)^2)
*)

• An interesting way to speed this up considerably is to avoid the definite integrals by doing this: r[n_, t_] := Simplify[Subtract @@ Table[Integrate[r[0, t - td]*r[n - 1, td], td] /. td -> i, {i, {t, 0}}]. It seems that by doing the indefinite integral first, I can circumvent additional convergence checks or simplifications that slow things down. But I'm surprised by how much faster it gets, e.g., for r[4,t]. – Jens May 28 '13 at 15:28
• @Jens thanks; strange that it speeds it up so much. Do you mind if I add this piece of information to the answer? – acl May 28 '13 at 22:03
• Of course you can add it to the answer! – Jens May 28 '13 at 22:13

As already mentioned, this is a convolution. Luckily, there's a more natural function to use for this problem than Integrate[], and that function is called, appropriately enough, Convolve[]. Now, since Convolve[] assumes an infinite integration region, we need a UnitStep[] multiplier in both the functions being convolved to limit the integration region to a finite interval. Here is one such implementation, making use of Convolve[], as well as a caching technique described here:

r[0, k_, t_] := Exp[-k t] Cos[t];
r[n_Integer, k_, t_] :=
Module[{kl, tl, y},
Set @@ Hold[r[n, kl_, tl_],
Simplify[
Convolve[UnitStep[y] r[n - 1, kl, y], UnitStep[y] r[0, kl, y],
y, tl], tl >= 0]];
r[n, k, t]];


Note that I had already taken the liberty to add k as an additional parameter. You can do the same thing for the $\omega_0$ factor within the cosine. The advantage of using Leonid's version of caching is that effort expended for computing, say, r[10, 4, t] is still usable for computing r[7, 8, t], since the caching remembers for a generic, as opposed to a specific, k value.

You can write the definition down like that in Mathematica:

r[0, t_] := Exp[-k t] Cos[ω0 t]
r[n_, t_] := r[n, t] =
Integrate[r[0, t - tt] r[n - 1, tt], {tt, 0, t},
Assumptions -> k > 0 && ω0 > 0]


... or is that what you considered to be brute force, and you're looking for a general solution of the recursion relation?

Anyway, the code above is awfully slow, even though I used the r := r = trick. It gets much faster when setting $k$ and $\omega_0$ to $1$.

• Try Expand before Integrate for a moderate speedup. – Szabolcs Feb 20 '12 at 18:32
• what is the r:=r= trick doing? why is it speeding things up? PS how do I 'grey' text when commenting? – BeauGeste Feb 20 '12 at 18:34
• @BeauGeste the r[t_] := r[t] = trick is called memoization. Essentially, it saves the result of r[t] at the value of t for later look-up. (This exploits the difference between OwnValues and DownValues in the evaluation loop.) – rcollyer Feb 20 '12 at 18:47
• @BeauGeste the greyed out text is just using the shorthand code markup, i.e. wrap the code in ` to get it greyed out. – rcollyer Feb 20 '12 at 18:48