Approximating the Digamma fucntion near 1

Question

Peace be upon you,

I had the following system of equations to be solved \begin{align*} \begin{cases} \psi(\alpha)-\psi(\alpha+\beta)=c_1\\ \psi(\beta)-\psi(\alpha+\beta)=c_2 \end{cases} \end{align*} where $c_1$ and $c_2$ are constants (I should solve the above, for miscellaneous constant couples). First of all, I wanted to solve it analytically so I tested some ideas and then asked for the solution in the question "A late-diverging "approximating solution" for a system of functional equations"; but after visiting no answer I proceeded to solve it numerically, using the popular approximation of the Digamma function \begin{align*} \psi(x)= \ln(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}}+O\left(\frac{1}{x^{16}}\right) \end{align*} as you can see from the approximation order, the error in 1 can even be 1! and for the values less that 1 the story is worse; while the answer of my system ($\alpha$ and $\beta$) maybe in this interval ($(0,1]$). I wonder if someone can present a better approximation of the Beta function which can cope my task for such values.

Any lighting up ideas?

Answer 1

An approximation close to $1$ (i.e. $x=1+\epsilon$ ) is given below (from WolframAlpha) :

For more accuracy, you can use the Taylor's expansion : $$\psi(1+\epsilon)=\Sigma_{n=0}^{\infty} \frac{\psi^{(n)}(1)}{n!}\epsilon^n$$ with numerical values of $\psi^{(n)}(1)$

Answer 2

I dare to post a second answer because it is important to clearly separate the answers to the first and the second question raised. The subject of one question was the analytical properties of the digamma function. The subject of the other question is how to compute approximate solutions of a system of equations.

Since the inverse function of the digamma is not related to a referenced special function, there is no hope to find the solution on a closed form. So, I think that numerical calculus is required : There are several methods for numerical solving of non-linear equations, based on more or less sophisticated iterative process.

These processes require a first estimate to start : Sometimes, this is a major difficulty on the practical viewpoint.

In the present case, I propose a method to compute a rough first estimate ( principle shown in attachment). Since the unknowns $\alpha$ and $\beta$ are in the range $0$ to $1$, the Taylor's expansions are made around $\frac{1}{2}$. Of course, if it was known that the expected values are closer to $1$ (for example) the formulas should be rewritten according to it.