4
$\begingroup$

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[
  MapThread[
   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]"}}]]

enter image description here

enter image description here

$\endgroup$
10
$\begingroup$

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:

figure

| improve this answer | |
$\endgroup$
  • 1
    $\begingroup$ wonderfully creative, J. M. ... can DiscretePlot also be used in some way to show the points along the curve? $\endgroup$ – PRG Mar 29 at 15:25
  • $\begingroup$ Sure, just add the option setting PlotMarkers -> {"Point", Medium} to DiscretePlot[]. $\endgroup$ – J. M.'s technical difficulties Mar 29 at 15:28
  • $\begingroup$ Thank you very much, JM, very clever coding!!!! $\endgroup$ – PRG Mar 29 at 15:36
4
$\begingroup$
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} *}

For your specific parameters,

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

enter image description here

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

enter image description here

| improve this answer | |
$\endgroup$
  • $\begingroup$ 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? $\endgroup$ – PRG Mar 29 at 18:23
  • $\begingroup$ sol[k_] := table[[k]] $\endgroup$ – Bob Hanlon Mar 29 at 18:31
  • $\begingroup$ but, what if I started with the RecurrenceTable instead of RSolve, as RSolve doesn't always produce a solution? $\endgroup$ – PRG Mar 29 at 18:47
  • $\begingroup$ Bob: Many thanks for the edit -- very helpful and I've learned much. Thanks for sharing your talents. $\endgroup$ – PRG Mar 30 at 1:26
  • $\begingroup$ Corrected index in definitions of L[k] and M[k] $\endgroup$ – Bob Hanlon Mar 30 at 2:04

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.