Inverse of difference of two digamma functions

Question

I recently encountered the expression below for which I was interested in solving for $x$:

\begin{equation} \psi(x+n+1) - \psi(x+1) =y \end{equation}

$\psi$ is the digamma function, $n$ is a positive integer and $x,y>0$. Initially, I solved this numerically but I was interested to see if there was the possibility to solve this, either in closed form or as a series. I tried to find some way by using the relationship to harmonic numbers, e.g. rewriting as

\begin{equation} \psi(x+n+1) - \psi(x+1) = \sum_{k=1}^{n} \frac{1}{k+x} = \int_0^1 \frac{t^x(1-t^n)}{1-t}dt \end{equation}

to see if there was something which I could do, but I didn't notice anything apparent. The only thing I've come up with is that for large values of $x$, $y$ is small and using $\psi(x) \sim\ln(x)$, that $x \approx \dfrac{n}{e^y - 1}$.

I'm happy for any suggestions on how to proceed!

Answer 1

$\def\F{\operatorname F}\def\H{\operatorname H}$Derivation: After different Lagrange reversion setups, this one had the simplest derivation without needing to factor. Additionally, less series are possible with Tschirnhausen transformations, but applying them to a degree $m$ polynomial in general is complicated. The harmonic number H$_m$ set up includes $(0,0)$: $$e^{\H_m-\sum\limits_{k=1}^m\frac1{y+k}}-1=x\implies y=x+\sum_{n=1}^\infty \frac1{n!}\frac{d^{n-1}}{dx^{n-1}}\left(x+1-e^{H_m}\prod_{k=1}^me^{-\frac1{x+k}}\right)^n$$ shown here, so here is a solution directly without reducing the equation. Using the binomial theorem on the derivative: $$\frac{d^{n-1}}{dx^{n-1}}\left(x+1-e^{\H_m}\prod_{k=1}^me^{-\frac1{x+k}}\right)^n=\sum_{j=0}^n\binom nj\left(-e^{\H_m}\right)^j\frac{d^{n-1}}{dx^{n-1}}\left((x+1)^{n-j} \prod_{k=1}^me^{-\frac j{x+k}}\right)$$ Applying the general Leibniz rule and factorial power $a^{(b)}$: $$\frac{d^{n-1}}{dx^{n-1}}\left((x+1)^{n-j} \prod_{k=1}^me^{-\frac j{x+k}}\right)=\sum_{v=0}^{n-1} \binom{n-1}v\frac{d^{n-1-v}}{dx^{n-1-v}}((x+1)^{n-j}) \frac{d^v}{dx^v}\left(\prod_{k=1}^me^{-\frac j{x+k}}\right)= \sum_{n_0=0}^{n-1} \binom{n-1}{n_0} (x+1)^{n_0-j+1}(n-j)^{(n-1-n_0)}\frac{d^{n_0}}{dx^{n_0}}\left(\prod_{k=1}^me^{-\frac j{x+k}}\right)$$ Also, we use the generalized product rule: $$\frac{d^{n_0}}{dx^{n_0}}\left(\prod_{k=1}^me^{-\frac j{x+k}}\right)= \sum_{n_1=0}^{n_0}\sum_{n_2=0}^{n_1}\dots\sum_{n_{m-1}=0}^{n_{m-2}}\frac{{n_0}!f_1^{(n_0-n_1)}f_m^{(n_{m-1})}}{(n_0-n_1)!n_{m-1}!}\prod_{k=2}^{m-1}\frac{f_k^{(n_{k-1}-n_k)}}{(n_{k-1}-n_k)!} $$ A formula with the confluent hypergeometric $_1\F_1(a;b;z)$ function and Kronecker $\delta_x,\delta_{x,y}$: $$\frac{d^r}{dx^r}e^{-\frac j{x+k}}=\sum_{u=0}^\infty\frac{(-j)^u}{u!}\frac{d^r}{dx^r}(x+k)^{-u}=\sum_{u=0}^\infty \frac{(-j)^u (-u)^{(r)}(x+k)^{-r-u}}{u!}=(x+k)^{-r}0^{(r)}-j(x+k)^{-r-1}(-1)^{(r)}\,_1\text F_1\left(r+1;2;-\frac j{x+k}\right)=\delta_r-j(-1)^r r!(x+k)^{-r-1}\,_1\text F_1\left(r+1;2;-\frac j{x+k}\right) $$ Simple Verification: Unfortunately, Wolfram Alpha cannot evaluate the full sums, but does for a partly written sum. Let’s verify that the above steps work for $m=2$. Testing $n=3,x=0.3$ gives matching series coefficients. Thereafter, testing $\frac{d^v}{dx^v}\left(e^{-\frac j{x+1}-\frac j{x+2}}\right),j=4,v=3,x=0.3$ gives the same result. Similarly when $m=3$, we get matching series coefficients. Testing $\frac{d^v}{dx^v}\left(e^{-\frac j{x+1}-\frac j{x+2}-\frac j{x+3}}\right),v=4,j=5,x=0.3$ gives the same result. Wolfram Alpha cannot evaluate the derivatives as confluent hypergeometric functions, but the derived $\frac{d^r}{dx^r} e^{-\frac j{x+k}}$ formula is true presents no convergence issues. Other values of $v,j,x$, including $v,j=0$, work when $m=2,3$ and combining verification links give evidence for the Lagrange reversion result working when $m>3$ Formula: Therefore we should get: $$\sum_{k=1}^m\frac1{y+k}=x\implies y=\H_m-\ln\left(1+\sum_{n=1}^\infty\sum_{j=0}^n\binom nj\left(-e^{\H_m}\right)^j \sum_{n_0=0}^{n-1}(x+1)^{n_0-j+1}(n-1)^{(n_0)}(n-j)^{(n-1-n_0)} \sum_{n_1=0}^{n_0}\sum_{n_2=0}^{n_1}\dots\sum_{n_{m-1}=0}^{n_{m-2}}\left(\frac{\delta_{n_0,n_1}} {(n_0-n_1)!}-j(-1)^{n_0-n_1}(x+1)^{n_1-n_0-1}\,_1\text F_1\left(n_0-n_1+1;2;-\frac j{x+1}\right)\right) \left(\frac{\delta_{n_{m-1}}} {n_{m-1}!}-j(-1)^{n_{m-1}}(x+m)^{-n_{m-1}-1}\,_1\text F_1\left(n_{m-1}+1;2;-\frac j{x+m}\right)\right)\prod_{k=2}^{m-1}\left(\frac{\delta_{n_{k-1},n_k}} {(n_{k-1}-n_k)!}-j(-1)^{n_{k-1}-n_k}(x+k)^{n_k-n_{k-1}-1}\,_1\text F_1\left(n_{k-1}-n_k +1;2;-\frac j{x+k}\right)\right) \right)= \boxed{\H_m-\ln\left(1+\sum_{n=1}^\infty\sum_{j=0}^n\sum_{n_0=0}^{n-1}\frac{\Gamma(n)n! (x+1)^{n_0-j+1} \left(-e^{\H_m}\right)^j}{(v-j+1)!\Gamma(n-v)j!}\sum_{0\le n_{m-1}\le\dots\le n_0}\left(\frac{\delta_{n_{m-1}}} {n_{m-1}!}-j(-1)^{n_{m-1}}(x+m)^{-n_{m-1}-1}\,_1\text F_1\left(n_{m-1}+1;2;-\frac j{x+m}\right)\right)\prod_{k=1}^{m-1}\left(\frac{\delta_{n_{k-1},n_k}} {(n_{k-1}-n_k)!}-j(-1)^{n_{k-1}-n_k}(x+k)^{n_k-n_{k-1}-1}\,_1\text F_1\left(n_{k-1}-n_k +1;2;-\frac j{x+k}\right)\right) \right)}$$ Closed form?: Since solving $\sum\limits_{k=1}^m\frac1{y+k}=x$ is solving a polynomial, you use the Riemann theta function, but this is too advanced for some.

Answer 2

$$\psi(x+n+1) - \psi(x+1) =y$$

For large values of $x$, the lhs write $$\frac{n}{x}-\frac{n (n+1)}{2 x^2}+\frac{n (n+1) (2 n+1)}{6 x^3}-\frac{n^2 (n+1)^2}{4 x^4}+O\left(\frac{1}{x^5}\right)$$ Using series reversion $$x_{(1)}=\frac{n}{y}-\frac{n+1}{2}+\frac{\left(n^2-1\right)}{12 n}y+O\left(y^3\right)$$

Comparing with your approximation $$x_{(2)}= \dfrac{n}{e^y - 1}=\frac{n}{y}-\frac{n}{2}+\frac{n }{12}y+O\left(y^3\right)$$ It is really good since $$x_{(2)}-x_{(1)}=-\frac{1}{2}-\frac{y}{12 n}+O\left(y^3\right)$$

Using the above for $n=1000$ and $y=1$, the exact solution is $581.477$ while $x_{(2)}=582.833$ and $\frac{1000}{e-1}=581.977$.

Edit

We could do better using more terms in the series expansion, then series reversion and transformation of the result into a $[n,n+1]$ Padé approximant.

This would give $$x_{(3)}=\frac{60 (n-1) n^3 -24 n^2 \left(n^2-4\right)y+3 n \left(n^3-3 n^2-9 n-13\right)y^2} {60 (n-1) n^2 y+6 n \left(n^2+11\right)y^2+(n-1) \left(n^2+11\right)y^3 }$$

For the worked example, this gives $x_{(3)}=581.590$.