6
$\begingroup$

The first two terms of a sequence are $a_1 = 1$ and $a_2 = \frac {1}{\sqrt3}$. For $n\ge1$,

$$a_{n + 2} = \frac {a_n + a_{n + 1}}{1 - a_na_{n + 1}}.$$ What is $|a_{2009}|$?

So I tried (thanks to Paul Wellin):

fun[n_] := Module[{prev = 1, this = 1/Sqrt[3], next},
  Do[next = Simplify[(prev + this)/(1 - prev this)];
   prev = this;
   this = next, {n - 1}];
  prev]

And I got:

fun[2009]

(* 0 *)

Which is the correct answer. Kind of amazed it went so quickly. I then tried:

Table[fun[i], {i, 1, 100}]

Which produces the first 100 terms of the sequence. Again, amazingly quick. Now, this turns out to be a periodic sequence. I counted with my finger, and determined 24.

Anyone have a good way to use Mathematica to determine if the sequence is periodic and find the period?

$\endgroup$
  • $\begingroup$ mathematica.stackexchange.com/questions/80163/… $\endgroup$ – Marvin Dec 5 '15 at 9:31
  • 1
    $\begingroup$ Entirely different method: Recognize the recursion as the sum of angles formula for tangent. Set $\tan \theta_{n} = a_n$ so $\theta_1 = \pi/4$, $\theta_2 = \pi/6$, and $\theta_{n+2} = \theta_n + \theta_{n+1}$. This is a much easier recurrence... $\endgroup$ – Eric Towers Dec 5 '15 at 19:08
  • $\begingroup$ @EricTowers Could you add this as an answer? I'd love to see what you do. $\endgroup$ – David Dec 6 '15 at 3:33
6
$\begingroup$

One way of doing this is starting with your list of the first 100 terms.Then ask for

zzz = SequencePosition[lst, Take[lst, -2]]

(* {{3, 4}, {27, 28}, {51, 52}, {75, 76}, {99, 100}} *)

The length of this result is larger than 1, so your list must be periodic. The period is

zzz[[2, 1]] - zzz[[1, 1]]

(* 24 *)
$\endgroup$
  • $\begingroup$ thank you introducing me to SequencePosition $\endgroup$ – ubpdqn Dec 5 '15 at 10:56
  • $\begingroup$ @FredSimons I agree, really simple explanation to help students with exploration. Thanks. $\endgroup$ – David Dec 5 '15 at 18:19
7
$\begingroup$

Started as a comment to the OP. Turning this into an answer was requested. The OP asked for Mathematica help solving the recurrence. I suggested recognizing that the recurrence was the sum of angles formula for tangent, then going from there. Consequently, this answer uses Mathematica to solve this problem using only one "external" ingredient right at the beginning. (At the end, we explore a claim made at the beginning.)

Problem: Find the period of $$\begin{align} a_1 &= 1 \\ a_2 &= \frac{1}{\sqrt{3}} \\ a_{n+2} &= \frac{a_n + a_{n+1}}{1-a_n a_{n+1}}, n \geq 1 \end{align}$$

(It would be seriously convenient if FunctionPeriod[] just did this. We could "punt" and use the accepted answer here, applying it to the straightforward implementation of the above recursion and initial values. However, that's not what I'm here to write.)

"External" fact: The recursion is the sum of angles formula for tangent.

Claim at the beginning: We could replace the recursion formula with the analog for any other sum of angles formula and get exactly the same sequence of angles. I.e., there are at least two other recursions with this same solution.

First, we check this external fact with Mathematica and get the initial values.

Simplify[Tan[a+b] == (Tan[a] + Tan[b])/(1-Tan[a]Tan[b])
(*  True  *)
ArcTan[1]
(*  Pi/4  *)
ArcTan[1/Sqrt[3]]
(*  Pi/6  *)

So we may reformulate the question as $$\begin{align} \theta_1 &= \frac{\pi}{4} \\ \theta_2 &= \frac{\pi}{6} \\ \theta_{n+2} &= \theta_1 + \theta_2 \pmod{2\pi}, n \geq 1 \end{align}$$

Sums of multiples of $\pi/4$ and $\pi/6$ have denominators dividing $12$.

LCM[4,6]
(*  12  *)

There are $24$ twelfths of $\pi$ in a circle.

2 Pi / (Pi/12)
(*  24  *)

Consequently, the period of this recurrence is at most $24^2$. (The recurrence is order $2$, so it is perhaps better to think of each adjacent pair of sequence members producing the next adjacent pair of sequence members.) (Although, as is frequently the case with homogeneous recurrences, if $\dots 0, 0 \dots$ occurs in our sequence, every subsequent angle is zero. So it is possible for the period to be less than $24$ and for there to be an initial aperiodic part, say if $\theta_1 = \theta_2 = \pi$, although it will turn out that this is not the case in this problem.)

So we generate the first $2 \times 24^2$ entries, which will contain two or more periods.

Remove[q];
q[1] = Pi/4;
q[2] = Pi/6;
q[n_ /; n >= 3] := q[n]= Mod[q[n-1]+q[n-2],2 Pi];
(* We memoize members of the sequence to avoid redundant recomputation. *)

Generate the $2 \times 24^2$ entries and glance at some of them to see if periodicity seems likely.

Table[ {k, q[k]} , {k, 24, 2*24^2, 24} ]
(*  {{ 24, 23 Pi/12 }, { 48, 23 Pi/12 }, ... , {1152, 23 Pi/12 }}  *)

Seems rather likely to be periodic. Test this with Fourier[]:

periods = Flatten[
  Position[
    Chop[Fourier[Table[q[k], {k, 1, 2*24^2}]]],
    n_ /; n != 0
  ]
]
(*  {1, 49, 97, 145, 193, 241, 289, 337, ... , 1009, 1057, 1105}  *)
GCD @@ Select[DeleteDuplicates[Flatten[
    Outer[Abs[#1 - #2] &, periods, periods]
  ]], # > 0 &]
(*  48  *)

So the period divides $48$. (Actually, we could know that it divides $24$. We ran to $2 \times 24^2$. Try it again with $3 \times 24^2$ and get $72$. In fact, if we run to $N \times 24^2$ for $N$ up to $24$, we'll get $24N$ for the GCD above.)

Let's have Mathematica try every divisor of $48$ as a candidate period.

{#, Simplify[q[3 + #] == q[3]]} & /@ Divisors[48]
(*  {{1, False}, {2, False}, {3, False}, {4, False}, 
     {6, False}, {8, False}, {12, False}, {16, False}, 
     {24, True}, {48, True}}  *)

so the only candidate periods are $24$ and $48$. (Note that starting at q[3], we're finding out about the recursion moreso than the initial values.) Then we test the shortest candidate:

q /@ Range[24] == q /@ Range[25, 2*24]

The above method works, and is valid, but perhaps one would be happier seeing that the recurrence actually recurs after $24$ steps. Starting with the symbolic first two angles $a,b$, what is the resulting angle after (a divisor of $24$) applications of the recursion rule. We frame this as working on pairs in Nest[].

TableForm[
  {#[[1]], #[[2]] === {a, b}, #[[2]]} & /@ (
    {#, PolynomialMod[
      Nest[
        (* Thinking in pairs, the new first angle is the old second 
           angle and the new second angle is the sum of the two old 
           angles. *)
        {#[[2]], #[[1]] + #[[2]]} &, 
        {a, b}, 
        #], 
    #]} & /@ 
  Divisors[24]),
  TableDepth -> 1
]
(*  { { 1, False, {0, 0}},
      { 2, False, {a + b, a}},
      { 3, False, {a + 2 b, 2 a}},
      { 4, False, {2 a + 3 b, 3 a + b}},
      { 6, False, {5 a + 2 b, 2 a + b}},
      { 8, False, {5 a + 5 b, 5 a + 2 b}},
      {12, False, {5 a, 5 b}},
      {24,  True, {a, b}} }  *)

So the recursion has a natural period of $24$. It's possible for the choice of initial conditions to cause the period to be shorter. (Recall that if the first two angles are $0,0$, then the sequence is all zeroes.) So we check for this. The resulting period still has to divide $24$, so we only check those (although if we were not certain of this, we could check every period up to $24$).

TableForm[
  { #[[1]] , #[[2]] == #[[3]] , #[[2]] , #[[3]] } & /@ (
    { # , {q[1], q[2]}, {q[1 + #], q[2 + #]} } & /@ 
      Divisors[24]
  ),
TableDepth -> 1]
(*  { { 1, False, {Pi/4, Pi/6}, {  Pi/ 6,  5 Pi/12}},
      { 2, False, {Pi/4, Pi/6}, {5 Pi/12,  7 Pi/12}},
      { 3, False, {Pi/4, Pi/6}, {7 Pi/12,    Pi   }},
      { 4, False, {Pi/4, Pi/6}, {  Pi   , 19 Pi/12}},
      { 6, False, {Pi/4, Pi/6}, {7 Pi/12,    Pi/ 6}},
      { 8, False, {Pi/4, Pi/6}, {3 Pi/ 4, 11 Pi/12}},
      {12, False, {Pi/4, Pi/6}, {  Pi/ 4,  5 Pi/ 6}},
      {24,  True, {Pi/4, Pi/6}, {  Pi/ 4,    Pi/ 6}} }  *)

And so the period is $24$. Then it's easy to compute the desired number:

Tan[q[Mod[2009,24,1]]]  (*Don't want to inadvertently try to evaluate q[0]. *)
  (* This is Tan[ q[17] ] = Tan[ Pi ]. *)
(*  0  *)

The original problem statement asked about $|a_{2009}|$ but since $0 = |0|$, there's no need to apply the absolute value function.

Returning to the claim at the beginning. Since adding angles is pretty straightforward, we can use any other sum of angles formula. This means there are (at least) two other recursion problems looking wildly different from this one (by having a wildly different looking general recursion formula), but having the same underlying computation. The sum of angles formulas for sine and cosine are, with the corresponding recursion: $$\begin{align} \sin(\theta_1 + \theta_2) &= \sin(\theta_1) \cos(\theta_2) + \cos(\theta_1) \sin(\theta_2) \\ a_{n+2} &= a_n c(a_{n+1}) \sqrt{1-a_{n+1}^2} + a_{n+1} c(a_n) \sqrt{1-a_n^2} \\ \\ \cos(\theta_1 + \theta_2) &= \cos(\theta_1)\cos(\theta_2) - \sin(\theta_1)\sin(\theta_2) \\ a_{n+2} &= a_n a_{n+1} - s(a_n)s(a_{n+1}) \sqrt{1-a_n^2}\sqrt{1-a_{n+1}^2} \end{align}$$ where the sign-tracking functions $s$ and $c$ are $s(\theta)$ is $1$ in quadrants I and II and the angle $\pi/2$ and is $-1$ otherwise, and $c$ is $1$ in quadrants I and IV and the angle $0$ and is $-1$ otherwise. We could eliminate these two sign-tracking functions by having two sequences, $q_n$ and $r_n$: $a_n \mapsto (q_n, r_n) = (\sin(\theta_n), \cos(\theta_n))$ or by switching to exponentials $a_n \mapsto \mathrm{e}^{\mathrm{i} \theta_n} = \cos(\theta_n) + \mathrm{i}\sin(\theta_n)$, which the final equality shows is really the same as switching to sines and cosines, but using real and imaginary parts of complex numbers to hold the two sequences, instead of pairs of real numbers.

Taking $\theta_1 = \pi/4$ and $\theta_2 = \pi/6$ and using either the sum of angles for sine or for cosine, we get the same sequence of angles. Consequently, we will get a final, final answer of $a_{2009} = \tan^{-1}(\sin(\theta_{2009})) = 0$ for the sine recursion and similarly for the cosine recursion. (Whether you take the simple recurrence $\theta_{n+2} = \theta_n + \theta_{n+1}$ and obscure it via sines, cosines, tangents, or, I suppose, other periodic functions, like Weierstrass's $\wp$-functions, you still get the same sequence of fractions of the functions's period(s), and the same period of the recurrence.)

$\endgroup$
  • $\begingroup$ +1, but the recurrence you display for the sine and cosine is in fact numerically unstable, and thus cannot be recommended in general. $\endgroup$ – J. M. will be back soon Dec 10 '15 at 5:18
  • $\begingroup$ @J.M. : I agree in part. If we're going to talk about stability, the initial tangent recursion isn't stable either. Replace $a_{n+1}$ by $a_{n+1}+\varepsilon$ and expand in powers of $\varepsilon$ and find that the error is magnified by $\frac{1+a_n^2}{(1-a_n a_{n+1})^2}$. This is harder to interpret than doing the same thing to the tangent version, yielding $1+\tan^2(a_n + a_{n+1})$, which is necessarily $\geq 1$. So the instability is baked into the problem and the alternate recursions are merely reflecting it. $\endgroup$ – Eric Towers Dec 10 '15 at 17:11
4
$\begingroup$

You can use NestWhileList to show the same pair of consecutive elements appear and then determine the period, e.g.

res = NestWhileList[{#[[2]], 
     Simplify[(#[[1]] + #[[2]])/(1 - #[[1]] #[[2]])]} &, {1, 
    1/Sqrt[3]}, Unequal, All];
plt = (Most@(First /@ res))~Join~(Style[#, Red, Bold] & /@ res[[-1]])
Length[plt] - 2
ListPlot[plt~Join~res[[1]], Joined -> True]

enter image description here

$\endgroup$
4
$\begingroup$

I add this as a separate answer to amplify another answer,in relation to $\tan(x+y)$.

This sequence could be pursued as Fibonacci or linear recurrence of angles,e.g.

mp[n_] := 
  Mod[MatrixPower[{{1, 1}, {1, 0}}, n].{Pi/6, Pi/4}, Pi, -Pi/2];
an[1] := mp[0][[2]];
an[2] := mp[0][[1]];
an[n_?(# > 2 &)] := First@mp[n - 2]
os[n_] := Tan[an[n]];
SequencePosition[an /@ Range[1, 26], {Pi/4, Pi/6}]
os /@ Range[26]

I exploited Fred Simons SequencePosition (yielding: {{1, 2}, {25, 26}}) and the last line reproduces the original sequence:

$\left\{1,\frac{1}{\sqrt{3}},\sqrt{3}+2,-\sqrt{3}-2,0,-\sqrt{3}-2,-\sqrt{3}-2,\frac{1}{\sqrt{3}},-1,\sqrt{3}-2,-\sqrt{3},-\sqrt{3}-2,1,-\frac{1}{\sqrt{3}},2-\sqrt{3},\sqrt{3}-2,0,\sqrt{3}-2,\sqrt{3}-2,-\frac{1}{\sqrt{3}},-1,-\sqrt{3}-2,\sqrt{3},\sqrt{3}-2,1,\frac{1}{\sqrt{3}}\right\}$

$\endgroup$
3
$\begingroup$

A quick visual inspection helps to see how a pattern with period 24 pops up

n = 100;   

tab = Table[fun[i], {i, 1, n}];

pos = Flatten @ Position[tab, First @ tab]

{1, 13, 25, 37, 49, 61, 73, 85, 97}

tra = Transpose[{pos, Array[First @ tab &, Length @ pos]}];

ListPlot[{tab, tra}, GridLines -> {pos[[;; ;; 2]], None}]

enter image description here

$\endgroup$
2
$\begingroup$

I rewrote the code to produce a list of the current "state" of the recursion, which is given by the last two terms $(a_n,a_{n-1}$. The system is autonomous, so the state is completely determined by this pair. Therefore if a pair ever occurs again, there is a loop.

Some utility functions:

step =           (* basic recursive step *)
  {Last[#], Simplify[Plus @@ #/(1 - Times @@ #)]} &;
steps[n_] :=     (* equivalent to fun[n], except it returns a list of all states *)
  NestList[step, {1, 1/Sqrt[3]}, n];
findloop[] :=    (* like steps[n], except it runs until a loop is found; caveat forever *)
  NestWhileList[step, {1, 1/Sqrt[3]}, UnsameQ, All];

First proof:

gr = Rule @@@ Partition[findloop[], 2, 1] // Graph;
cycle = First@FindHamiltonianCycle[gr];
Length@cycle
HighlightGraph[gr, PathGraph[cycle]]

Mathematica graphics

Less fussy proof, assuming you can make a good guess (i.e. 24):

Rule @@@ Partition[steps[24], 2, 1] // Graph
(* graph similar to above, but without highlighting *)

Even less fussy, assuming you can make a sufficient guess (i.e. 100):

gr = Rule @@@ Partition[steps[100], 2, 1] // Graph
Length@VertexList[gr]
(*
  graph similar to above, but with multiple edges
  24  <-- number of vertices == period
*)

Here is a translation of Brent's algorithm for finding the length of a cycle. See Brent, doi:10.1007/BF01933190 (1980). This could be compiled when step can be compiled for cases where a loop is long or occurs a long way down the sequence.

brentcyclelength[step_, x0_, ρ_: 2, u_: 0] := 
 Module[{x, y, r, j, k, done = False},
  y = x0;
  r = ρ^u;
  k = 0;
  While[! done,
   x = y;
   j = k;
   r = ρ*r;
   While[! done && k < r,
    k++;
    y = step[y];
    done = x == y
    ];
   ];
  k - j
  ]

brentcyclelength[step, {1, 1/Sqrt[3]}]
(*  24  *)
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.