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

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.

• You need memoization. Check out a similar problem here. And, does you parameters, $\alpha$ etc. have any numeric values? Oct 3 '19 at 0:45

As mentioned by corey979, memorization will speed up the calculations

Clear["Global*"]

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

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

AbsoluteTiming[p;]

(* {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;]

(* {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)) *)

• Thanks a lot for your help. Oct 3 '19 at 21:51

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

With[{m = 3},
Clear[p, q];
p = -2 β/α;
p = (α^2 + 4 β^2)/α;
q = -1;
q = 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.

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}
`