# How to plot a system of recurrence equations

I have the following code working well ... but I would then like to plot M[k] and L[k] on the same axes with k along the horizontal and M[k] and L[k] values along the vertical. The graph I'm looking to automatically generate from the code is below. I've tried various constructs of DiscretePlot and ListPlot but to no avail ... I'd appreciate suggestions!!

Clear[a, b, k, M, L, M0, L0, NN]
NN = 25; a = 2/10; b = 1/10; M0 = 30; L0 = 50;
N[TableForm[
Prepend, {RecurrenceTable[{M[k + 1] == (1 - a)*M[k] + a*L[k],
L[k + 1] == (1 - b)*L[k] + b*M[k], M[0] == M0, L[0] == L0}, {M,
L}, {k, 0, NN}], Range[0, NN]}],
TableHeadings -> {{}, {"k", "M[k]", "L[k]"}}]]


Straightforward method:

With[{NN = 25, a = 2/10, b = 1/10, M0 = 30, L0 = 50},
ListLinePlot[Transpose[RecurrenceTable[{M[k + 1] == (1 - a) M[k] + a L[k],
L[k + 1] == (1 - b) L[k] + b M[k],
M[0] == M0, L[0] == L0},
{M, L}, {k, 0, NN}]],
DataRange -> {0, NN}, PlotRange -> All]]


Slick method:

With[{NN = 25, a = 2/10, b = 1/10, M0 = 30, L0 = 50},
ListLinePlot[Transpose[NestList[{{1 - a, a}, {b, 1 - b}}.# &, {M0, L0}, NN]],
DataRange -> {0, NN}, PlotRange -> All]]


Even slicker method:

With[{NN = 25, a = 2/10, b = 1/10, M0 = 30, L0 = 50},
DiscretePlot[MatrixPower[{{1 - a, a}, {b, 1 - b}}, k, {M0, L0}] // Evaluate,
{k, 0, NN}, Filling -> None, Joined -> True,
PlotRange -> All]]


All three versions produce the following figure:

• wonderfully creative, J. M. ... can DiscretePlot also be used in some way to show the points along the curve? Mar 29, 2020 at 15:25
• Sure, just add the option setting PlotMarkers -> {"Point", Medium} to DiscretePlot[]. Mar 29, 2020 at 15:28
• Thank you very much, JM, very clever coding!!!! Mar 29, 2020 at 15:36
Clear["Global*"]

eqns = {
M[k + 1] == (1 - a)*M[k] + a*L[k],
L[k + 1] == (1 - b)*L[k] + b*M[k],
M[0] == M0, L[0] == L0};


RSolve provides the exact solution to the recurrence equations

sol = RSolve[eqns, {L, M}, k][[1]]

{* {L -> Function[{k}, (a L0 + (1 - a - b)^k b L0 + b M0 - (1 - a - b)^k b M0)/(
a + b)], M ->
Function[{k}, -((-a L0 + a (1 - a - b)^k L0 - a (1 - a - b)^k M0 - b M0)/(
a + b))]} *}


Verifying,

eqns /. sol // Simplify

{* {True, True, True, True} *}


solEx[k_] = {L[k], M[k]} /. sol /.
{a -> 2/10, b -> 1/10, M0 -> 30, L0 -> 50} // Simplify

{* {1/3 10^(1 - k) (2 7^k + 13 10^k), 1/3 10^(1 - k) (-4 7^k + 13 10^k)} *}


The functions share a common limit

lim = Limit[solEx[k], k -> Infinity]

(* {130/3, 130/3} *)


Plotting

With[{NN = 25},
Plot[Evaluate@solEx[k], {k, 0, NN},
PlotRange -> All,
PlotLegends -> Placed[{"L", "M"}, {0.5, 0.3}],
Prolog -> {Gray, Dashed,
Line[{{0, lim[[1]]}, {NN, lim[[1]]}}]}]]


EDIT: If RSolve is unable to solve:

Clear["Global*"]

NN = 25; a = 2/10; b = 1/10; M0 = 30; L0 = 50;
rt = RecurrenceTable[{M[k + 1] == (1 - a)*M[k] + a*L[k],
L[k + 1] == (1 - b)*L[k] + b*M[k], M[0] == M0, L[0] == L0}, {M, L}, {k, 0,
NN}];

M[k_Integer] := rt[[k + 1, 1]] /; 0 <= k <= NN

L[k_Integer] := rt[[k + 1, 2]] /; 0 <= k <= NN

DiscretePlot[{L[k], M[k]}, {k, 0, NN},
PlotRange -> All, Filling -> None,
PlotLegends -> Placed["Expressions", {0.5, 0.3}]]


• Hi Bob: This is very nice and clever. The convergence line is a nice feature; when RSolve can't solve a system, I resort to the RecurrenceTable. How would sol = work when storing values from a recurrence table? Mar 29, 2020 at 18:23
• sol[k_] := table[[k]] Mar 29, 2020 at 18:31
• but, what if I started with the RecurrenceTable instead of RSolve, as RSolve doesn't always produce a solution? Mar 29, 2020 at 18:47
• Bob: Many thanks for the edit -- very helpful and I've learned much. Thanks for sharing your talents. Mar 30, 2020 at 1:26
• Corrected index in definitions of L[k] and M[k] Mar 30, 2020 at 2:04