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$.

Answer 2

Since you need an implementation for graphing, a simple approximation with relative low accuracy and optimized for absolute error seems suitable. A short search with computer program shows that $$\psi(x) \approx ln (x + a) - \frac{1}{b x}, \space a = 0.4849142940227510, \space b=1.0271785180163817$$ approximates with an absolute error of less than $0.00123$ on $[\frac{1}{2}, 2^{1024})$ when evaluated in IEEE-754 double precision. For arguments $x \lt \frac{1}{2}$, the function can be computed using the reflection formula $\psi(x) = \psi(1-x) + \pi \space cot \space(\pi (1-x))$.