Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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

share|improve this question
3  
As this is a convolution, might I ask why you aren't trying to use the convolution theorem? –  rcollyer Feb 20 '12 at 18:22
    
Another question related to recursive integrals is: Solving a Volterra integral equation numerically –  Jens Jun 15 '12 at 22:10
add comment

3 Answers

up vote 10 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)
*)
share|improve this answer
    
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
add comment

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

share|improve this answer
    
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
1  
@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
add comment

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.

share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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