# Speed up iteration of pair of fully coupled recursion equations

I am trying to do several hundred iterations of these fully coupled recursion equations (but will use twenty iterations for the example).

ClearAll["Global*"]

x1[t_] := x1[t] = (1 - m) x1[t - 1] + m x2[t - 1]
x2[t_] := x2[t] = (1 - m) x2[t - 1] + m x1[t - 1]

x1[0] := 1
x2[0] := 0

ListPlot[Table[{{t, x1[t]}, {t, x2[t]}} /. m -> 0.01, {t, 0,
20}]] // AbsoluteTiming


Memoization sped it up from 27.933541 seconds to 7.681030 seconds but I think it takes so long because it is still recursively calculating all the values from the other equation. I have read several posts about partially coupled recursion equations (e.g. Solve pair of recurrence relations, How do I use RSolve to solve a system of recurrence relations?) but ultimately, I won't be able to use RSolve because my equations will be too complicated and I cannot find a solution to speed up iterations of fully coupled recursion equations.

Is there a way to link x1[t] and x2[t] together for memoization? Or is a For loop the way to go? Thanks!

Adding a Print on the most right-hand side of your definitions:

ClearAll[x1, x2]

x1[t_] := x1[t] = (
Print["Computation of x1 for t = ", t];
(1 - m) x1[t - 1] + m x2[t - 1]
);

x2[t_] := x2[t] = (
Print["Computation of x2 for t = ", t];
(1 - m) x2[t - 1] + m x1[t - 1]
);

x1[0] := 1
x2[0] := 0


and evaluating the table with an iterator that goes up to 3, we get

Table[{{t, x1[t]}, {t, x2[t]}} /. m -> 0.01, {t, 0, 3}]

During evaluation of In[..]:= Computation of x1 for t = 1
During evaluation of In[..]:= Computation of x2 for t = 1
During evaluation of In[..]:= Computation of x1 for t = 2
During evaluation of In[..]:= Computation of x2 for t = 2
During evaluation of In[..]:= Computation of x1 for t = 3
During evaluation of In[..]:= Computation of x2 for t = 3

(* {
{{0, 1}, {0, 0}}, {{1, 0.99}, {1, 0.01}},
{{2, 0.9802}, {2, 0.0198}}, {{3, 0.970596}, {3, 0.029404}}
}
*)


so the computation is done once for each of the x's at each step, which shows that values at previous steps are indeed stored and used.

A possibility to speed up the computation would be to avoid using ReplaceAll and define your functions with m as argument:

ClearAll[x1, x2]

x1[t_, m_] := x1[t, m] = (1 - m) x1[t - 1, m] + m x2[t - 1, m]
x2[t_, m_] := x2[t, m] = (1 - m) x2[t - 1, m] + m x1[t - 1, m]

x1[0, m_] := 1
x2[0, m_] := 0

ListPlot[Table[{{t, x1[t, 0.01]}, {t, x2[t, 0.01]}}, {t, 0, 20}]] // AbsoluteTiming


On my computer, the timing goes from about 7 seconds, to less than 0.05 second.

• Thanks! The following also speeds it up on a similar order of magnitude. Subbing in a series of parameters in the recursion definition, rather than taking a parameter value as a function argument, will be helpful to me (and possibly other users) down the road as I get into more complicated equations. parameters = {m->0.01} but can be parameters = {m->0.01, s->0.1, alpha ->0.5} x1[t_] := x1[t] = (1 - m) x1[t - 1] + m x2[t - 1]/.parameters x2[t_] := x2[t] = (1 - m) x2[t - 1] + m x1[t - 1]/.parameters x1[0] := 1 x2[0] := 0 – biologyUser Sep 15 '16 at 22:00
• I agree, this code avoids one to rewrite all the definitions. It is better. In case you have many parameters, and more complicated expressions (with more leaves), it is possible that the argument-approach performs better, but I doubt there will be an important difference in efficiency with your new code. (This could make another series of tests! :-)) – user31159 Sep 15 '16 at 22:09
• Faster, shorter, and easier to read: DiscretePlot[{x1[t, 0.01], x2[t, 0.01]}, {t, 0, 20}, Filling -> None] – Karsten 7. Sep 16 '16 at 8:03

RSolve will provide an exact closed-form solution that is faster

Clear[soln]

soln[m_, t_] = RSolve[{
x1[t] == (1 - m) x1[t - 1] + m x2[t - 1],
x2[t] == (1 - m) x2[t - 1] + m x1[t - 1],
x1[0] == 1, x2[0] == 0}, {x1[t], x2[t]}, t][[1]]


Plot[
Evaluate[
{x1[t], x2[t]} /. soln[1/100, t]],
{t, 0, 10},
PlotLegends -> "Expressions"]


Could you not use RecurrenceTable?

m = 0.01;
res = RecurrenceTable[{
x1[t] == (1 - m) x1[t - 1] + m x2[t - 1],
x2[t] == (1 - m) x2[t - 1] + m x1[t - 1],
x1[0] == 1, x2[0] == 0}, {x1, x2}, {t, 20}]


is superfast. Here's an ugly way to recreate your plot:

ListPlot[{
Transpose[{Table[t, {t, 0, 20}], Transpose[res][[1]]}],
Transpose[{Table[t, {t, 0, 20}], Transpose[res][[2]]}]
}]
`