# Solving a system of recurrence relations using Mathematica

Can Mathematica solve (for b[k] and c[r]) the following system of recurrence relations? b[k] is defined for only odd natural numbers, and c[r] is defined for only even non-negative integers.

For n greater than or equal to 3;

b[2n+1] = c[2n] + c[2n-2] + 2*a[n],
c[2n] = b[2n-1] + c[2n-2],
a[n] = .0794367(.460811)^n + .255972(-.675131)^n + .664591(3.21432)^n;

c[0]=1, b[1]=4, c[2]=5, b[3]=10, c[4]=15, b[5]=34.


I seemed to have tried "everything" in Mathematica 11, but nothing seems to work.... Any help is greatly appreciated!

UPDATE (September 4, 2016)

After consolidating everything, here is the question.

Find a "closed-form" solution (formula) to the following piecewise defined recurrence relation:

For integers k greater than or equal to 3:

b[k] = b[k-1] + b[k-3] + 2*[.0794367(.460811)^((k-1)/2) + \
.255972(-.675131)^((k-1)/2) + .664591(3.21432)^((k-1)/2)]


If k is odd:

b[k] = b[k-1] + b[k-2], if k is even; b[0] = 1, b[1] = 4, b[2] = 5


This piecewise defined recurrence relation is the one that I am interested in. I verified (using this piecewise defined recurrence relation, by hand calculations) the correct values:

b[0]=1, b[1]=4, b[2]=5, b[3]=10, b[4]=15, b[5]=34, b[6]=49,
b[7]=108, b[8]=157, b[9]=348, b[10]=505, ...


Sincerely, Richard M. Low

• Have you tried RSolve? What syntax did you use? Many problems are caused by subtle errors in syntax that someone could spot easily. Commented Sep 3, 2016 at 7:24
– Wjx
Commented Sep 3, 2016 at 8:50
• Please cross reference your posts here with the corresponding ones on Wolfram Community: community.wolfram.com/groups/-/m/t/917888 Just edit your post and include a link (on both websites). This is to prevent duplication of effort (i.e. two people posting the same answer in two places and not being aware of it). Commented Sep 6, 2016 at 10:00

Eliminating b from the first two expressions in the question and redefining n yields

c[n] == 2 c[n - 1] + c[n - 2] + 2 a[n - 1]


Although RSolve can handle this expression with a as given in the question, it seems better to Rationalize it.

794367 10^-7 (460811 10^-6)^n+255972 10^-6 (-675131 10^-6)^n+664591 10^-6 (321432 10^-5)^n


Then,

Simplify[RSolve[{c[n] == 2 c[n - 1] + c[n - 2] + 2 (794367 10^-7 (460811 10^-6)^(n - 1) +
255972 10^-6 (-675131 10^-6)^(n - 1) + 664591 10^-6 (321432 10^-5)^(n - 1)),
c[0] == 1, c[1] == 5}, c[n], n], n ∈ Integers] // First

(* c[n] -> *)


The first several values are

Table[Simplify[c[n] /. %], {n, 0, 8}]

(* {1, 5, 74999988159517/5000000000000, 244999948333156036527/5000000000000000000,
784999809145104350538564757/5000000000000000000000000,
2524999383279806241504389904195767/5000000000000000000000000000000, (4^(-4 - 3 n) 5^(-7 - 6 n) (10^(6 n) (1172409346315586312019869378750000 -
1088292064380910194582788537056291 Sqrt[2]) (1 - Sqrt[2])^n + 10^(6 n) (1 + Sqrt[2])^
n (1172409346315586312019869378750000 + 1088292064380910194582788537056291 Sqrt[2]) +
1500000 (130809072583670204544034435447 240^n 13393^(1 + n) -
110652256667034006265872216 460811^(1 + n) - 756091031029961199112699840 (-1)^n 675131^
(1 + n))))/89285790105620659634534766727097
8114998166092225878868708271269265029997/5000000000000000000000000000000000000,
26084994986029714868240982411237137250977319807/
5000000000000000000000000000000000000000000,
83844987673045565789165365810827032895161630505820037/
5000000000000000000000000000000000000000000000000} *)


or, numericized,

N[%, 15]

(* {1.00000000000000, 5.00000000000000, 14.9999976319034, 48.9999896666312,
156.999961829021, 504.999876655961, 1622.99963321845, 5216.99899720594,
16768.9975346091} *)


That these numbers are very close to being integers suggests that a more precise expression for a would yield integers directly.

Edit: Computation of b

Obtaining b from c is straightforward. For convenience, define

fc[n_] := (4^(-4 - 3 n) 5^(-7 - 6 n) (10^(6 n) (1172409346315586312019869378750000 -
1088292064380910194582788537056291 Sqrt[2]) (1 - Sqrt[2])^n + 10^(6 n) (1 + Sqrt[2])^
n (1172409346315586312019869378750000 + 1088292064380910194582788537056291 Sqrt[2]) +
1500000 (130809072583670204544034435447 240^n 13393^(1 + n) -
110652256667034006265872216 460811^(1 + n) - 756091031029961199112699840 (-1)^n 675131^
(1 + n))))/89285790105620659634534766727097


Then

fb[n_] := Simplify[fc[n] + fc[n - 1] + 2 (794367 10^-7 (460811 10^-6)^n +
255972 10^-6 (-675131 10^-6)^n + 664591 10^-6 (321432 10^-5)^n)]


The two can be interleaved, as requested, by

fh[n_] := If[OddQ[n], fb[(n - 1)/2], fc[n/2]]
Table[fh[n], {n, 0, 11}] // N

(* {1., 4., 5., 10., 15., 34., 49., 108., 157., 348., 505., 1118.} *)

ListLogPlot[Table[fh[n], {n, 0, 30}] // N, PlotRange -> All]


Another approach is to leave a undefined at first.

Simplify[RSolve[{c[n] == 2 c[n - 1] + c[n - 2] + 2 a[n], c[0] == 1, c[1] == 5},
c[n], n], n ∈ Integers] // First;
Collect[% /. {K[1] -> k, K[2] -> k}, {a[0], a[1], _Sum}, Simplify[#, n ∈ Integers] &];
ReleaseHold[Hold[Evaluate[%]] /. {k, i1_, i2_} -> {i, i1 + 2, i2 + 2} /.
k -> k - 2 /. i -> k];
Simplify[% /. z1_ Sum[z2_, z3_] :> z1 Sum[z2, z3] + (z1 z2 /. k -> 0) + (z1 z2 /. k -> 1)]
/. {k, 0, n + 1} -> {k, 2, n + 1}

(* {c[n] -> ((1 - Sqrt[2])^n - 2*Sqrt[2]*(1 - Sqrt[2])^n + (1 + Sqrt[2])^n +
2*Sqrt[2]*(1 + Sqrt[2])^n + 2*(1 + Sqrt[2])^n*
Sum[((-1)^(2 - k)*(1 - Sqrt[2])^(-2 + k)*(-1 + Sqrt[2])*a[k])/Sqrt[2], {k, 2, 1 + n}]
+ 2*(1 - Sqrt[2])^n*Sum[((-1)^(2 - k)*(1 + Sqrt[2])^(-1 + k)*a[k])/
Sqrt[2], {k, 2, 1 + n}])/2} *)


Results can be tested at each step of the simplification, in all cases yielding

Table[Simplify[c[n] /. % /. n -> i], {i, 0, 5}]

(* {1, 5, 11 + 2 a[2], 27 + 4 a[2] + 2 a[3], 65 + 10 a[2] + 4 a[3] + 2 a[4],
157 + 24 a[2] + 10 a[3] + 4 a[4] + 2 a[5]} *)


It is natural to ask whether the coefficients of a represent a known series, and they do.

2 FindSequenceFunction[Cases[Last@%, z_ a[_] -> z/2] // Reverse][n]

(* 2 Fibonacci[n, 2] *)


Thus, c[n] can be represented as

f[n_] := Simplify[1/2 ((1 - Sqrt[2])^n - 2 Sqrt[2] (1 - Sqrt[2])^n + (1 + Sqrt[2])^n +
2 Sqrt[2] (1 + Sqrt[2])^n)] + Plus @@ (2 Fibonacci[n + 1 - #, 2] a[#] & /@ Range[2, n])


One might hope that a[n], rounded to the nearest integer, also represents a known series, but, FindSequenceFunction returns unevaluated.

Table[Round@a[n], {n, 0, 6}] // FindSequenceFunction

(* FindSequenceFunction[{1, 2, 7, 22, 71, 228, 733}] *)


Second Addendum: a, b, and c in closed form

Vaclav Kotesovec made the very useful observation in a comment, that a[n], rounded to integers, has the recurrence relation,

a[n] == 3 a[n - 1] + a[n - 2] - a[n - 3]


It turns out that Mathematica can produce this result without difficulty, once it is realized from Kotesovec's comment that it is feasible to do so. With a[n] as given above,

Table[FullSimplify[a[n]], {n, 0, 11}] // N // Round // FindLinearRecurrence
(* {3, 1, -1} *)


This can be solved to obtain a in closed form.

Clear[a]
a[n] /. First@RSolve[{a[n] == 3 a[n - 1] + a[n - 2] - a[n - 3],
a[0] == 1, a[1] == 2, a[2] == 7}, a[n], n]
(* Root[1 - #1 - 3 #1^2 + #1^3 &, 2]^n Root[-1 + 18 #1 - 74 #1^2 + 74 #1^3 &, 1] +
Root[1 - #1 - 3 #1^2 + #1^3 &, 1]^n Root[-1 + 18 #1 - 74 #1^2 + 74 #1^3 &, 2] +
Root[1 - #1 - 3 #1^2 + #1^3 &, 3]^n Root[-1 + 18 #1 - 74 #1^2 + 74 #1^3 &, 3] *)


Sample results are as expected,

Table[FullSimplify[%], {n, 0, 11}]
(* {1, 2, 7, 22, 71, 228, 733, 2356, 7573, 24342, 78243, 251498} *)


This exact result can be numericized to whatever precision is desired using N. For \$MachinePrecision,

%% // N
(* 0.255972 (-0.675131)^n + 0.0794367 0.460811^n + 0.664591 3.21432^n *)


which is the approximate expression for a shown in the question. Now, all this would be of little interest, except that the same can be done for b and c!

Table[Simplify[fc[n]], {n, 0, 11}] // N // Round // FindLinearRecurrence
(* {3, 1, -1} *)
c[n] /. First@RSolve[{c[n] == 3 c[n - 1] + c[n - 2] - c[n - 3],
c[0] == 1, c[1] == 5, c[2] == 15}, c[n], n]
(* Root[1 - #1 - 3 #1^2 + #1^3 &, 1]^n Root[-1 - 25 #1 - 37 #1^2 + 37 #1^3 &, 1] +
Root[1 - #1 - 3 #1^2 + #1^3 &, 2]^n Root[-1 - 25 #1 - 37 #1^2 + 37 #1^3 &, 2] +
Root[1 - #1 - 3 #1^2 + #1^3 &, 3]^n Root[-1 - 25 #1 - 37 #1^2 + 37 #1^3 &, 3] *)
Table[FullSimplify[%], {n, 0, 11}]
(* {1, 5, 15, 49, 157, 505, 1623, 5217, 16769, 53901, 173255, 556897} *)
%% // N
(* -0.428786 (-0.675131)^n - 0.0428314 0.460811^n + 1.47162 3.21432^n *)

Table[Simplify[fb[n]], {n, 0, 11}] // N // Round // FindLinearRecurrence
(* {3, 1, -1} *)
b[n] /. First@RSolve[{b[n] == 3 b[n - 1] + b[n - 2] - b[n - 3],
b[0] == 4, b[1] == 10, b[2] == 34}, b[n], n]
(* Root[1 - #1 - 3 #1^2 + #1^3 &, 2]^n Root[-2 + 90 #1 - 148 #1^2 + 37 #1^3 &, 1] +
Root[1 - #1 - 3 #1^2 + #1^3 &, 1]^n Root[-2 + 90 #1 - 148 #1^2 + 37 #1^3 &, 2] +
Root[1 - #1 - 3 #1^2 + #1^3 &, 3]^n Root[-2 + 90 #1 - 148 #1^2 + 37 #1^3 &, 3] *)
Table[FullSimplify[%], {n, 0, 11}]
(* {4, 10, 34, 108, 348, 1118, 3594, 11552, 37132, 119354, 383642, 1233148} *)
N // %%
(* 0.718273 (-0.675131)^n + 0.0230942 0.460811^n + 3.25863 3.21432^n *)


Thus, we have been able to provide in a straightforward manner both exact and MachinePrecision expressions for a, b, and c.

• bbgodfrey, Your answer is very helpful! Thank you very much. Commented Sep 5, 2016 at 22:32

If numerical values are desired:

a[n_] := .0794367 (.460811)^n + .255972 (-.675131)^
n + .664591 (3.21432)^n;
b[j_] := c[j - 1] + c[j - 3] + 2 a[(j - 1)/2]
c[j_] := b[j - 1] + c[j - 2]
c[0] := 1
b[1] := 4
c[2] := 5
b[3] := 10
c[4] := 15
b[5] := 34.


Presentation:

fun[j_] := If[EvenQ[j], c[j], b[j]];
bcf[n_?OddQ] :=
With[{r = {0 -> Red, 1 -> Blue}},
Framed[TableForm[{Style[fun@#, Mod[#, 2] /. r], a@#} & /@
Range[0, n],
Subscript[#1,
Style[ToString[#2[[1]] - 1], Mod[#2[[1]] - 1, 2] /. r]] &,
Flatten[Table[{Style["c", Red], Style["b", Blue]}, (n + 1)/
2]]], {"value", "\!$$\*SubscriptBox[\(a$$, $$n$$]\)"}}]]]


This does not scale well. I have not tried FindSequenceFunction etc. Others may wish to.

Apologies for any errors.

• Sequence a(n) is oeis.org/A030186 Commented Sep 3, 2016 at 13:29
• Everybody has been quite helpful thus far! Thank you. My original question can be formulated in the following manner: Can Mathematica solve (ie., find a "closed-form" formula for b[k]) for the following "piecewise" defined recurrence relation? For integers k greater than or equal to 3, b[k] = b[k-1] + b[k-3] + 2*[.0794367(.460811)^((k-1)/2) + .255972(-.675131)^((k-1)/2) + .664591(3.21432)^((k-1)/2)], if k is odd; b[k] = b[k-1] + b[k-2], if k is even; b[0] = 1, b[1] = 4, b[2] = 5. Sincerely, Richard M. Low Commented Sep 5, 2016 at 8:25
• ubpdqn, Thank you for your reply. It is helping me to get to the final goal. Commented Sep 5, 2016 at 22:38
• @VaclavKotesovec Thank you for identifying the sequence a[n]. Knowing its recurrence relation, I then was able to determine the recurrence relations with Mathematica for b[n] and c[n]`. Very helpful. Commented Sep 7, 2016 at 13:42