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

$R_n(t) = \int_0^t 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 would like to be able to do it automatically for any n.

link|improve this question

56% accept rate
3  
As this is a convolution, might I ask why you aren't trying to use the convolution thereom? – rcollyer Feb 20 at 18:22
feedback

3 Answers

up vote 7 down vote accepted

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)
*)
link|improve this answer
feedback

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 does not remember for a generic, as opposed to a specific, k value.

link|improve this answer
feedback

You can write the definition down like that in Mathematica:

r[0, t_] := Exp[-k t] Cos[\[Omega]0 t]
r[n_, t_] := r[n, t] = Integrate[
    r[0, t - tt] r[n - 1, tt],
    {tt, 0, t}, 
    Assumptions -> k > 0 && \[Omega]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$.

link|improve this answer
Try Expand before Integrate for a moderate speedup. – Szabolcs Feb 20 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 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 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 at 18:48
feedback

Your Answer

 
or
required, but never shown

Not the answer you're looking for? Browse other questions tagged or ask your own question.