Is there a closed form for the polygamma function? Mathematica gave me that
$$
\sum_{k=n}^\infty \frac1{k^2} = \texttt{PolyGamma[1,n]}.
$$
However, in all my attempts to simplify and approximate the number as a decimal, it kept leaving it in terms of the "PolyGamma" function. I understand that
\begin{align}
\texttt{PolyGamma[1,1]} &= \psi^{(1)}(1) \\
&= (-1)^2 \int_0^\infty \frac{te^{-t}}{1-e^{-t}}\ \mathsf dt\\
&= \frac{\pi^2}6,
\end{align}
and in general
\begin{align}
\texttt{PolyGamma[1,n]} &= \psi^{(1)}(n)\\ &=  \int_0^\infty \frac{te^{-nt}}{1-e^{-t}}\ \mathsf dt\\ &= -\int_0^1 \frac{t^{n-1}}{1-t}\log t\ \mathsf dt.
\end{align}
Is this just not an integral that cannot be solved analytically? Wolfram Alpha gives me the series expansion
$$
\frac1{n^2} + \frac{\pi^2}6 + n \psi^{(2)}(1) + \frac{\pi^4n^2}{30} + O(n^3),
$$
which is fairly useless considering $n$ as a continuous variable, as this is the series expansion about $0$.
If this integral cannot be computed analytically, then what is a good approximation for $\texttt{PolyGamma[1,n]}$ as a function of $n$? I see that it can be computed exactly for any fixed $n$. But I am interested in the map $n\mapsto \psi^{(1)}(n)$.
Edit: From Wikipedia it appears that an easy way to approximate the trigamma function is to take the derivative of the series expansion of the digamma function:
$$
\psi^{(1)}(n) = \frac1n + \frac1{2x^3} + \frac1{6x^3} - \frac1{30x^5} + \frac1{42x^7} - \frac1{30x^9} + \frac5{66x^{11}} - \frac{691}{2730x^{13}} + \frac7{6x^{15}} + O(n^{17})
$$
However I see no clear pattern as to the coefficients in this Laurent series.

The Euler-Maclaurin Sum Formula gives the asymptotic approximation
$$
\begin{align}
\sum_{k=n}^\infty\frac1{k^2}
&\sim\sum_{k=0}^\infty\frac{B_{2k}}{n^{2k+1}}\\
&=\frac1n+\frac1{2n^2}+\frac1{6n^3}-\frac1{30n^5}+\frac1{42n^7}-\frac1{30n^9}+O\!\left(\frac1{n^{11}}\right)
\end{align}
$$
where $B_n$ are the Bernoulli numbers.
Note that this is an asymptotic approximation; that is, the series does not converge.

If you need an exact answer, then the best you can do would be
$$
\begin{align}
\sum_{k=n}^\infty\frac1{k^2}
&=\frac{\pi^2}6-\sum_{k=1}^{n-1}\frac1{k^2}\\
&=\frac{\pi^2}6-H^{(2)}_{n-1}
\end{align}
$$
where $H^{(2)}_n$ are the Generalized Harmonic Numbers.
