Approximate the first solution of $\sum_{k=0}^n \frac 1{x-k}=0$

Question

I would like to obtain formally a "good" approximation of the first root of function $$F_n(x)=\sum_{k=0}^n \frac 1{x-k} $$ where $n$ can be extremely large.

For this problem, as earlier when solving it numerically, I introduce two auxiliary functions $$G_n(x)=x(x-1)f_n(x) \qquad \text{and} \qquad J_n(x)=xf_n(x)$$

Using the first iteration of Newton method with $x_0=0$, we easily generated the bounds $$x_g^{(n)}= \frac{1}{H_n+1} < x^{(n)} < \frac{1}{H_n}=x_j^{(n)}$$

A few numbers for illustration

$$\left( \begin{array}{cccc} n & x_g^{(n)} & \text{solution} & x_j ^{(n)}\\ 100 & 0.161619 & 0.180829 & 0.192776 \\ 200 & 0.145390 & 0.161805 & 0.170125 \\ 300 & 0.137312 & 0.152321 & 0.159168 \\ 400 & 0.132102 & 0.146205 & 0.152209 \\ 500 & 0.128323 & 0.141774 & 0.147214 \\ 600 & 0.125392 & 0.138339 & 0.143370 \\ 700 & 0.123016 & 0.135557 & 0.140272 \\ 800 & 0.121029 & 0.133232 & 0.137695 \\ 900 & 0.119329 & 0.131244 & 0.135498 \\ 1000 & 0.117848 & 0.129513 & 0.133592 \\ \end{array} \right)$$

What we have is $\forall n$ $$G_n\Big(x_g^{(n)}\Big) <0 \qquad \text{and} \qquad G_n\Big(x_g^{(n)}\Big)\,\,G''_n\Big(x_g^{(n)}\Big) >0$$

$$J_n\Big(x_h^{(n)}\Big) <0 \qquad \text{and} \qquad J_n\Big(x_j^{(n)}\Big)\,\,J''_n\Big(x_j^{(n)}\Big) >0$$ which, by Darboux theorem, confirm that $x_g^{(n)}$ is an underestimate and $x_h^{(n)}$ is an overestimate of the searched solution.

The problem is that the second iterate of Newton method is extremely expensive in terms of computer resources.

My question is

Would it be possible to define a function $f(n)$ such that the solution would be $$x^{(n)} \sim \frac{1}{H_n+f(n)}\qquad \text{with} \qquad 0 < f(n) < 1$$ without resorting to curve fitting procedure ? I am not looking for any extravagant accuracy (one percent maximum relative error would be sufficient).

Edit

An empirical work seems to show that a decent estimate could be $$x^{n}_{1}=\frac{1}{H_n+\frac{16}{9\,H_n}}$$

Even if it does not means too much, over the $n$ used in the table

$$\sum G_n\Big(x_g^{(n)}\Big) =-0.863239$$ $$\sum G_n\Big(x_j^{(n)}\Big) =+0.363728$$ $$\sum G_n\Big(x^{(n)}_{1}\Big)=-0.0012461$$

Repeating the calculations $$\left( \begin{array}{ccc} n & x_1^{(n)} & \text{solution}\\ 100 & 0.180829 & 0.180829 \\ 200 & 0.161800 & 0.161805 \\ 300 & 0.152308 & 0.152321 \\ 400 & 0.146188 & 0.146205 \\ 500 & 0.141753 & 0.141774 \\ 600 & 0.138315 & 0.138339 \\ 700 & 0.135531 & 0.135557 \\ 800 & 0.133205 & 0.133232 \\ 900 & 0.131215 & 0.131244 \\ 1000 & 0.129484 & 0.129513 \\ \end{array} \right)$$

Could a continued fraction be a way to go ? If it is, how to build it?

Answer 1

Not an answer, but:

We have $$x = \frac{1}{g(n,x)} \tag 1$$

for $x \in (0,1)$ where

$$ \begin{align} g(n,x) &= \frac{1}{1-x} + \frac{1}{2-x} + \cdots \frac{1}{n-z}\\ &= H_n + x \, \zeta(2,n) + x^2 \, \zeta(3,n) + x^3 \, \zeta(4,n)\cdots \tag 2 \end{align}$$

where $\zeta(k,n) = \sum_{j=1}^n j^{-k}$, so eg $\zeta(2,n)\approx \frac{\pi^2}{6}-\frac{1}{n}$.

This does not look very useful. But at least it gives some support for the pseudo continous fraction approximation:

$$ x \approx \frac{1}{H_n + \frac{\pi^2/6}{H_n} \cdots} \tag 3$$

(Notice that $\pi^2/6 = 1.645\cdots$ is not very far from $16/9=1.777$)

We can also write

$$g(n,x) = \psi(n+1-x)-\psi(1-x) \tag 4$$

where $\psi()$ is the digamma function. (Sanity check).

We have, for large $n$ and small $x$ (more precisely, for $x \sim 1/\log n$):

$$ \psi(n+1-x) = \log(n) +\frac{1-2x}{2n} + o(1/n)$$

$$ \psi(1-x) = -\gamma - \frac{\pi^2}{6} x +O(x^2) $$

Hence, to first order, we get something similar:

$$x \approx \frac{1}{\log n + \gamma + \frac{\pi^2}{6} x} \tag 5$$

---

Update (still not a real answer): a numerical recipe. Using the following approximations for the digamma function, for $z\to \infty$ and for $z \to 1$:

$$g_1(z) = \log(z) - \frac{1}{2z} - \frac{1}{12 z^2} + \frac{1}{120 z^4} \tag 5$$

$$g_0(z) = - \gamma + \frac{\pi^2}{6}(z-1)+ \alpha (z-1)^2 + \frac{\pi^4}{90} (z-1)^3+ \beta (z-1)^4 \tag 6$$

with $\gamma = 0.57721566490153286$ (Euler's constant), $\alpha = \frac{1}{2}\psi^{(2)}(1) = -1.20205690315959$ and $\beta = \frac{1}{24}\psi^{(4)}(1) = -1.03692775514337$, the recursion

$$ x \leftarrow \frac{1}{g_1(n+1-x) - g_0(1-x)} \tag 7$$

starting with $x = \frac{1}{\gamma + \log n}$ converges in few ($\sim 10$) iterations to good approximations, especially for large $n$:

    n      approx     exact 
    10    0.285602   0.285421
    30    0.225494   0.225460
   100    0.180836   0.180829
   300    0.152323   0.152321
  1000    0.129514   0.129514
 10000    0.100341635  0.100341522
100000    0.081744233  0.081744206

Of course, if one has a reliable and efficient implementation of the digamma function for all the real range, one can just use in place of $g_0$ and $g_1$ in $(7)$.

Answer 2

After @Sharky Kesa's answer

$$f(x)=\psi(1+x) + \pi \cot(\pi x) + \log(1 - \frac{x}{n})$$ Expanding as series around $x=0$

$$f(x)=\frac 1x-\gamma-\sum_{k=1}^\infty \Big(\frac 1{k\,n^k}+\zeta(k+1)\Big)\,x^k$$ we can use power series reversion to obtain $$x=\sum_{k=1}^\infty \frac {a_k}{\log^k(n)}$$ and we know all the coefficients $a_k$ using the explicit formula for the $k^{\text{th}}$ term using Morse and Feshbach explicit formula.

Limiting the initial expansion to $O(x^{11})$ , the results $$\left( \begin{array}{ccc} n & \text{estimate} & \text{solution}\\ 100 & 0.180868 & 0.180829 \\ 200 & 0.161826 & 0.161805 \\ 300 & 0.152335 & 0.152321 \\ 400 & 0.146216 & 0.146205 \\ 500 & 0.141782 & 0.141774 \\ 600 & 0.138346 & 0.138339 \\ 700 & 0.135562 & 0.135557 \\ 800 & 0.133237 & 0.133232 \\ 900 & 0.131248 & 0.131244 \\ 1000 & 0.129517 & 0.129513 \\ \end{array} \right)$$

Answer 3

Not a full answer but just some of the work I have on this question, not sure how to proceed further:

Observe that \begin{align*} F_n(x) &= -\sum_{k = 0}^n \frac{1}{k - x}\\ &= -\sum_{k = 0}^n \frac{\frac{1}{k}}{1 - \frac{x}{k}}\\ &= -\sum_{k = 0}^n \sum_{j = 0}^\infty \frac{x^j}{k^{j+1}}\\ &= -\sum_{j = 0}^\infty x^j H_n^{(j+1)}\\ \end{align*} Setting this equal to $0$, we get \begin{align*} 0 &= \sum_{j = 0}^\infty x^j H_n^{(j+1)}\\ &= \sum_{j = 0}^\infty x^j \left (\sum_{k = 0}^\infty \frac{1}{k^{j+1}} - \sum_{k = n+1}^\infty \frac{1}{k^{j+1}} \right )\\ &\approx \sum_{j = 0}^\infty x^j \left (\zeta(j + 1) - \int_n^\infty \frac{1}{t^{j+1}} \ \mathrm{d}t \right )\\ &= H_n + \sum_{j = 1}^\infty x^j \left (\zeta(j + 1) - \frac{1}{jn^j}\right )\\ &= H_n - \gamma - \psi(- x) + \log\left (1 - \frac{x}{n} \right )\\ \log(n - x) &\approx \psi(-x)\\ &= \psi(1 + x) + \pi \cot(\pi x)\\ &\approx \pi \cot(\pi x) - \log\left (e^{\frac{1}{x+1}}-1\right )\\ \log(n) &\approx \pi \cot(\pi x) + \log \left (\frac{1 - \frac{x}{n}}{e^{\frac{1}{x + 1}} - 1} \right ) \end{align*} For large $n$, the $\log$ term is equal $-\log(e - 1) + O\left (\frac{x}{n}\right )$. Hence, an approximation for $x$ could be \begin{align*} x &\approx \frac{1}{\pi} \cot^{-1}\left (\frac{\log((e - 1) n)}{\pi} \right ) \end{align*} This doesn't seem to be very accurate for $n = 100000$ though, it gets $x = 0.081153030\dots$, but the exact solution is $x = 0.081744206\dots$