6
$\begingroup$

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$.

$\endgroup$
2

3 Answers 3

10
$\begingroup$

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)
*)
$\endgroup$
3
  • $\begingroup$ 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]. $\endgroup$
    – Jens
    May 28, 2013 at 15:28
  • $\begingroup$ @Jens thanks; strange that it speeds it up so much. Do you mind if I add this piece of information to the answer? $\endgroup$
    – acl
    May 28, 2013 at 22:03
  • $\begingroup$ Of course you can add it to the answer! $\endgroup$
    – Jens
    May 28, 2013 at 22:13
10
$\begingroup$

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.

$\endgroup$
3
$\begingroup$

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$.

$\endgroup$
4
  • $\begingroup$ Try Expand before Integrate for a moderate speedup. $\endgroup$
    – Szabolcs
    Feb 20, 2012 at 18:32
  • $\begingroup$ what is the r:=r= trick doing? why is it speeding things up? PS how do I 'grey' text when commenting? $\endgroup$
    – BeauGeste
    Feb 20, 2012 at 18:34
  • 1
    $\begingroup$ @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.) $\endgroup$
    – rcollyer
    Feb 20, 2012 at 18:47
  • $\begingroup$ @BeauGeste the greyed out text is just using the shorthand code markup, i.e. wrap the code in ` to get it greyed out. $\endgroup$
    – rcollyer
    Feb 20, 2012 at 18:48

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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