7

I have a set of coupled recursion equations and want to express the final solution in terms of one of the variables:

RSolve[
 {d[n] == a[n] + b[n] + c[n],
  a[n] == a[n - 1] + c[n - 1],
  b[n] == a[n - 1],
  c[n] == b[n - 1] + c[n - 1]},
 {d[n]}, n]

I know the result can be expressed as:

d[n] = d[n-1] + d[n-2] + d[n-4]

but Mathematica gives a variety of error messages, depending upon the variables. (I get the error even if I include the following conditions: a[1] == a[2] == b[1] == b[2] == c[1] == c[2] == 0.) I've also included obvious equivalent conditions such as d[n-1] = a[n-1]+b[n-1]+c[n-1] and so on.

How can Mathematica solve these recursion relations and express the final answer solely in terms of d[]?

David G. Stork
  • 41,180
  • 3
  • 34
  • 96

3 Answers3

6

I apologize if this is unhelpful (and it is not efficient):

m = {{1, 0, 1}, {1, 0, 0}, {0, 1, 1}}
d[n_] := MatrixPower[m, n].{a, b, c}.{1, 1, 1}
df = FindSequenceFunction[d /@ Range[10]] /. {a -> 1, b -> 1, c -> 1};

yields a DifferenceRoot.

Testing one instance:

FullSimplify[d[n] - d[n - 1] - d[n - 2] - d[n - 4]]

yields 0.

Obviously, this is only 1 instance but perhaps it prompts better ideas.

ubpdqn
  • 60,617
  • 3
  • 59
  • 148
6

Without setting initial conditions. Looking for linear recurrences:

rs = First@RSolve[{
      a[n] == a[n - 1] + c[n - 1],
      b[n] == a[n - 1],
      c[n] == b[n - 1] + c[n - 1],
      a[1] == a1, c[1] == c1, b[1] == b1}, {a@n, b@n, c@n},  n] ;

ss = Simplify[Table[a[n] + b[n] + c[n] /. rs, {n, 0, 9}] // N // Chop] // Chop;
ss1 = ss /. x_Real :> Round[x]
FindLinearRecurrence@ss1

(* {2, -1, 1} *)

So we found a different linear relationship than yours ;)

d[n] == 2 d[n-1] - d[n-2] +  d[n-3]
Dr. belisarius
  • 115,881
  • 13
  • 203
  • 453
  • 1
    ,belisarius if you look at the DifferenceRoot in my silly answer it matches yours, it just happens to also be true that d[n-1]+d[n-3]=d[n-4] for this initial condition:) – ubpdqn Feb 24 '16 at 07:09
  • @ubpdqn See update. I got rid of the initial conditions generalizing the result – Dr. belisarius Feb 24 '16 at 15:16
  • 1
    very nice...a useful exercise which I will contemplate after some sleep, Btw I already voted. To paraphrase John 'Hannibal' Smith: 'I love it when a code does what I want to do'...the ramblings of sleep deprived – ubpdqn Feb 24 '16 at 15:22
3

A general solution, making no assumptions about initial conditions, is as follows. (The derivation proceeds much like the corresponding derivation for a set of first order ODEs.) First, eliminate b.

eqs1 = {a[n] == a[n - 1] + c[n - 1], b[n] == a[n - 1], c[n] == b[n - 1] + c[n - 1]};
bsolv = Solve[d[n] == a[n] + b[n] + c[n], b[n]][[1, 1]];
eqs2 = eqs1 /. {bsolv, (bsolv /. n -> n - 1)};
eqs3 = Equal @@@ First@Solve[eqs2, {a[n], d[n], c[n]}];
(* {a[-1 + n] + c[-1 + n] == a[n], 
    a[-1 + n] + c[-1 + n] + d[-1 + n] == d[n], 
    a[-1 + n] + c[n] == d[-1 + n]} *)

Second, Eliminate a and c.

Eliminate[Join[{eqs3[[2]]}, (eqs3 /. n -> n - 1), (eqs3 /. n -> n - 2)], 
    {a[n - 1], a[n - 2], a[n - 3], c[n - 1], c[n - 2], c[n - 3]}];
ans = Equal @@ Solve[%, d[n]][[1, 1]]
(* d[n] == d[-3 + n] - d[-2 + n] + 2 d[-1 + n] *)

which is the result obtained by Dr. Belisarius by a different approach. As expected, it is a third order difference equation. The expression in the question is obtained by

Equal @@ Solve[(Subtract @@ ans) + (Subtract @@ ans /. n -> n - 1) == 0, d[n]][[1, 1]]

This difference equation is an order higher than necessary and, therefore, has a solution with four arbitrary constants instead of three.

More compact derivation

An alternative, more compact derivation can be achieved by treating all four equations in parallel.

{d[n] == a[n] + b[n] + c[n], a[n] == a[n - 1] + c[n - 1], 
    b[n] == a[n - 1], c[n] == b[n - 1] + c[n - 1]};

Equal @@@ First@Solve[%, {d[n], a[n], b[n], c[n]}];
Eliminate[Join[{First@%}, Flatten@Array[% /. n -> n - # &, 3]], 
    Flatten@Array[{a[n - #], b[n - #], c[n - #]} &, 4]];
Equal @@ Solve[%, d[n]][[1, 1]]

(* d[n] == d[-3 + n] - d[-2 + n] + 2 d[-1 + n] *)
bbgodfrey
  • 61,439
  • 17
  • 89
  • 156
  • 1
    thanks for a systematic/general approach. +1 of course, :) – ubpdqn Feb 24 '16 at 07:06
  • Of the several very helpful solutions (+1), I prefer this More compact derivation because it is compact and relies on Solve but not FindSequenceFunction. (One could generate values of d[i] and use that function, as below, but it wouldn't be elegant.) I thought that RSolve would be necessary (and spent much time exploring that), but apparently it is not necessary. Thanks to all! – David G. Stork Feb 24 '16 at 16:47