Numerical Methods – Approximating the Digamma Function

Question

I am trying to approximate the digamma function in order to graph it in latex. I've already found an approximation of the gamma function from the Tex SE:
$$
\Gamma(z)=
(2.506628274631\sqrt{(1/z)} + 0.20888568(1/z)^{(1.5)} + 0.00870357(1/z)^{(2.5)} – \frac{(174.2106599*(1/z)^{(3.5)})}{25920} – \frac{(715.6423511(1/z)^{(4.5)})}{1244160)}\exp((-ln(1/z)-1)z)
$$
and for the digamma, $\psi$, ive been using an approximation of the derivative:
$$
\psi(x)=\frac{\ln \Gamma(x+0.0001)-\ln \Gamma(x-0.0001)}{0.0002}
$$
But im wondering if there is a better way

Answer 1

Yes, there is a formula $$ψ(x) \approx \log(x) - \frac{1}{2x} - \frac{1}{12x^2} +\frac{1}{120x^4} -\frac{1}{252x^6} +\frac{1}{240x^8} -\frac{5}{660x^{10}}+\frac{691}{32760x^{12}} -\frac{1}{12x^{14}}$$ This is especially accurate for larger values of $x$. If $x$ is small, you can shift $x$ to a higher value using the relation $$ψ(x+1) = \frac{1}{x} + ψ(x)$$

You can find my source here, under the "Computation and Approximation" section.

Edit: Accuracy of Approximation
This series is quite accurate for $x$ in the interval $(1,\infty)$. Given the OP's bounds of $x = (0,20)$, we can find the difference between the expected and calculated results (note that $S(x)$ is the series in question) $$ψ(1) - S(1) = 0.0674022$$ $$ψ(20) - S(20) = 4.44089×10^{-16}$$ From these results we find that we get exactly what we expected: the series is decent in this interval, but is MUCH better for larger $x$. However, we get a problem for $x = \frac{1}{2}$: $$ψ(\frac{1}{2}) - S(\frac{1}{2}) = 1285.81$$ Obviously we are way off here. However, some context is useful. Our series blows up for $x<1$ towards $-\infty$, and so we can't expect anything useful here. Nevertheless, we can use the shift formula above to mitigate this issue, rewriting it as $$ψ(x) = ψ(x+1) - \frac{1}{x}$$ Applying this just once for $x=\frac{1}{2}$, we get a difference of $1.9816604...$, an approximation $\approx 650$ times better. Thus, we see that we can use the series provided to get a decent approximation of the graph of $ψ(x)$ for all $x>0$.