Solve the recurrence relation $a_{n+1}=\frac{2a_n^2}{1-2a_n^2}$ How do we solve the recurrence relation $a_{n+1}=\frac{2a_n^2}{1-2a_n^2}$?
I found this problem on page 56 of Carl M. Bender and Steven A. Orszag's book Advanced Mathematical Methods for Scientists and Engineers, in the section that deals with nonlinear difference equations.
I tried using the substitution $b_n=\frac{1}{2a_n}$ to get $b_{n+1}=b_n^2-\frac12$. I did some searching here and found that there does not seem to be a (known) closed-form solution to this quadratic map.
In this section of the book, the author only mentioned a few examples of solving nonlinear difference equations using substitutions and generating functions, and does not mention anything about quadratic maps, so I assume this problem can be solved by appropriate substitutions.
Could you give me some hints? The solution can involve summation, etc., if necessary.
 A: This answer to your question may not have a closed form, but it is not that complicated. Define the real constants
$$ w :=\! \sqrt{3}, \;
p :=\! (-1\!+\!w)/2, \;
q :=\! (-1\!-\!w)/2. \tag{1} $$
Define the recurrence function
$$ f(x) := 2x^2/(1-2x^2). \tag{2} $$
Notice that $\,p\,$ is a repelling while $\,q\,$ is
an attracting fixed point of $\,f(x).\,$
Define the convergent power series
$$ s(x) \!:=\! x \!+\!\! \left(\frac12\!-\!\frac56w\right)\!x^2
 \!+\!\! \left(\frac73\!-\!\frac23w\right)\!x^3 \!\!+\! O(x^4).\! 
\tag{3} $$
Then the following functional equation holds
$$ f(q+s(x)) = q+s(-2px) \tag{4} $$
where the the coefficients of $\,s(x)\,$ are
uniquely determined so that the functional
equation holds.
Thus the convergent ($\;|\!-\!2p|<1\;$) sequence 
for any $\,x_0\,$ defined by
$$ a_n := q + s((-2p)^n x_0) \tag{5} $$
satisfies $\,a_n \to q\,$ as $\,n\to\infty\,$
and the recursion
$$ a_{n+1} = f(a_n). \tag{6} $$
For your information here is a PARI/GP
code to calculate the power series for $\,s(x)\,$
w = quadgen(12); p = (-1+w)/2; q = (-1-w)/2;
f(x) = 2*x^2/(1-2*x^2); 
nxt(n) = {my(s, c='c);
  s = truncate(stx + O(x)^n) + c*x^n*(1 + O(x));
  s = truncate(f(q+s) - (q+subst(s, x, -2*p*x)))/x^n;
  sx -= x^n * polcoeff(s, 0, c)/polcoeff(s, 1, c)};
sx = x ; for(n=2, M=9, nxt(n)); print(sx + x*O(x)^M);

