Memoization (or otherwise) for a recurrent definition of a (double) sequence of functions

Question

Beside I've found answers to similar questions (e.g. this), I can't adapt them to my case.

I am interested in the following sequence of functions $p_{m,n}(t)$, where $m$ and $n$ are integer:

$$ \begin{aligned}p_{m,n}(t)&:=\begin{cases} 0, & \text{if $m<1$ or $n<1$ or $n>m+1$,}\\[4mm]\displaystyle \int_0^t e^{b(t-s)n} q_{m,n}(s) \,ds, & \text{otherwise, where}\end{cases}\\[3mm] q_{m,n}(s)&:= \chi_{m=1}\chi_{n=2}\, b + b n p_{m,n+1}(s)\\[2mm]&\qquad +\chi_{m\geq 2} b n(n-1)\bigl(p_{m-1,n-2}(s)+p_{m-1,n-1}(s)+p_{m-1,n}(s)\bigr)\end{aligned} $$

Here $b>0$ is a parameter and $\chi$ is equal to $1$ if the subscript happens and to $0$ otherwise.

As you may see, the definition is indeed recurrent, since $p_{m,m+1}(\cdot)$ does not really contain the term with $p_{m,m+2}(\cdot)\equiv 0$, and hence can be found recurrently, after this one can find $p_{m,m}(\cdot)$ which depends on $p_{m,m+1}(\cdot)$, and so on down to $p_{m,1}(\cdot)$.

I coded this as follows:

ClearAll[p, b];
$Assumptions = {b > 0};
p[m_, n_] := p[n, m] = If[ m < 1 || n < 1 || n > m + 1, 0 # &,
            Module[{s}, Integrate[Exp[(# - s) b n]
                (KroneckerDelta[m, 1] KroneckerDelta[n, 2] b  
                 + b n p[m, n + 1][s] + If[m < 2, 0,
                   b n (n - 1) (p[m - 1, n - 2][s] + p[m - 1, n - 1][s] 
                                +p[m - 1, n][s])]),
                {s, 0, #}] &]];

It works, but very slowly:

In[5]:= p[5, 6][t] // AbsoluteTiming

Out[5]= {83.3974, 
 15 + 86 E^(b t) - 47350 E^(5 b t) + 5034 E^(6 b t) + 
  3240 E^(3 b t) (34 + 15 b t) + 10 E^(2 b t) (361 + 30 b t) + 
  5 E^(4 b t) (-14311 + 21216 b t)}

How can I speed up this?

(Remark: for me the integration variable $s$ looks weird in the module description, but otherwise, I got a mistake, probably because of names conflict at some stage.)

Answer 1

Memoize p properly (you had a typo), and make sure you use Evaluate inside of function definitions:

Clear[p];
q[{m_, n_}, s_] := 
  KroneckerDelta[m, 1] KroneckerDelta[n, 2] b + b n p[m, n + 1][s] +
  If[m < 2, 0, b n (n - 1) (p[m - 1, n - 2][s] + p[m - 1, n - 1][s] + p[m - 1, n][s])]
p[m_, n_] /; m < 1 || n < 1 || n > m + 1 = Function[t, 0];
p[m_, n_] := p[m, n] =
  Function[t, Evaluate[Integrate[Exp[(t - s) b n] q[{m, n}, s], {s, 0, t}]]]

Test:

p[5, 6][t] // AbsoluteTiming

{12.9964,
 -(2/5) E^(b t) (2743 + 118125 E^(4 b t) - 12618 E^(5 b t) + 540 b t + 
 900 E^(2 b t) (-377 + 42 b t) - 100 E^(3 b t) (-758 + 2301 b t) + 
 450 E^(b t) (345 + 174 b t + 20 b^2 t^2))}

---

If you want much much higher speed, you need to use a custom integration subroutine. Here I've tried to assemble one that does the integral you need, but not much more:

myint[n_, X_] := Function[t, Evaluate[Expand[
  Expand[F[0, 0]*TrigToExp[X[s]]] /. {
    c_. s^α_. E^(β_. b s) -> c*F[α, β]/F[0, 0], 
    c_. E^(β_. b s) -> c*F[0, β]/F[0, 0], 
    c_. s^α_. -> c*F[α, 0]/F[0, 0]} /. 
  F[α_, β_] :> If[β == n, (E^(b n t) t^(α + 1))/(α + 1), 
    E^(b n t)/(b (n - β))^(α + 1)*(Gamma[α + 1] - 
      FunctionExpand[Gamma[α + 1, b t (n - β)]])]]]]

In terms of this we can get superfast recursion:

Clear[p];
q[{m_, n_}, s_] := 
  KroneckerDelta[m, 1] KroneckerDelta[n, 2] b + b n p[m, n + 1][s] + 
  If[m < 2, 0, b n (n - 1) (p[m - 1, n - 2][s] + p[m - 1, n - 1][s] + p[m - 1, n][s])]
p[m_, n_] /; m < 1 || n < 1 || n > m + 1 = Function[t, 0];
p[m_, n_] := p[m, n] = myint[n, q[{m, n}, #] &]

Test: roughly 100 times faster than before!

p[5, 6][t] // AbsoluteTiming

{0.093304,
 -((5486 E^(b t))/5) - 62100 E^(2 b t) + 
 135720 E^(3 b t) - 30320 E^(4 b t) - 47250 E^(5 b t) + 
 25236/5 E^(6 b t) - 216 b E^(b t) t - 31320 b E^(2 b t) t - 
 15120 b E^(3 b t) t + 92040 b E^(4 b t) t - 3600 b^2 E^(2 b t) t^2}