# Nonlinear recurrence relation

I entered the following recurrence relation into RSolve and it just gave the question back to me. $$a_{n+1} = \frac{(1+a_n +(a_{n−1})^3)} 3$$

Is there another way to get a formula for the nth term?

• indeed RSolve[{a[n + 1] == (1 + a[n] + a[n - 1]^3 )/3, a[1] == 0, a[2] == 0}, a[n], n] Apr 14 '13 at 6:44
• It is generally difficult to obtain explicit solutions for nonlinear recurrences. Do you have any reason why you're supposed to expect a closed form solution? Apr 14 '13 at 6:51
• There is an undocumented function, SequenceLimit[], which you might be interested in, if what you want is to numerically estimate the limit of your sequence. Apr 14 '13 at 7:09
• Well you can always turn it into a system of first order difference equations and get the equilibrium points and their stable and unstable manifolds around them. This will give you a better picture of the dynamics of the equation. Apr 14 '13 at 7:18
• @Spawn1701D. Thanks for the tip. That sounds promising. Can you tell me where to find an explanation of such methods.
– jim
Apr 14 '13 at 12:37

In my process of learning, I am always encountering many recursive formular. The answer as shown below is my summarization about dealing with these recursive formular in Mathematica's functional paradigm.

For the general formular $$a_{n+1}=func(a_{n-1},a_n)$$ you can use the a general solution

# & @@@NestList[{#2, $func$[#1,#2]} & @@ # &, {$a_1$, $a_2$}, $n$]

where $n$ is the $n$th number

So in your question :$$func(x,y)=\frac{x^3+y+1}{3}$$

# & @@@
NestList[{#2, 1/3 (#1^3 + 1 + #2)} & @@ # &, {1, 2}, 5]

{1, 2, 4/3, 31/9, 184/81, 32176/2187, 14579713/1594323}


### Other example

Fibonacci

$a_{n+1}=a_{n-1}+a_n$

# & @@@
NestList[{#2, #1+#2} & @@ # &, {1, 1}, 10]

 {1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89}


Double sequences

$$x_n=\frac{1}{2}\left(x_{n-1}+\sqrt{1-y_{n-1}{}^2}\right) \\ y_n=\frac{1}{2}\left(y_{n-1}+\sqrt{1-x_n{}^2}\right)$$

So you can do like this:

func[x_, y_] := 1/2 (Sqrt[1 - y^2] + x);
NestList[
Function[{x, y},
{func[x, y], func[y, func[x, y]]}] @@ # &, {.5, Sqrt[3]/4}, 10]

{{0.5, Sqrt[3]/4}, {0.700694, 0.573237}, {0.760042, 0.611556},
{0.775621, 0.621377}, {0.779567, 0.623848}, {0.780556, 0.624467},
{0.780804, 0.624622}, {0.780865, 0.624661}, {0.780881, 0.62467},
{0.780885, 0.624673}, {0.780886, 0.624673}}


Or

NestList[
{1/2 (Sqrt[1 - #2^2] + #1),
1/2 (#2 + Sqrt[1 - (1/2 (Sqrt[1 - #2^2] + #1))^2])} & @@ # &, {.5, Sqrt[3]/4}, 10]


Related

### Reference

One way to approach this is to write the equation recursively:

a[n_] := a[n] = (1 + a[n - 1] + a[n - 2]^3)/3;
a[1] = a1;
a[0] = a0;


This leaves the initial conditions in terms of two generic parameters a0 and a1. For example, a[3] gives

1/3 (1 + a1^3 + 1/3 (1 + a0^3 + a1))


FullSimplify[a[4]] is:

1/81 (27 + (1 + a0^3 + a1)^3 + 3 (4 + a0^3 + a1 + 3 a1^3))


and so on. This gives an expression for a[n] expanded in terms of the initial conditions. If you wish a numerical answer, assign a0 and a1 to specific values.

Just for fun: exploring critical points

f[x_, y_] := (1 + y + x^3)/3
nf[a0_, a1_, n_] :=
NestList[{#[[2]], f[#[[1]], #[[2]]]} &, {a0, a1}, n]
g = {#1, #2, f[#1, #2]} &;
cp = {a, b} /. Solve[nf[a, b, 1] == {a, b}, {a, b}];


A limited exploration:

vis[a0_, a1_, n_] := With[{pt = g @@@ nf[a0, a1, n]},
Show[Plot3D[f[x, y], {x, -2, 2}, {y, -2, 2}, PlotRange -> {-10, 10},
Mesh -> None, PlotStyle -> Opacity[0.4]],
Graphics3D[{Black, PointSize[0.03], Point[{a0, a1, f[a0, a1]}],
PointSize[0.01], Red, Point[pt], Line[pt], Green,
PointSize[0.03], Point[g @@@ cp]}]]
]