What is the fastest way to calculate series from recurrence solutions with no closed form?

Question

I am not very experienced with Mathematica. I want to calculate the series from recurrence solutions. These solutions have no closed-form. I want to improve the approximation and hence need a good number of terms. I know I can use RecurrenceTable command in Mathematica, but it is not computational efficient(calculations are going on from 2 days). I know this is because of the reason that every-time it has to calculate one new term, it has to go for all the previous terms.
The function looks like:

$$\text{RecurrenceTable}\left[\left\{p(n+1)=\frac{2 \beta p(n)+2 p(n-1)-\alpha q(n)}{n+1},q(n+1)=\frac{\alpha p(n)+2 \beta q(n)+2 q(n-1)}{n},p(0)=-\frac{2 \beta }{\alpha },p(1)=\frac{\alpha ^2+4 \beta ^2}{\alpha },q(0)=-1,q(1)=0\right\},\{p(n),q(n)\},\{n,0,100\}\right] $$

I want to know what would be the best way to calculate the same.

Answer 1

As mentioned by corey979, memorization will speed up the calculations

Clear["Global`*"]

p[0] = -2 β/α; p[1] = (α^2 + 4 β^2)/α;
p[n_Integer?Positive] := p[n] =
   (2 β*p[n - 1] + 2 p[n - 2] - α*q[n - 1])/n // Simplify;

q[0] = -1; q[1] = 0;
q[n_Integer?Positive] := q[n] =
   (α*p[n - 1] + 2 β*q[n - 1] + 2 q[n - 2])/(n - 1) // Simplify;

AbsoluteTiming[p[50];]

(* {36.9799, Null} *)

Almost all of the time was spent doing the simplifications.

Length@DownValues[p]

(* 52 *)

Since all but the last of the corresponding q values were already calculated in evaluating the p values, q values go much faster.

AbsoluteTiming[q[50];]

(* {0.853231, Null} *)

Length@DownValues[q]

(* 52 *)

Looking at examples

p@10

(* (1/(1625702400 α))β (30563 α^10 - 
   2 α^8 (962581 + 480430 β^2) + 
   8 α^6 (4280439 + 3642776 β^2 + 205756 β^4) + 
   64 α^4 (-2925591 - 2408022 β^2 + 870180 β^4 + 
      283784 β^6) - 
   256 α^2 (-980385 + 250644 β^2 + 2610624 β^4 + 
      1030352 β^6 + 88976 β^8) + 
   28672 (-945 + 2130 β^2 + 8520 β^4 + 4368 β^6 + 
      688 β^8 + 32 β^10)) *)

q@10

(* (1/3657830400)(4096 α^10 - α^8 (335872 + 
     556135 β^2) + α^6 (8454144 + 26321138 β^2 + 
     5401420 β^4) + 
  8 α^4 (-9424896 - 39546855 β^2 - 9479064 β^4 + 
     1439780 β^6) - 
  64 α^2 (-3283968 - 15113973 β^2 + 5114814 β^4 + 
     7943308 β^6 + 1088920 β^8) + 
  256 (-483840 - 2446233 β^2 + 6587220 β^4 + 7425312 β^6 + 
     1797968 β^8 + 114064 β^10)) *)

Answer 2

I recommend using Do in place of RecurrenceTable. It should be faster and now memoization is automatic.

With[{m = 3},
  Clear[p, q];
  p[0] = -2 β/α;
  p[1] = (α^2 + 4 β^2)/α;
  q[0] = -1;
  q[1] = 0;
  Do[
    p[n + 1] = (2 β p[n] + 2 p[n - 1] - α q[n])/(n + 1); 
    q[n + 1] = (α p[n] + 2 β q[n] + 2 q[n - 1])/n,
    {n, 1, m}]]

Then

Definition[p]

Definition[q]

Going up yo n = 100 symbolically will take a lot memory and fair amount of time.

Answer 3

One more possibility is to recast your coupled difference equations as repeated matrix-vector multiplication:

$$\begin{pmatrix}p(n+1)\\q(n+1)\\p(n)\\q(n)\end{pmatrix}=\begin{pmatrix}\tfrac{2\beta}{n+1}&-\tfrac{\alpha}{n+1}&\tfrac2{n+1}&0\\\tfrac{\alpha}{n}&\tfrac{2\beta}{n}&0&\tfrac2{n}\\1&0&0&0\\0&1&0&0\end{pmatrix}\begin{pmatrix}p(n)\\q(n)\\p(n-1)\\q(n-1)\end{pmatrix}$$

A Mathematica implementation would use NestList[]:

n = 0; nmax = 11;
res = Take[NestList[(n++;
                     Simplify[{{2 β/(1 + n), -α/(1 + n), 2/(1 + n), 0},
                               {α/n, 2 β/n, 0, 2/n},
                               {1, 0, 0, 0},
                               {0, 1, 0, 0}}.#]) &,
                    {(α^2 + 4 β^2)/α, 0, -2 β/α, -1}, nmax - 1], All, 2]

Note that the eleventh term is already quite complex:

LeafCount /@ Last[res]
   {154, 118}