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]

Addendum
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
.