A more natural solution to finding the general terms of a recurrence relation in $2$ variables A high school contest math problem in a problem book:

Find the general terms of
$$a_{1}=a,\quad b_{1}=b,\quad a_{ n + 1 }=\frac { 2 a _ { n } b _ { n } } { a _ { n } + b _ { n } },\quad b_{ n + 1 }=\sqrt{ a_{n+1} b _ {n}}$$

The easiest method is to use the following trigonometric identities:

$$2 \tan \frac { \theta } { 2 } = \frac { 2 \tan \theta \sin \theta } { \tan \theta + \sin \theta } , \quad 2 \sin \frac { \theta } { 2 } = \sqrt { 2 \tan \frac { \theta } { 2 } \sin \theta }$$

We make the substitution

$$a _ { 1 } = a = c \cdot \tan \theta , \quad b _ { 1 } = b = c \cdot \sin \theta$$

Observe that

$$a_{2}=\frac{2a_1b_1}{a_1+b_1}=\frac{2c\tan\theta\sin\theta}{\tan\theta+\sin\theta}=2c\tan\frac{\theta}{2}$$
$$a_{3}=\frac{2a_2b_2}{a_2+b_2}=\frac{2(2c)\tan\frac{\theta}{2}\sin\frac{\theta}{2}}{\tan\frac{\theta}{2}+\sin\frac{\theta}{2}}=2^2c\tan\frac{\theta}{2^2}$$
$$...$$
$$a_n= 2^{n-1} c \tan \frac { \theta } { 2 ^ { n - 1 } }$$

It's possible to express the above formula in terms of $a$ and $b$, just notice

$$\theta=\cos^{-1}\left(\frac{b}{a}\right), c=\frac{a}{\tan\theta}$$

Similarly,

$$b_n=2 ^ { n - 1 } c\sin \frac { \theta } { 2 ^ { n - 1 } }$$

Although I consider the above solution as elegant, it's unnatural at first glance, as the $2$ trigonometric identities are rarely used in practice, at least for those who are mediocre in trigonometry. So my question is:

Is there an alternative solution (non-trigonometric) to the above recurrence relation that is more natural in a certain sense?

Edit: The above method has a defect: complex numbers may pop up for some real values of $a$ and $b$ in the substitution step.
There are several possible methods. A general method that I have
used for similar compound means is the following. First, we want
to find a recurrence for each sequence by itself. Thus, given
$\,a_0 = a, \,  b_0 = b, \, b_1 = \sqrt{a_1b_0}\,$ 
we solve for $\,a_1\,$ getting
$\,a_1 = b_1^2/b_0.\,$ Similarly $\,a_2 = b_2^2/b_1.\,$ By using
$\,a_2 = 2a_1b_1/(a_1+b_1)\,$
 we get a quadratic in $b_2$ given $b_0$ and $b_1$
 $$ b_2^2(b_0+b_1)-2 b_1^2 = 0. \tag{1} $$
Now we need an Ansatz. A little bit of numerical work suggests that the convergence of $b_n$ is linear. Thus, assume for some $\,0&lt;r&lt;1,\,$
without loss of generality, that $\,b_n = f(x r^n)\,$ for $\,n\ge 0\,$
where $$\,f(x) := 1 + x + c_2 x^2 + c_3 x^3 + c_4 x^4 + O(x^5). \tag{2}$$
Substitute this in equation $(1)$. The left side is
$\,(r-1)(4r-1)x + O(x^2).\,$ Solving for $\,r\,$ by the hypothesis,
$\,r=1/4.\,$ Solving for the other coefficients of $\,f(x)\,$ using equation
$(1)$ we get
$$ f(x) = 1 + x + 3/10 x^2 + 3/70 x^3 + 1/280 x^4 + 3/15400 x^5 + O(x)^6.
\tag{3}$$
Looking for the denominators leads to OEIS sequence A007019 


  Denominators of coefficients of the Taylor series of sinh(sqrt(2*x))/(sqrt(2*x))


This gives us the solution
 $$ f(x) = \sinh(\sqrt{6x})/\sqrt{6x}). \tag{4}$$
Combining our results gives us
$$ b_n = c\sinh(t/2^n)/(t/2^n) \tag{5}$$

Now we do something similar for $\,a_n.\,$ Thus, given
$\,a_1 = 2 a_0 b_0/(a_0+b_0)\,$ we solve for $\,b_0\,$ getting
$\,b_0 = a_0 a_1/(2a_0 - a_1).\,$ Similarly,
$\,b_1 = a_1 a_2/(2a_1 - a_2)\,$ but since $\,b_1^2 = a_1 b_0\,$
we finally get
 $$ a_0(a_2^2 -4a_1a_2-4a_1^2)-a_1a_2^2 = 0. \tag{6}$$
This is a quadratic in $a_2$ given $a_0$ and $a_1$.

A little bit of numerical work suggests that the convergence of
 $a_n$ is linear. Thus, assume for some $\,0&lt;r&lt;1,\,$
without loss of generality, that $\,a_n = f(x r^n)\,$ for $\,n\ge 0\,$
where $$\,f(x) := 1 + x + c_2 x^2 + c_3 x^3 + c_4 x^4 + O(x^5). \tag{7}$$
 Substitute this in equation $(6)$. The left side is
$\, -(1-r)(1-4r)x + O(x^2).\,$ Solving for $\,r\,$ by the hypothesis,
$\,r=1/4.\,$ Solving for the other coefficients of $\,f(x)\,$ using equation
$(1)$ we get
$$ f(x) = 1 + x + 6/5 x^2 + 51/35 x^3 + 62/35 x^4 + O(x^5).\tag{8}$$
The numerators are all divisible by $3$ so this suggests looking at
$$ f(x/3) = 1 + 1/3 x + 2/15 x^2 + 17/315 x^3 + O(x^4). \tag{9}$$
A lookup of the numerators leads to OEIS squence A002430


  Numerators in Taylor series for tan(x). Also from Taylor series for tanh(x).


which gives $$ f(x/3) = \tanh(\sqrt{x})/\sqrt{x}. \tag{10}$$
Combining our results gives us
$$ a_n = c\tanh(t/2^n)/(t/2^n). \tag{11}$$

To summarize, both $\,c := \lim a_n = \lim b_n\,$ as $\,n\to\infty\,$
and $\,t\,$ depend on the initial values $\,a,b.\,$
Let $\,d := \sqrt{b^2-a^2},\,$ and $\, t := \log((b+d)/a). \,$ Then
$\, c = t\ a\ b/d.\,$

Note that if $\,0&lt;a&lt;b\,$ then $\,t&gt;0\,$ is real while if $\,0&lt;b&lt;a\,$ then $\,t\,$ is pure imaginary. Here are numerical examples for both cases: 

c=1.520691, d=1.732050, t=1.316957
 n    a_n       b_n       a_n/c    b_n/c    (t/2^n)^2
 0  1.000000  2.000000  0.657595  1.315191  1.734378
 1  1.333333  1.632993  0.876793  1.073849  0.433594
 2  1.468027  1.548315  0.965368  1.018165  0.108398
 3  1.507103  1.527570  0.991063  1.004523  0.027099
..  ........  ........  ........  ........  ........
10  1.520691  1.520692  0.999999  1.000000  0.000001

c=1.209199, d=-1.732050, t=1.047197 i
 n    a_n       b_n       a_n/c    b_n/c    (t/2^n)^2
 0  2.000000  1.000000  1.653987  0.826993  -1.096623
 1  1.333333  1.154701  1.102658  0.954929  -0.274155
 2  1.237604  1.195434  1.023459  0.988615  -0.068538
 3  1.216154  1.205759  1.005751  0.997146  -0.017134
..  ........  ........  ........  ........  .........
10  1.209200  1.209199  1.000000  0.999999  -0.000001
An algorithmic solution would be a “natural solution” in the sense that the algorithm performs the computations specified by the recurrence relations. For an example of a non-recursive model for the computation please check out the following: 

"Computing the N-th Term of a Recursive Relation Using a Non-Recursive Function – A Reply to a Question at Mathematics Stack Exchange"