Delay difference equation

Question

How can I correctly specify and solve a system of delay difference equations in mathematica similar to NDSolve with regard to delay differential equations?

For example if I want to "see" the values explicitly with RecurrenceTable and my system of equations has the following form:

a[n+1] =b[n]+c[n-1], b[n+1]= a[n-1]+c[n], c[n+1]= b[n-1]-a[n]

I am not sure how to specify the initial conditions like it is done with a initial history function with delay differential equations like the following system for example:

NDSolve[{a'[t] == -a[t] + c[t - 1], b'[t] == -b[t] + a[t-1], 
  c'[t] == -c[t] + b[t - 1], a[t /; t <= 0] == 1, b[t /; t <= 0] == 2,
   c[t /; t <= 0] == 3}, {a, b, c}, {t, 0, 20}]

I hope it is clear what I mean and hope for some helpful hints :)

Answer 1

You can just write the equations directly:

a[n_] := a[n] = b[n - 1] + c[n - 2];
b[n_] := b[n] = a[n - 2] + c[n - 1];
c[n_] := c[n] = b[n - 2] - a[n - 1];
a[1] = 1; a[0] = 0.5;
b[1] = 1; b[0] = 0.25;
c[1] = 2; c[0] = 1.5;

I've shifted the time indices, and arbitrarily assigned some initial conditions. To see individual values you can just type:

a[20]

To see a range of values:

a[#] & /@ Range[20]

Answer 2

The solution provided by bill s is analogous to NDSolve. A solution analogous to DSolve is

First@RSolve[{a[n] == b[n - 1] + c[n - 2], 
    b[n] == a[n - 2] + c[n - 1], 
    c[n] == b[n - 2] - a[n - 1], 
    a[1] == 1, a[0] == 0.5, 
    b[1] == 1, b[0] == 0.25, 
    c[1] == 2, c[0] == 1.5}, {a[n], b[n], c[n]}, n];
s = Chop[FullSimplify[%, n ∈ Integers]
(* {a[n] -> 1.25 - 0.416667 (-1)^n - 1.58333 Cos[(n π)/3] + 
        1.25 Cos[(2 n π)/3] + 1.29904 Sin[(n π)/3] - 0.433013 Sin[(2 n π)/3], 
    b[n] -> 1.25 - 1. Cos[(2 n π)/3] - 0.866025 Sin[(2 n π)/3], 
    c[n] -> -0.416667 (-1)^n + 1.91667 Cos[(n π)/3] + 0.721688 Sin[(n π)/3]} *)

which gives analytic solutions. Note that, if the initial conditions (chosen to match those of bill s) were omitted, the solution would be in terms of six arbitrary constants, similar to the solution of three coupled second order ODEs.

To obtain numerical values, define

fa[m_] := a[n] /. s /. n -> m; 
fb[m_] := b[n] /. s /. n -> m; 
fc[m_] := c[n] /. s /. n -> m;
{fa[20], fb[20], fc[20]}
(* {2.5, 2.5, -0.75} *)

Answer 3

Just another approach:

mat = {
   {0, 0, 1, 0, 0, 1},
   {1, 0, 0, 0, 0, 0},
   {0, 1, 0, 0, 1, 0},
   {0, 0, 1, 0, 0, 0},
   {-1, 0, 0, 1, 0, 0},
   {0, 0, 0, 0, 1, 0}
   };
in = {1, 0.5, 1, 0.25, 2, 1.5};
abc[0] := in[[{2, 4, 6}]];
abc[n_] := Nest[mat.# &, in, n - 1][[{1, 3, 5}]];
a[n_] := abc[n][[1]];
b[n_] := abc[n][[2]];
c[n_] := abc[n][[3]];

Visualizing:

ListPlot[Transpose[abc[#] & /@ Range[20]], Joined -> True, 
 Frame -> True, PlotLegends -> {"a", "b", "c"}]