I am looking for the first positive solution $x_n$ of the equation $$f_n(x)=\sum_{k=1}^n\tan(kx)=0 \qquad \qquad (n\geq 2)$$ It is simple to show that $$\frac \pi{2n}<x_n <\frac \pi{2(n-1)}$$ The numerical solution is very easy (and inexpensive in terms of computer resources) to obtain using Newton, Halley or Householder methods but the analytical solution does not seem to be possible as soon as $n\geq 7$ (up to $n=6$, using $x=\tan^{-1}(t)$, the equation factors as products of polynomials of low degree). For example solving $f_6$ reduces to the problem of finding the first positive root of $$t \left(t^2-3\right) \left(t^6-21 t^4+35 t^2-7\right) \left(23 t^8-202 t^6+136 t^4-22 t^2+1\right)=0$$ (which is feasible even if tedious) while for $f_7$, the root which is looked for is solution of $$44 t^{16}-763 t^{14}+4375 t^{12}-10067 t^{10}+9199 t^8-3593 t^6+589 t^4-41 t^2+1=0$$ (the other factors were not printed).
Looking at the numerical results for $2\leq n\leq 10000$, my surprise has been to notice that the solution is $$x_n\approx \frac \pi{2n-1}\tag 1$$ which the reciprocal of the harmonic mean of the bounds ( this is exactly true for $n=2$ ).
For $x_4$, the approximation gives $\approx 0.448799$ while the exact solution is $\cos ^{-1}\left(\frac{\sqrt{9+\sqrt{17}}}{4}\right)\approx 0.437896$.
For $x_{10000}$, the approximation gives $0.000157087$ while the "exact" solution is $0.000157097$.
Using $y_n=\frac \pi{2n-1}$, polynomial regressions $$x_n=\sum_{i=1}^m a_i y^i$$ lead to extremely good fits with highly significant parameters. As an example, for $m=4$ $$\begin{array}{clclclclc} \text{} & \text{Estimate} & \text{Standard Error} & \text{Confidence Interval} \\ a & +0.997756 & 0.0000157 & \{+0.997726,+0.997786\} \\ b & -0.125612 & 0.0001336 & \{-0.125874,-0.125350\} \\ c & +0.202807 & 0.0002973 & \{+0.202224,+0.203390\} \\ d & -0.077174 & 0.0001766 & \{-0.077520,-0.076828\} \\ \end{array}$$ Even with a single term in the regression, $R^2=0.999886$.
My attempts based on Taylor series around the approximation were totally unsuccessful.
My questions are :
1.Is there any way to justify $(1)$ and/or to obtain a better approximation ?
2.Could an asymptotics of the solution be derived ?
Edit
After @Did's very interesting comment, I generated the solutions $x_n$ for $1000\leq n\leq 250000$ $(\Delta n =1000)$ and made a curve fit to the model $$x_n=\frac{\pi }{a n+b+\frac{c}{\log (n)}}$$ The fit is extremely good (all residuals being smaller than $2\times 10^{-12}$) and the parameters are highly significant $$\begin{array}{clclclclc} \text{} & \text{Estimate} & \text{Standard Error} & \text{Confidence Interval} \\ a & +2.00000 & 5.8\times 10^{-9} & \{+2.00000,+2.00000\} \\ b & -0.02216 & 0.000108935 & \{-0.02237,-0.02194\} \\ c & -1.74434 & 0.000719344 & \{-1.74576,-1.74293\} \\ \end{array}$$
This is almost exactly what @Did suggested !!
Update
After Daniel Fischer's answer and Did's comment, I restarted the fitting work for the range $5\times 10^4\leq n\leq 5\times 10^5$ $(\Delta n =10^3)$ and used, as a model, $$x_n=\frac{\pi }{a n+\frac{b}{\log (n)}}$$ $$\begin{array}{clclclclc} \text{} & \text{Estimate} & \text{Standard Error} & \text{Confidence Interval} \\ a & +2.00000 & 6.02 \times 10^{-9}& \{+2.00000,+2.00000\} \\ b & -1.95234 & 0.00545 & \{-1.96306,-1.94162\} \\ \end{array}$$ I also considered separately $y_n=\frac{2n_n}{2n x_n-\pi}$ for which $$\begin{array}{clclclclc} \text{} & \text{Estimate} & \text{Standard Error} & \text{Confidence Interval} \\ c & 1.02127 & 0.00008 & \{1.02112,1.02142\} \\ \end{array}$$ This confirms the previous fit. In order to double check, I also minimized the sum of the squares of relative errors : this lead to the same result $(a=2.00000,b=-1.96355)$.
I also consider (as Did asked after Daniel Fischer's answer) $$x_n\approx\frac{\pi}{2n}\left(1+\frac{c}{n\log n}\right)$$ for which was obtained $$\begin{array}{clclclclc} \text{} & \text{Estimate} & \text{Standard Error} & \text{Confidence Interval} \\ c & 0.97571 & 0.00123 & \{0.97329,0.97813\} \end{array}$$
Now, considering the two models $$x_n^{(1)}=\frac{\pi }{2 \left(n-\frac{1}{\log (n)}\right)} \qquad x_n^{(2)}=\frac{\pi}{2n}\left(1+\frac{1}{n\log n}\right)$$ over the set of data points the sum of squared errors are respectively $1.24\times 10^{-20}$ and $6.94\times 10^{-23}$.
A friend of mine challenged me asking to show the iterations of the root finding method for $n=10^6$. Here they are
$$x_0= 1.570796440492926 \times 10^{-6}$$ $$x_1= 1.570796438443414 \times 10^{-6} $$ $$x_2= 1.570796438479340 \times 10^{-6} $$ $$x_3= 1.570796438479351 \times 10^{-6}$$
Many thanks to all of you for your contributions.