# 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. – mikado Sep 3 '16 at 7:24
• please show your efforts, tell us what you've tried. :) – Wjx Sep 3 '16 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). – Szabolcs Sep 6 '16 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. – richmlow Sep 5 '16 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 – Vaclav Kotesovec Sep 3 '16 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 – richmlow Sep 5 '16 at 8:25
• ubpdqn, Thank you for your reply. It is helping me to get to the final goal. – richmlow Sep 5 '16 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. – bbgodfrey Sep 7 '16 at 13:42