How was this approximation of $\pi$ involving $\sqrt{5}$ arrived at?

Ramanujan discussed this and many similar approximations in his monumental paper Modular Equations and Approximations to $\pi$. Let me describe his idea in brief which is based on deep interconnection between the theory of elliptic integrals and theta functions.

Let $0 < k < 1$ and $k' = \sqrt{1 - k^{2}}$ then we define two elliptic integrals $$K(k) = \int_{0}^{\pi/2}\frac{dx}{\sqrt{1 - k^{2}\sin^{2}x}},\,E(k) = \int_{0}^{\pi/2}\sqrt{1 - k^{2}\sin^{2}x}\,dx\tag{1}$$ When the number $k$ (called modulus) is evident from context we use the symbols $K, E, K', E'$ for $K(k), E(k), K(k'), E(k')$. These numbers have a curious relation to $\pi$ which is not difficult to prove $$KE' + K'E - KK' = \frac{\pi}{2}\tag{2}$$ which goes by the name Legendre's Identity. This is the identity which forms the link between $\pi$ and elliptic integrals and it was exploited by Ramanujan to the fullest extent.

The formulas $(1)$ can be inverted in the sense that if we are given the values of $K, K'$ then we can obtain the corresponding value of $k$. But this happens in almost magical way via theta functions. Let us define a new variable $q$ called nome by $q = e^{-\pi K'/K}$ and we have the following theta functions (in simplified form) of Jacobi: \begin{align} \vartheta_{2}(q) &= \sum_{n = -\infty}^{\infty}q^{(n + 1/2)^{2}} = 2(q^{1/4} + q^{9/4} + q^{25/4} + \cdots )\notag\\ &= 2q^{1/4}\prod_{n = 1}^{\infty}(1 - q^{2n})(1 + q^{2n})^{2}\tag{3a}\\ \vartheta_{3}(q) &= \sum_{n = -\infty}^{\infty}q^{n^{2}} = 1 + 2q + 2q^{4} + 2q^{9} + \cdots\notag\\ &= \prod_{n = 1}^{\infty}(1 - q^{2n})(1 + q^{2n - 1})^{2}\tag{3b}\\ \vartheta_{4}(q) &= \sum_{n = -\infty}^{\infty}(-1)^{n}q^{n^{2}} = 1 - 2q + 2q^{4} - 2q^{9} + \cdots\notag\\ &= \prod_{n = 1}^{\infty}(1 - q^{2n})(1 - q^{2n - 1})^{2}\tag{3c} \end{align} and the values $k, k', K$ are obtained via formulas: $$k = \frac{\vartheta_{2}^{2}(q)}{\vartheta_{3}^{2}(q)},\, k' = \frac{\vartheta_{4}^{2}(q)}{\vartheta_{3}^{2}(q)},\, \frac{2K}{\pi} = \vartheta_{3}^{2}(q)\tag{4}$$

It is interesting to study the dependence of the nome $q = e^{-\pi K'/K}$ on modulus $k$. Equivalently one can also analyze the behavior of the ratio $K'/K$ as $k$ varies in the interval $(0, 1)$. It can be easily proved using the definitions $(1)$ that $K'/K$ is strictly decreasing as $k$ increases and $K'/K \to \infty$ as $k \to 0^{+} $ and $K'/K \to 0$ as $k \to 1^{-} $. Thus the ratio $K'/K$ as a function of $k$ establishes a bijection between interval $(0, 1)$ and the set of positive real numbers.

Next idea is to start with a specific value of $k$ and a positive number $n$ and then the number $nK'/K$ is also positive and hence (by the observation in previous paragraph) there is a unique number $l \in (0, 1)$ such that $nK'/K = K'(l)/K(l)$. It is customary to denote $K(l), K'(l)$ by $L, L'$ respectively and thus starting with a modulus $k \in (0, 1)$ and a positive real number $n$ we have found a unique $l \in (0, 1)$ such that $$\frac{L'}{L} = n\frac{K'}{K}$$ Jacobi extensively studied the integral transformations related to the elliptic integral $K(k)$ and proved that:

Jacobi's Theorem: If $n$ is a positive rational number and $k\in (0, 1)$ there there is a unique $l \in (0, 1)$ such that $$\frac{L'}{L} = n\frac{K'}{K}$$ and the relation between $k, l$ is algebraic. In other words in such a case there exists a polyomial $P(x, y)$ in two variables with integer coefficients such that $P(k, l) = 0$.

The algebraic relation between $k, l$ induced by the equation $$\frac{L'}{L} = n\frac{K'}{K}$$ is called a modular equation of degree $n$. And Jacobi gave such modular equations for $n = 3, 5$. Finding modular equations using integral transformations as described by Jacobi is really tough and impractical beyond $n = 7$.

Ramanujan was in love with these modular equations and he gave many modular equations for different positive integer values of $n$. But he had a much deeper insight into these topics and he asked:

What would happen if we add a further constraint $l=k'=\sqrt{1-k^{2}}$ in the modular equation connecting $k, l$?

With this additional constraint we now have an algebraic equation $P(k, k') = 0$ and thus both $k, k'$ are algebraic numbers. Also since $l = k'$ we have $L = K', L' = K$ so that $L'/L = nK'/K$ leads us to $L'/L = \sqrt{n}, K'/K = 1/\sqrt{n}$. The corresponding nomes for $k, l$ are now seen to be $e^{-\pi/\sqrt{n}}$ and $e^{-\pi\sqrt{n}}$. And thus we get the following deep theorem:

Theorem: If $n$ is a positive rational number and $q = e^{-\pi\sqrt{n}}$ then the numbers $k, k'$ given by equation $(4)$ in terms of $q$ are algebraic.

Ramanujan's great ability was manipulation of radicals and if there was an algebraic number of any mathematical significance then Ramanujan could express it in the form of an explicit radical. Ramanujan introduced the functions $P, Q, R$ given by \begin{align} P(q) &= 1 - 24\sum_{i = 1}^{\infty}\frac{iq^{i}}{1 - q^{i}}\tag{5a}\\ Q(q) &= 1 + 240\sum_{i = 1}^{\infty}\frac{i^{3}q^{i}}{1 - q^{i}}\tag{5b}\\ R(q) &= 1 - 504\sum_{i = 1}^{\infty}\frac{i^{5}q^{i}}{1 - q^{i}}\tag{5c} \end{align} and obtained the formulas \begin{align} P(q^{2}) &= \left(\frac{2K}{\pi}\right)^{2}\left(\frac{3E}{K} + k^{2} - 2\right)\tag{6a}\\ Q(q^{2}) &= \left(\frac{2K}{\pi}\right)^{4}\left(1 - k^{2} + k^{4}\right)\tag{6b}\\ R(q^{2}) &= \left(\frac{2K}{\pi}\right)^{6}(1 + k^{2})(1 - 2k^{2})\left(1 - \frac{k^{2}}{2}\right)\tag{6c} \end{align} Ramanujan was not happy with equation $(6a)$ and wanted it to contain only an expression in $k$ apart from the factor $(2K/\pi)^{2}$ and he found that for each positive rational number $n$ the expression $nP(q^{2n}) - P(q^{2})$ could be written as $(2K/\pi)^{2}F_{n}(k)$ where $F_{n}(k)$ is a complicated but algebraic function of $k$. Ramanujan obtained explicit radical expressions for $F_{n} (k) $ for some integer values of $n$. For small values of $n$ like $2,3,5$ it is not difficult to get these expressions, but no one knows how Ramanujan obtained the expressions for larger values of $n$. The method described by Ramanujan to evaluate these expressions is very tiresome and requires great skill in manipulation of radicals. Using this result he proved the following significant theorem:

Ramanujan's Theorem: If $n$ is a positive rational number and $q = e^{-\pi\sqrt{n}}$ then $$P(q^{2}) - \frac{3}{\pi\sqrt{n}} = \left(\frac{2K}{\pi}\right)^{2}A_{n}$$ where $A_{n}$ is an algebraic number dependent on $n$.

If $n$ is large then $q = e^{-\pi\sqrt{n}}$ is small and the functions $P, Q, R$ are close to $1$ and hence can be approximated by $1$. Also note that if $n$ is positive rational then $k$ is algebraic and in the light of above theorem and equations $(6b), (6c)$ we can see that $$\left(P(q^{2}) - \frac{3}{\pi\sqrt{n}}\right)Q(q^{2}) = B_{n}R(q^{2})\tag{7}$$ where $B_{n}$ is an algebraic number. Using approximation $P, Q, R \approx 1$ for large $n$ we get the approximation for $\pi$ $$\pi \approx \frac{3}{(1 - B_{n})\sqrt{n}}\tag{8}$$ The approximation in question is for $n = 25$. It is difficult to calculate the expression $B_{n}$ for generic $n$ but for $n = 25$ we can get the value the value of $B_{n}$ (with significant computational effort) so that the desired approximation for $\pi$ is obtained. Ramanujan also gives the estimate for error as $24\pi(10\pi\sqrt{n} - 31)e^{-2\pi\sqrt{n}}$.

Ramanujan also uses the identity $$C(q)=1-24\sum_{i=1}^{\infty}\frac{(2i-1)q^{2i-1}}{1-q^{2i-1}}=\left(\frac{2K}{\pi}\right)^{2}(1-2k^{2})\tag{9}$$ to get the equation $$P(q^{2})-\frac{3}{\pi\sqrt{n}}=C_{n}C(q)$$ where $C_{n} $ is algebraic and this leads to the approximation $$\pi\approx\frac{3}{(1-C_{n})\sqrt{n}}\tag{10}$$ and for $n=25$ this gives the much simpler but less accurate approximation $$\pi\approx \frac{9}{5}+\sqrt{\frac{9}{5}}\tag{11}$$ which is mentioned in another answer here.

Proof for Ramanujan's theorem can be found in my blog post. These techniques are further used by Ramanujan to give many famous series for $1/\pi$ (also covered in the blog post).

Update: Finally all the ingredients necessary to evaluate $B_n$ above for $n=25$ are available on math.se. I perform the desired evaluation below.

To set up things let us assume that $l, \gamma, k$ are elliptic moduli corresponding to the nomes $e^{-5\pi},e^{-\pi},e^{-\pi/5}$ and let $q=e^{-\pi/5}$. Let the corresponding elliptic integrals be denoted by $L, \Gamma, K$. It is well known that $$\gamma =\frac{1}{\sqrt {2}}, 2l^2=1-\sqrt{1-\phi^{-24}},2k^2=1+\sqrt{1-\phi^{-24}} \tag{12}$$ and we have $k'=l, l'=k, \gamma'=\gamma$ ($\phi$ is the golden ratio).

Also from this answer we have $$\frac{L} {\Gamma} =\frac{\sqrt{5}+2}{5}\tag{13}$$ Clearly from the transformation formula of theta functions it is obvious that $K/L=5$.

We now use the formula for $F_n(k) $ for $n=5$ given by Ramanujan and obtain $$5P(q^{10})-P(q^2)=\frac{4K\Gamma}{\pi^2}(3+k\gamma+k'\gamma')\sqrt{\frac{1+k\gamma+k'\gamma'}{2}}=25(\sqrt{5}-2)\left(\frac{2L}{\pi}\right)^2(3+(k+k')\gamma)\sqrt{\frac{1+(k+k')\gamma}{2}}\tag{14a}$$ and $$5P(q^{50})-P(q^{10})=\frac{4\Gamma L} {\pi^2}(3+l\gamma+l'\gamma')\sqrt{\frac{1+l\gamma+l'\gamma'}{2}}=5(\sqrt{5}-2)\left(\frac{2L}{\pi}\right)^2(3+(l+l')\gamma)\sqrt{\frac{1+(l+l')\gamma}{2}}\tag{14b}$$ Multiplying $(14b)$ by $5$ and adding the result to $(14a)$ we get $$25P(q^{50})-P(q^2)=50(\sqrt{5}-2)\left(\frac{2L}{\pi}\right)^2(3+(l+l')\gamma)\sqrt{\frac{1+(l+l')\gamma}{2}}$$ We have $$l+l'=\sqrt{l^2+l'^2+2ll'}=\sqrt{1+\phi^{-12}}$$ and thus $$(l+l') \gamma=\sqrt{\frac{1+\phi^{-12}}{2}}=3\phi^{-3}=3(2\phi-3)$$ and we finally arrive at $$25P(q^{50})-P(q^2)=50(\sqrt{5}-2)\left(\frac{2L}{\pi}\right) ^2(6\phi-6)\sqrt{3\phi-4}$$ On the other hand using logarithmic differentiation of transformation formula for Dedekind eta function allows us to prove that $$25P(q^{50})+P(q^2)=\frac{30}{\pi}$$ And adding last two equations we get $$P(q^{50})-\frac{3}{5\pi}=\left(\frac{2L}{\pi}\right) ^26(5-3\phi)\sqrt{3\phi-4}\tag{15}$$ And we have the standard relations $$Q(q^{50})=\left(\frac{2L}{\pi}\right)^4(1-l^2l'^2)=\left(\frac{2L}{\pi}\right)^4\left(\frac{4-\phi^{-24}}{4}\right)$$ and $$R(q^{50})=\left(\frac{2L}{\pi}\right)^6(1-2l^2)(1+(l^2l'^2/2))=\left(\frac{2L}{\pi}\right)^6\sqrt {1-\phi^{-24}}\left(\frac{8+\phi^{-24}}{8}\right)$$

The expression $B_n$ equals $$\frac{(P(q^{50})-(3/5\pi))Q(q^{50})} {R(q^{50})} $$ and using previous equations we get $$B_n=\frac{12(5-3\phi) \sqrt{3\phi-4}}{\sqrt{1-\phi^{-24}}}\cdot\frac{4-\phi^{-24}}{8+\phi^{-24}}$$ It turns out that the first fraction evaluates to $\phi$ (after some effort) and hence $$B_n=\frac{4\phi-\phi^{-23}}{8+\phi^{-24}}$$ and $$1-B_n=\frac{8-4\phi+\phi^{-23}+\phi^{-24}}{8+\phi^{-24}}=\frac{8-4\phi+\phi^{-22}}{8+\phi^{-24}}$$ The approximation for $\pi$ is now given by $$\pi\approx\frac{3}{5(1-B_n)}=\frac{3}{5}\cdot \frac{8+\phi^{-24}} {8-4\phi+\phi^{-22}}$$ Let $a_m$ denote Fibonacci sequence with $a_0=0,a_1=1$ and then $a_{-m} =(-1)^{m+1}a_m$ and we have $$\phi^m=a_m\phi+a_{m-1}$$ Using this relation we get $$\pi\approx\frac{3}{5}\cdot \frac{8+a_{25} -a_{24} \phi}{8+a_{23} -(4+a_{22})\phi}$$ Using a table of Fibonacci sequence we get $$a_{22} =17711,a_{23} =28657,a_{24}=46368,a_{25}=75025$$ and then $$\pi\approx\frac{3}{5}\cdot\frac{75033-46368\phi}{28665-17715\phi}=\frac{63}{25}\cdot\frac{3573-2208\phi}{5733-3543\phi}=\frac{63}{25}\cdot\frac{1646-736\sqrt{5}}{2641-1181\sqrt{5}} $$ and by this answer the second fraction equals $$\frac{17+15\sqrt{5}}{7+15\sqrt{5}}$$ Thus we obtain the famous approximation by Ramanujan $$\pi\approx \frac{63}{25}\cdot\frac{17+15\sqrt{5}}{7+15\sqrt{5}}$$ The approximation in equation $(11)$ can also be obtained in similar manner (with much less effort) by evaluating $C_n$ for $n=25$.

Ramanujan obtained most of his results (including the above) using algebraic manipulation combined with processes of calculus (differentiation and integration) and rarely used any sophisticated tools. Modern approach to Ramanujan's mathematics is so unlike his methods and is based on sophisticated techniques and theories (modular forms) which fail to generate any interesting results like those given by Ramanujan and are almost always dependent on software like Maple, Macsyma or Mathematica for any serious calculation.

When I first came across your question, I thought it was a modern-day approximation by somebody using a computer. But when d125q pointed out it was by Ramanujan, then I figured out he must have used a systematic method.

One way is to use a Ramanujan-Sato pi formula like,

$$\frac{1}{\pi} = \frac{1}{16}\sum_{n=0}^\infty \frac{(2n)!^3}{n!^6}\frac{(42\phi-6)n+(5\phi-3)}{(2^{12}\phi^8)^n}\tag1$$

where $\phi=\frac{1+\sqrt{5}}{2}$, and truncate it as for finite number of terms. For example, using just $n=0\;\text{to}\;1$, and getting the reciprocal, it yields,

$$\pi \approx \frac{2^{13}}{3(-383+560\sqrt{5})}$$

It is only good for $10^{-7}$, and the next is $10^{-10}$, but there is an infinite choice of $n$.

There are three formulas in Mathworld that use a $\sqrt{5}$, including a version of $(1)$. And there is also a fourth. However, Ramanujan must have known still another formula because I can't get the approximation in your post by truncating any of the four.

This is not yet a complete answer, but may be useful.

The largest root of the polynomial

$$269x^2-503x+209$$ is $$\frac{17+15\sqrt{5}}{7+15\sqrt{5}}$$

Changing the polynomial to $$(25)^2\times269x^2-25\times63\times 503 x+63^2\times 209$$

modifies the root to the approximation given. $$\pi\approx\frac{63}{25}\times\frac{17+15\sqrt{5}}{7+15\sqrt{5}}$$

In terms of the golden ratio, this is $$\pi\approx\frac{63}{25}\left(1+\dfrac{1}{3\phi-\dfrac{4}{5}}\right)$$

This seems related to the simpler polynomial $$25x^2-90x+36$$ that can have its coefficients factored as $$5^2x^2-5\times 6 \times 3x+6^2$$

and has the largest root $$\pi\approx \frac{9}{5}+\sqrt{\frac{9}{5}},$$

which is another approximation by Ramanujan that has a similar expression in terms of the golden ratio.

$$\pi\approx \dfrac{6}{5}\left(1+\dfrac{1}{\phi-1}\right)$$

This suggests an intermediate approximation of the form $$\pi \approx r_2\left(1+ \dfrac{1}{2\phi+d_2}\right)$$

and hopefully better precision with an expression similar to $$\pi \approx r_4\left(1+ \dfrac{1}{4\phi+d_4}\right)$$

or $$\pi \approx r_5\left(1+ \dfrac{1}{5\phi-d_5}\right)$$ if even multiples of $\phi$ are of no use. Fractions $r_n$ and $d_n$ are yet to be found.