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

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

• Your memoization p[m_, n_] := p[n, m] = ... has a typo. – Roman Apr 19 '19 at 17:53

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}