# Asymptotic expansion of $u_{n + 1} = \frac12 \arctan(u_n)$

(I'm aware of Asymptotic expansion of $v_n = 2^nu_n$ where $u_{n+1} = \dfrac{1}{2}\arctan(u_n)$ but it has no answers…)

Let be $$u_0 \in \mathbb{R}$$ and the sequence $$(u_n)_n$$ defined by: $$u_{n + 1} = \frac12 \arctan(u_n)$$.

I define also: $$v_n = 2^n u_n$$, so I can show that: $$\lim (u_n)_n = 0$$ (by studying $$x \mapsto \frac12 \arctan(x)$$), thus, I can show that $$(v_n)_n$$ is monotone and converges because it is bound.

Now, I conclude: $$u_n \sim \dfrac{l}{2^n}$$, I'd like to determine $$l$$ more precisely.

Here is what I tried, I suspect $$l$$ to be something like $$f(\pi)$$ for some $$f$$ :

• push the asymptotic expansion of $$\arctan$$ to the 2nd order and reinject it ;
• use $$\arctan(u_n) + \arctan(1/u_n) = \dfrac{\pi}{2}$$ ;
• use series techniques to look for $$\sum v_{n + 1} - v_n$$, maybe conclude using Cesaro summation
• Have you tried to use a computer to obtain an approximate value of $\ell$ ? May 31 '19 at 19:13
• @Somos Not really for this problem. Note that the answer will depend on $u_0$. May 31 '19 at 19:32
• I can't find a closed form, but can find the first terms in a power-series $$f(\ell) = \ell +\frac{4 \ell^3}{9}+ \frac{176 \ell^5}{675}+ \frac{142144 \ell^7}{893025} + \frac{67031296\ell^9}{683164125} + \frac{777200229376 \ell^{11}}{12812743164375} + \ldots$$ which is such that $u_n = f(\ell/2^n)$ and $\ell$ is determined via $f(\ell) = u_0$. For small values of $u_0 \ll 1$ we have $\ell \approx u_0$. The general term in this power-series is $\sim (0.78\ell)^{2k+1}$. May 31 '19 at 22:04
• @JeanMarie Yes but it looks like a bit unstable, I tried to look at $\log (u_n)$ but was not able to get any fixed value, maybe I should use convergence acceleration techniques Jun 2 '19 at 3:21
• @Winther You should transform your comment into an answer. Jun 2 '19 at 7:36

The iteration has the form $$u_{n+1}=a_1u_n+a_3u_n^3+...$$ As usual in such situations (See the answer in Convergence of $\sqrt{n}x_{n}$ where $x_{n+1} = \sin(x_{n})$ with citation of de Bruijn: "Asymptotic Methods ..."), one can try a Bernoulli-like approach and examine the dynamics of $$u_n^{-2}$$. There one finds $$\frac1{u_{n+1}^2}=\frac4{u_n^2(1-\frac13u_n^2+\frac15u_n^4\mp...)^2} =\frac4{u_n^2}+\frac83-\frac4{15}u_n^2+O(u_n^4)\tag1$$ Thus for a first approximation use $$x_{n+1}=4x_n+\frac83\iff x_{n+1}+\frac89=4(x_n+\frac89)$$ so that $$u_n^{-2}\sim x_n=4^n(x_0+\frac89)-\frac89.\tag2$$

This gives as first approximation $$u_n\sim \frac{2^{-n}u_0}{\sqrt{1+\frac89u_0^2(1-4^{-n})}}.\tag3$$

For the next term use $$v_n=(u_n^{-2}+\frac89)^{-1}$$ and express (1) in terms of $$v_n$$ $$\frac1{v_{n+1}}=\frac4{v_n}-\frac4{15}\frac{v_n}{1-\frac89v_n}+O(v_n^2) \tag4$$ so that $$\frac1{v_{n+1}}-\frac{4}{15^2}v_{n+1}=\frac4{v_n}-\frac1{15} v_n - \frac{1}{15^2}v_n+O(v_n^2)=4\left(\frac1{v_{n}}-\frac{4}{15^2}v_{n}\right)+O(v_n^2) \tag5$$ and consequently $$\frac1{v_{n}}-\frac{4}{15^2}v_{n}=4^n\left(\frac1{v_{0}}-\frac{4}{15^2}v_{0}+O(v_0^2)\right) \tag6$$ As $$\frac1v-\frac{4}{15^2}v=\frac1v(1-\frac4{15^2}v^2)$$, leaving out the second term adds an error $$O(v_n^2)$$ which is a small fraction of $$O(v_0^2)$$. Thus $$\frac1{u_n^2}+\frac89=\frac1{v_n}=4^n\left(\frac1{v_{0}}-\frac{4}{15^2}v_{0}+O(v_0^2)\right)=4^n\left(\frac1{u_0^2}+\frac89-\frac{4}{15^2}\frac{u_0^2}{1+\frac89u_0^2}+O(u_0^4)\right)\tag7$$ so that the improved approximation is $$u_n=\frac{2^{-n}u_0}{\sqrt{1+\frac89u_0^2(1-4^{-n})-\frac{4}{25}\frac{u_0^4}{9+8u_0^2}+O(u_0^6)}} \tag8$$

• Is the result exact? Jul 12 '19 at 10:43
• No, in an exact formula there will be additional terms $O(8^{-n})$ and smaller, and possibly additional constants $O(u_0^2)$ under the square root. Note that the $x_n$ iteration only uses the first two terms of the previous formula. The summation of $4^ku_{n-k}^2$ is not arbitrarily small. I'm working on adding the next step. Jul 12 '19 at 10:52
• I have a feel that, as you work on, you will finally obtain an answer equivalent to mine. Looking forward to seeing a beautiful result, as mine is extremely ugly :( Jul 12 '19 at 10:57
• The equation right after the word ‘consequently’ should have an error term of $O(4^n v_0^2)$. Jul 12 '19 at 14:40
• Yes, obviously. The coefficients form the recursion are $1,4,...,4^{n-1}$ while the terms themselves go like $16^{1-n},16^{2-n},...,1$ which makes a sum $O(4^n)$. Jul 12 '19 at 15:10

(For easier discussion, I suggest you to read the introduction of Schroder's equation and the section on 'Conjugacy' of iterated function, in case you are not familiar with these topics.)

Let $$f(x)=\frac12\arctan x$$, and $$f_n(x)$$ be the $$n$$th iteration of $$f$$.

Let us reduce functional iteration to multiplication: if we can solve the corresponding Schroder's equation $$\Psi(f(x))=s\Psi(x)$$

then it is well known (and also straightforward) that $$f_n(x)=\Psi^{-1}(f'(a)^n\cdot\Psi(x))$$ where $$a$$ is a fixed point of $$f$$.

For the moment, let us focus on $$\Psi(f(x))=s\Psi(x)$$.

Clearly, in our case, $$a=0$$, and $$s=f'(a)=\frac12$$.

For $$a = 0$$, if $$h$$ is analytic on the unit disk, fixes $$0$$, and $$0 < |h′(0)| < 1$$, then Gabriel Koenigs showed in 1884 that there is an analytic (non-trivial) $$\Psi$$ satisfying Schröder's equation $$\Psi(h(x))=s\Psi(x)$$.

Thus, $$\Psi$$ is analytic.

A few more observations:

1. $$\Psi(0)=0$$.
2. $$\Psi'(0)$$ is up to our choice, since if a function $$\psi$$ is a solution to the Schroder's equation, then so is $$k\cdot \psi$$ for any constant $$k$$. For convenience, set $$\Psi'(0)=1$$.
3. All other Taylor series coefficients of $$\Psi$$ are then uniquely determined, and can be found recursively. (The method will be illustrated below.)
4. By Lagrange inversion theorem, $$\Psi$$ is invertible in a neighbourhood of $$0$$, and $$\Psi^{-1}(z)=0+\frac1{\Psi'(0)}z+o(z)\implies \Psi^{-1}(z)\sim z\quad(z\to 0)$$.
5. Therefore, $$f_n(x)=\Psi^{-1}(f'(a)^n\cdot\Psi(x))=\Psi^{-1}(2^{-n}\Psi(x))\sim 2^{-n}\Psi(x)$$ as $$n\to\infty$$.

Hence, for the limit the OP wanted to evaluate, $$\ell:=\lim_{n\to\infty}2^nf_n(x_0)=\Psi(x_0)$$

We shall now determine all the Taylor series coefficients of $$\Psi(x)$$ (valid only for $$|x|<1$$), since it can be assumed $$0\le x_0<1$$.

Obviously, $$\Psi$$ is an odd function. Let $$\Psi(x)=x+\sum^\infty_{k=2}\phi_{2k-1} x^{2k-1}$$

The basic idea is to repeatedly differentiate both sides of $$\Psi(f(x))=s\Psi(x)$$ and substitute in $$x=0$$, then recursively solve for the coefficients.

For example, differentiating both sides three times and substitute in $$x=0$$, we obtain $$-\Psi'(0)+\frac18\Psi'''(0)=\frac12\Psi'''(0)\implies\phi_3=-\frac49$$

Slightly modifying the notations of our respectable MSE user @Sangchul Lee, for $$\lambda=(\lambda_1,\lambda_2,\cdots,\lambda_n)$$ a $$n$$-tuple of non-negative integers:

• write $$\lambda \vdash n$$ if $$\sum^n_{i=1}(2i-1)\lambda_i=2n-1$$.
• write $$|\lambda| = \sum_{i=1}^{n} \lambda_i$$.
• define the tuple factorial as $$\lambda !=\frac{|\lambda|!}{\lambda_1!\cdot\lambda_2!\cdots\lambda_n !}$$.

I will state, without proof, Faà di Bruno's formula for odd inner function:

$$(\Psi\circ f)^{(2n-1)}=(2n-1)!\sum_{\lambda \vdash n}\lambda!\cdot\phi_{|\lambda|}\prod^n_{i=1}\left(\frac{f^{(2i-1)}(0)}{(2i-1)!}\right)^{\lambda_i}$$

$$\implies \frac12\phi_{2n-1}=\sum_{\lambda \vdash n}\lambda!\cdot\phi_{|\lambda|}\prod^n_{i=1}\left(\frac{(-1)^{i+1}}{2(2i-1)}\right)^{\lambda_i}$$

Further simplifications lead to the final result:

$$\ell=\Psi(x_0)=\sum^\infty_{k=1}\phi_{2k-1} x_0^{2k-1} \qquad{\text{where}}\qquad \phi_1=1$$

$$\phi_{2n-1}=\frac{(-1)^{n}}{2^{-1}-2^{1-2n}}\sum_{\substack{\lambda \vdash n \\ \lambda_1\ne 2n-1}}\phi_{|\lambda|}\frac{\lambda! (-1)^{(|\lambda|+1)/2}}{2^{|\lambda|}}\prod^n_{i=1}\frac1{(2i-1)^{\lambda_i}}$$

Yeah, I know it’s ugly. But that’s the best we can obtain.

• I have put a Magma script to solve the functional equation directly in a new answer. Jul 12 '19 at 15:57

For convenience we make a slight generalization of the problem. Let $$\,u_0\,$$ and $$\,y\,$$ be given numbers and suppose $$\,u_{n+1} = y \arctan(u_n)\,$$ for $$\,n\ge 0\,$$ where $$\ y=1/2\$$ in your original recursion. Define with power series the function $$F(x,y,z) := z\left(x + \frac{-1+z^2}{1-y^2}\frac{x^3}3 +\frac{(1-z^2)((3-2z^2)+y^2(2-3z^2)}{(1-y^2)(1-y^4)}\frac{x^5}{15} + O(x^7) \right)$$ which satisfies the equation $$\,F(x,y,yz) = \arctan(F(x,y,z))y.\,$$ Then we get the equation $$\, u_n = F(x,y,y^n)\,$$ where $$\, x = \lim_{n\to\infty} u_n/y^n.\,$$ I know some more terms in the power series expansion if you are interested. Thus we get the result $$\, u_n \approx y^n(x - (1-y^{2n})x^3/(3(1-y^2))).\,$$

Partial answer for $$u_0>0$$, then $$u_{n+1}-\frac{u_n}{2}=\frac{1}{2}(\arctan{u_n}-u_n)<0$$ because $$f(x)=\arctan{x}-x<0$$ for positive $$x$$, thus $$0 Using MVT, $$\exists z\in(u_{n+1},u_n)$$ s.t. $$u_{n+1}-u_n=\frac{1}{2}\left(\arctan{u_n}-\arctan{u_{n-1}}\right)= \frac{1}{2}\frac{u_{n}-u_{n-1}}{z^2+1}$$ or (because $$\color{red}{u_n-u_{n-1}<0}$$) $$\frac{1}{2}\cdot \frac{u_{n}-u_{n-1}}{u_{n+1}^2+1}< u_{n+1}-u_n< \frac{1}{2}\cdot \frac{u_{n}-u_{n-1}}{u_{n}^2+1}$$ or $$\frac{u_{0}-u_{1}}{2^n}\prod\limits_{k=1}^n\frac{1}{u_{k+1}^2+1}< u_{n+1}-u_n< \frac{u_{0}-u_{1}}{2^n}\prod\limits_{k=1}^n\frac{1}{u_{k}^2+1}$$ Considering $$u_{n+1}-u_n \sim -\frac{l}{2^{n+1}}$$ then $$\frac{u_{0}-u_{1}}{2^n}\prod\limits_{k=1}^n\frac{1}{u_{k+1}^2+1}> \frac{l}{2^{n+1}}> \frac{u_{0}-u_{1}}{2^n}\prod\limits_{k=1}^n\frac{1}{u_{k}^2+1}$$ or $$\frac{2(u_{0}-u_{1})}{\prod\limits_{k=1}^n\left(u_{k+1}^2+1\right)}= \frac{2(u_{0}-u_{1})\left(u_{1}^2+1\right)}{\prod\limits_{k=1}^{n+1}\left(u_{k}^2+1\right)}> l> \frac{2(u_{0}-u_{1})}{\prod\limits_{k=1}^n\left(u_{k}^2+1\right)}$$ or $$L_1>l>L_2$$ where $$L_2=\frac{2(u_{0}-u_{1})}{\lim\limits_{n\to\infty}\prod\limits_{k=1}^n(u_{k+1}^2+1)} \text{ and } L_1=L_2\left(u_{1}^2+1\right) \tag{2}$$

So, it looks like Robert (see the comments) was right, it depends on $$u_0$$.

Note: the following limit exists $$\lim\limits_{n\to\infty}\prod\limits_{k=1}^n(u_{k+1}^2+1)$$ because $$0<\sum\limits_{k=1}\ln(u_{k+1}^2+1)<\sum\limits_{k=1}u_{k+1}^2<\infty$$ by ratio test from $$(1)$$.

The following code is computing $$(2)$$ but with a $$\frac{1}{u_0}$$ factor. You will notice a certain stability for $$\frac{L_1}{2^n u_0 \cdot u_n}$$ and $$\frac{L_2}{2^n u_0 \cdot u_n}$$ for various $$u_0$$

from math import atan
from math import pow

N = 300
U_0 = 190.0

u = []

it = U_0
u.append(it)

for i in range(1, N):
it = 0.5 * atan(it)
u.append(it)

val = 1.0
for i in range(1, N):
val *= (u[i] * u[i] + 1.0)

L2 = (2.0 * (u - u) / val) / u
L1 = L2 * (u * u + 1.0)
MID = (L1 + L2) / 2.0

print "limit L1 =",L1
print "limit L2 =",L2
print "limit MID =",MID

for i in range(N-100, N):
Lp1 = L1 / pow(2, i)
Lp2 = L2 / pow(2, i)
MIDp = MID / pow(2, i)

r1 = Lp1 / u[i]
r2 = Lp2 / u[i]
rMID = MIDp / u[i]

print Lp2," vs ",u[i]," vs ",Lp1," --- ",MIDp
print r2," vs ",r1," --- ",rMID


Try it here.

## Complement to the answer of @Szeto

In many cases when you start tinkering with the Faà di Bruno formula, you will be better served computing with truncated Taylor series.

So we want to solve $$Ψ(x)=2Ψ(\tfrac12\arctan(x))$$ where $$Ψ(x)\sim x$$ for $$x\approx 0$$ by the scaling normalization. As $$\arctan(x)\sim x$$ for $$x\approx 0$$, the coefficient determination for $$Ψ(x)=x+c_2x^2+c_3x^3+...$$ is a finite problem for each coefficient, it is only influenced by lower degree coefficients. Thus assuming that the coefficients $$c_0=0,c_1=1,c_2,..c_{k-1}$$ are already determined, one gets the next coefficient from $$(1-2^{1-k})c_kx^k=A_k(x)-x+c_2(2^{-1}A_{k-1}(x)^2-x^2)+c_3(2^{-2}A_{k-2}(x)^3-x^3)+...+c_{k-1}(2^{2-k}A_{2}(x)^{k-1}-x^{k-1})$$ by comparing the coefficients of $$x^k$$ on both sides. The $$A_k(x)$$ are the $$k$$-th partial sums of the arcus tangent series, $$\arctan(x)=A_k(x)+O(x^{k+1})$$. This can be simplified, it is not necessary to subtract the lower powers, one can also take the odd nature of the series into account.

Using a CAS like Magma (online calculator one can extract the equation for the next coefficient directly from the unmodified equation with the following script:

A<a>:=FunctionField(Rationals());
PS<x>:=PowerSeriesRing(A);
Pol<z>:=PolynomialRing(Rationals());

Psi := x;
for k in [2..20] do
Psia := Psi+(a+O(x))*x^k;
eqn := Coefficient( Psia-2*Evaluate(Psia, 1/2*Arctan(x+O(x^(k+1))) ), k );
c := Roots(Pol!eqn)[1,1]; k,c;
Psi +:= c*x^k;
end for;


which when executed gives in the end for $$\Psi(x)+O(x^{21})$$

x - 4/9*x^3 + 224/675*x^5 - 51008/178605*x^7 + 25619968/97594875*x^9
- 91726170112/366078376125*x^11 + 45580629370863616/186023558824228125*x^13
- 171377650156414910464/703297837896306778125*x^15
+ 56540215172481124229054464/230453119032672323522109375*x^17
- 353563937806248194123298285027328/1417897477708832149477498284609375*x^19


The inverse function $$\Psi^{-1}$$ is then obtained as

x + 4/9*x^3+176/675*x^5 + 142144/893025*x^7 + 67031296/683164125*x^9
+ 777200229376/12812743164375*x^11
+ 76806067707244544/2046259147066509375*x^13
+ 7434789485314586820608/320000516242819584046875*x^15
+ 3317928226689969972317978624/230683572151704995845631484375*x^17
+ 30692357195871908183846360294096896/3446908768310170955379798329885390625*x^19


which is the series in the comment to the question by @Winther