0
$\begingroup$

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

$\endgroup$
  • 1
    $\begingroup$ Your memoization p[m_, n_] := p[n, m] = ... has a typo. $\endgroup$ – Roman Apr 19 at 17:53
5
$\begingroup$

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}
$\endgroup$

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.