I have $$s_1 = 1, s_n = (1-\frac{1}{4n^2})s_{n-1}.$$
I see that its limit exists, but cannot figure out what the limit is. How do I find its limit?
The only way I used (know) is: $$\lim s_n = \lim s_{n+1}$$ But here, I am just getting $L = L$.
Also, when I program this sequence on computer it's approaching to $0.636619...$ and Wolfram Alpha is giving me complicated recurrence equation solution.
Assume that $s_0=1$, then $$ s_n=\frac{((2n)!)^2\,(2n+1)}{2^{4n}\,(n!)^4}. $$ Stirling equivalent $$ k!\sim\sqrt{2\pi k}\,k^k\,\mathrm e^{-k}, $$ then yields $$ \lim\limits_{n\to\infty}s_n=\frac2\pi\approx0.63662. $$ If $s_1=1$ instead of $s_1=\frac34$, then $$\lim\limits_{n\to\infty}s_n=\frac43\cdot\frac2\pi.$$
There is a harmonic sum hiding here which we now examine in more detail.
Introduce $$P = \prod_{n\ge 1} \left(1-\frac{1}{4n^2}\right)$$ so that $$ S= \log P = \sum_{n\ge 1} \log\left(1-\frac{1}{4n^2}\right)$$
The sum term may be evaluated by inverting its Mellin transform. Put $$S(x) = \sum_{n\ge 1} \log\left(1-\frac{1}{4(xn)^2}\right)$$ so that we are interested in $S(1).$
Recall the harmonic sum identity $$\mathfrak{M}\left(\sum_{k\ge 1} \lambda_k g(\mu_k x);s\right) = \left(\sum_{k\ge 1} \frac{\lambda_k}{\mu_k^s} \right) g^*(s)$$ where $g^*(s)$ is the Mellin transform of $g(x).$
In the present case we have $$\lambda_k = 1, \quad \mu_k = k \quad \text{and} \quad g(x) = \log\left(1-\frac{1}{4x^2}\right).$$
We need the Mellin transform $g^*(s)$ of $g(x)$ which is $$\int_0^\infty \log\left(1-\frac{1}{4x^2}\right) x^{s-1} dx \\ = \left[\log\left(1-\frac{1}{4x^2}\right) \frac{x^s}{s}\right]_0^\infty + 2 \int_0^\infty \frac{1}{x(1-4 x^2)} \frac{x^s}{s} dx \\ = \frac{2}{s} \int_0^\infty \frac{1}{1-4 x^2} x^{s-1} dx.$$
The fundamental strip of the term from the integration by parts is $\langle0, 2\rangle,$ which is also the fundamental strip of the second Mellin transform, call it $h^*(s)$ of $h(x).$
This transform has poles at $\pm 1/2$ i.e. two poles on the real line. We will be using a semicircular contour of radius $R$ in the upper half plane with $R$ going to infinity to evaluate this transform. We put small semicircular indentations around the two poles. These are both counterclockwise and pick up half the residue at that point. The integral on the negative real line is $$\int_{-\infty}^0 h(x) x^{s-1} dx = - \int_\infty^0 h(-x) (-x)^{s-1} dx \\= e^{\pi i (s-1)} \int_0^\infty h(x) x^{s-1} dx = - e^{\pi i s} h^*(s).$$
This yields for $h^*(s)$ that $$h^*(s) \left(1-e^{\pi i s}\right) \\= \frac{1}{2}\times 2\pi i \left(\mathrm{Res}\left(\frac{1}{1-4 x^2} x^{s-1}; x=-1/2\right) +\mathrm{Res}\left(\frac{1}{1-4 x^2} x^{s-1}; x=1/2\right)\right) \\=\pi i \left(\frac{1}{4} (-1/2)^{s-1} - \frac{1}{4} (1/2)^{s-1}\right).$$ This gives for $h^*(s)$ that $$h^*(s) = \pi i \frac{1}{4} (1/2)^{s-1} \frac{e^{i\pi(s-1)}-1}{1-e^{\pi i s}} = -\pi i \frac{1}{4} (1/2)^{s-1} \frac{1+e^{i\pi s}}{1-e^{\pi i s}} \\ = -\pi (1/2)^{s+1} \frac{i(e^{-i\pi s/2}+e^{i\pi s/2})}{e^{-i\pi s/2}-e^{\pi i s/2}} = \pi (1/2)^{s+1} \cot(\pi s/2).$$
It follows that the Mellin transform $Q(s)$ of the harmonic sum $S(x)$ is given by
$$Q(s) = \frac{2}{s} \pi (1/2)^{s+1} \cot(\pi s/2) \zeta(s)$$ which is $$ \frac{\pi}{s}\frac{\cot(\pi s/2)}{2^s} \zeta(s) \quad\text{because}\quad \sum_{k\ge 1} \frac{\lambda_k}{\mu_k^s} = \sum_{k\ge 1} \frac{1}{k^s} = \zeta(s)$$ for $\Re(s) > 1.$
The Mellin inversion integral here is determined by the abscissa of convergence of the zeta function term and is $$\frac{1}{2\pi i} \int_{3/2-i\infty}^{3/2+i\infty} Q(s)/x^s ds$$ which we evaluate by shifting it to the right for an expansion about infinity. (In any case the pole of the zeta function is cancled by a zero of the cotangent term.)
Collecting the contributions from the poles at the even positive integers of the contangent term we get the sum (with a minus due to the shift to the right and setting $x=1$ for the value $S(1)$)
$$-\sum_{q\ge 1} \frac{\pi}{2q} \frac{2}{\pi} \frac{\zeta(2q)}{2^{2q}} = - \sum_{q\ge 1} \frac{\zeta(2q)}{q 2^{2q}}.$$ In terms of Bernoulli numbers this becomes $$ \sum_{q\ge 1} \frac{(-1)^q B_{2q} (2\pi)^{2q}}{2(2q)!} \frac{1}{q 2^{2q}} = \sum_{q\ge 1} \frac{i^{2q} B_{2q} \pi^{2q}}{2q(2q)!}.$$ This is $$ -\frac{i B_1\pi}{1\times 1!} + \sum_{q\ge 1} \frac{i^q B_q \pi^q}{q\times q!} \quad\text{or}\quad \frac{1}{2} i\pi + \sum_{q\ge 1} \frac{i^q B_q \pi^q}{q\times q!}.$$
Now recall that $$-1 + \frac{t}{e^t-1} = \sum_{q\ge 1} B_q \frac{t^q}{q!} \quad\text{so that}\quad -\frac{1}{t} + \frac{1}{e^t-1} = \sum_{q\ge 1} B_q \frac{t^{q-1}}{q!}.$$ Integration of the generating function yields $$-\log t + \log(e^t-1) - t = \sum_{q\ge 1} B_q \frac{t^q}{q\times q!}.$$
The left and the right are zero in the limit as $t$ goes to zero so there are no problems with the constant that appeared during the integration.
Putting $t=i\pi$ and collecting everything we finally obtain $$S(1) = \frac{1}{2}\pi i - \log(i\pi) + \log(-2) - i \pi \\= \frac{1}{2}\pi i - \left(\frac{1}{2} \pi i + \log \pi\right) + \left(\log 2 + i\pi\right) - i\pi = \log\left(\frac{2}{\pi}\right)$$ and therefore $$P = \exp\log\left(\frac{2}{\pi}\right) = \frac{2}{\pi}.$$
Observation. The reader may well contend that in our last step we have tacitly chosen a branch cut for the logarithm. The fact is however that the result does not depend on the sheet. Suppose we have a logarithm that produces arguments in the range $[2\pi, 4\pi).$ We get
$$ \frac{1}{2}\pi i - \left(\frac{5}{2} \pi i + \log \pi\right) + \left(\log 2 + 3i\pi\right) - i\pi = \log\left(\frac{2}{\pi}\right).$$
For a logarithm with the cut on the positive imaginary axis and argument in $[\pi/2, 5\pi/2)$ we get
$$ \frac{1}{2}\pi i - \left(\frac{1}{2} \pi i + \log \pi\right) + \left(\log 2 + i\pi\right) - i\pi = \log\left(\frac{2}{\pi}\right).$$
There is another "divergent" Mellin transform with poles on the positive real line at this MSE link.
Clearly, the limit of the sequence is equal to $\displaystyle \prod_{n = 2}^{\infty} \left(\frac{4n^2 - 1}{4n^2} \right)$. This is the reciprocal of the Wallis product (without the first term). The Wallis Product is equal to $\frac{\pi}{2}$ with the first term, $\frac{4}{3}$, so without it it is equal to $\frac{3 \pi}{8}$. Since this is the reciprocal, the desired answer is $\frac{8}{3 \pi}$. (By your program approximation, I suspect you meant to say that $s_0 = 1$ and not $s_1 = 1$, in which case the answer would just be $\frac{2}{\pi}$
One may prove the Wallis Product in many ways, originally by expressing $\frac{\sin x}{x}$ as a polynomial with roots at $\pm n \pi$ for $n \in \mathbb{Z}^+$. You may read about it here: http://en.wikipedia.org/wiki/Wallis_product
$\newcommand{\+}{^{\dagger}} \newcommand{\angles}[1]{\left\langle\, #1 \,\right\rangle} \newcommand{\braces}[1]{\left\lbrace\, #1 \,\right\rbrace} \newcommand{\bracks}[1]{\left\lbrack\, #1 \,\right\rbrack} \newcommand{\ceil}[1]{\,\left\lceil\, #1 \,\right\rceil\,} \newcommand{\dd}{{\rm d}} \newcommand{\down}{\downarrow} \newcommand{\ds}[1]{\displaystyle{#1}} \newcommand{\expo}[1]{\,{\rm e}^{#1}\,} \newcommand{\fermi}{\,{\rm f}} \newcommand{\floor}[1]{\,\left\lfloor #1 \right\rfloor\,} \newcommand{\half}{{1 \over 2}} \newcommand{\ic}{{\rm i}} \newcommand{\iff}{\Longleftrightarrow} \newcommand{\imp}{\Longrightarrow} \newcommand{\isdiv}{\,\left.\right\vert\,} \newcommand{\ket}[1]{\left\vert #1\right\rangle} \newcommand{\ol}[1]{\overline{#1}} \newcommand{\pars}[1]{\left(\, #1 \,\right)} \newcommand{\partiald}[3][]{\frac{\partial^{#1} #2}{\partial #3^{#1}}} \newcommand{\pp}{{\cal P}} \newcommand{\root}[2][]{\,\sqrt[#1]{\vphantom{\large A}\,#2\,}\,} \newcommand{\sech}{\,{\rm sech}} \newcommand{\sgn}{\,{\rm sgn}} \newcommand{\totald}[3][]{\frac{{\rm d}^{#1} #2}{{\rm d} #3^{#1}}} \newcommand{\ul}[1]{\underline{#1}} \newcommand{\verts}[1]{\left\vert\, #1 \,\right\vert} \newcommand{\wt}[1]{\widetilde{#1}}$ $\ds{s_{1} = 1\,,\qquad s_{n} = \pars{1 - {1 \over 4n^{2}}}s_{n - 1}\,,\quad n \geq 2 .\qquad\lim_{n \to \infty}s_{n}: {\large ?}}$
$$ s_{n} =\pars{1 - {1 \over 4n^{2}}}\bracks{1 - {1 \over 4\pars{n - 1}^{2}}} \ldots \pars{1 - {1 \over 4\times 2^{2}}}s_{1} $$
\begin{align} \ln\pars{s_{n}}&=\overbrace{\ln\pars{s_{1}}}^{\ds{0}}\ +\ \sum_{k = 2}^{n}\ln\pars{1 - {1 \over 4k^{2}}} =-\sum_{k = 2}^{n}\int_{0}^{1}{\dd x \over x + 4k^{2} - 1} \\[3mm]&=-\,{1 \over 4}\int_{0}^{1}\sum_{k = 0}^{n - 2} {1 \over \pars{k + 2 + \root{1 - x}/2}\pars{k + 2 - \root{1 - x}/2}} \end{align}
With $\ds{\quad t \equiv \root{1 - x}\quad\imp\quad x = 1 - t^{2}}$: \begin{align} \color{#c00000}{\lim_{n \to \infty}\ln\pars{s_{n}}}& =-\,{1 \over 4} \int_{0}^{1}{\Psi\pars{2 + \root{1 - x}/2} - \Psi\pars{2 - \root{1 - x}/2} \over \root{1 - x}}\,\dd x \\[3mm]&=-\,\half\int_{0}^{1}\bracks{% \Psi\pars{2 + {t \over 2}} - \Psi\pars{2 - {t \over 2}}}\,\dd t =-\left.\ln\pars{\Gamma\pars{2 + {t \over 2}}\Gamma\pars{2 - {t \over 2}}} \vphantom{\LARGE A}\right\vert_{0}^{1} \\[3mm]&=-\ln\pars{\Gamma\pars{5 \over 2}\Gamma\pars{3 \over 2}} +\ln\pars{\Gamma\pars{2}\Gamma\pars{2}} =-\ln\pars{{3 \over 2}\,\Gamma^{2}\pars{3 \over 2}} \\[3mm]&=-\ln\pars{{3 \over 2}\,{1 \over 4}\,\Gamma^{2}\pars{\half}} =-\ln\pars{3\pi \over 8}=\color{#c00000}{\ln\pars{8 \over 3\pi}} \end{align}
Then, $$ \color{#00f}{\large\lim_{n \to \infty}s_{n} = {8 \over 3\pi}} \approx 0.8488 $$
$\ds{\Gamma\pars{z}}$ and $\ds{\Psi\pars{z}}$ are the Gamma and Digamma Functions, respectively.
