About the first positive root of $\sum_{k=1}^n\tan(kx)=0$

Did is spot on. I'm not surprised, because he usually is, but I had to calculate to be sure that the used approximations were good enough to yield the right result.

For $n \in \mathbb{N}\setminus \{0\}$, let

$$f_n(x) = \sum_{k = 1}^n \tan (kx)$$

and $x_n$ the smallest positive zero of $f_n$. Then $f_n$ is an entire $\pi$-periodic meromorphic function with simple poles at the points of

$$P = \biggl\{ \frac{2m+1}{2k}\pi : m \in \mathbb{Z},\, 1 \leqslant k \leqslant n\biggr\}.$$

On $\mathbb{R}\setminus P$, $f_n$ is real valued, and $f_n' > 0$, so $f_n$ maps each interval between consecutive poles diffeomorphically onto $\mathbb{R}$. The case $n = 1$ is trivial ($f_1 = \tan$), and for $n \geqslant 2$ the two smallest positive poles are at $\frac{\pi}{2n}$ and $\frac{\pi}{2(n-1)}$. Since $f_n(0) = 0$ and $f_n' > 0$ on $\mathbb{R}\setminus P$, we have $f_n(x) > 0$ for $0 < x < \frac{\pi}{2n}$, so it follows that $x_n \in \bigl(\frac{\pi}{2n},\frac{\pi}{2(n-1)}\bigr)$, and $f_n$ has - for $n \geqslant 2$ - no other zeros in that interval. It is easily seen that $x_2 = \frac{\pi}{3}$, so in the following we assume $n \geqslant 3$. Since

\begin{align} f_n \biggl(\frac{\pi}{2n-1}\biggr) &= f_{n-2}\biggl(\frac{\pi}{2n-1}\biggr) + \tan \frac{\pi(n-1)}{2n-1} + \tan \frac{\pi n}{2n-1}\\ &= f_{n-2}\biggl(\frac{\pi}{2n-1}\biggr) + \tan \biggl(\frac{\pi}{2} - \frac{\pi}{4n-2}\biggr) + \tan \biggl(\frac{\pi}{2} + \frac{\pi}{4n-2}\biggr)\\ &= f_{n-2}\biggl(\frac{\pi}{2n-1}\biggr)\\ &> 0, \end{align}

it follows that $x_n \in \bigl(\frac{\pi}{2n}, \frac{\pi}{2n-1}\bigr)$. Write $x_n = \frac{1}{n}\bigl(\frac{\pi}{2} + \delta_n\bigr)$. Then $0 < \delta_n < \frac{\pi}{4n-2}$, in particular $\delta_n < \frac{3}{5} x_n$.

We use $\tan \bigl(\frac{\pi}{2} - x\bigr) = \cot x$, and

$$\frac{1}{x} - \frac{x}{2} < \cot x < \frac{1}{x} - \frac{x}{3}\tag{1}$$

for $0 < x < \frac{\pi}{2}$ in our calculations. First we have

$$-\tan (n x_n) = -\tan \biggl(\frac{\pi}{2} + \delta_n\biggr) = \cot \delta_n = \frac{1}{\delta_n} + O(\delta_n).\tag{2}$$

To determine the asymptotic behaviour of $\delta_n$, we next note that (since $f_n(x_n) = 0$)

\begin{align} -\tan (nx_n) &= \sum_{k = 1}^{n-1} \tan (k x_n)\\ &= \sum_{m = 1}^{n-1} \tan \bigl((n-m)x_n\bigr)\\ &= \sum_{m = 1}^{n-1} \tan \biggl(\frac{\pi}{2} - (mx_n - \delta_n)\biggr)\\ &= \sum_{m = 1}^{n-1} \cot (mx_n - \delta_n). \end{align}

The inequalities $(1)$ now yield

$$\sum_{m = 1}^{n-1} \frac{1}{m x_n - \delta_n} - \frac{1}{2} \sum_{m = 1}^{n-1} (m x_n - \delta_n) < -\tan (n x_n) < \sum_{m = 1}^{n-1} \frac{1}{m x_n - \delta_n} - \frac{1}{3} \sum_{m = 1}^{n-1} (m x_n - \delta_n).\tag{3}$$

We find

$$\sum_{m = 1}^{n-1} (m x_n - \delta_n) = \frac{n(n-1)}{2}x_n - (n-1)\delta_n = (n-1)\frac{\pi - 2\delta_n}{4}$$

and

\begin{align} \sum_{m = 1}^{n-1} \frac{1}{m x_n -\delta_n} &= \sum_{m = 1}^{n-1} \frac{1}{m x_n} + \frac{\delta_n}{x_n^2} \sum_{m = 1}^{n-1} \frac{1}{m\bigl(m - \frac{\delta_n}{x_n}\bigr)}\\ &= \frac{\log n + \gamma + O(n^{-1})}{x_n} + O(\delta_n x_n^{-2}). \end{align}

With $x_n^{-1} \sim \frac{2}{\pi} n$ it follows that

$$\frac{1}{\delta_n} \sim - \tan (n x_n) = \frac{2}{\pi}n \log n + O(n),$$

or

$$\delta_n = \frac{\pi}{2n\log n}\bigl( 1 + O\bigl((\log n)^{-1}\bigr)\bigr).\tag{4}$$

With some tedious work, we can get some bounds on the $O\bigl((\log n)^{-1}\bigr)$ term in $(4)$, but since $\frac{\pi}{12} < \frac{2\gamma}{\pi} < \frac{\pi}{8}$, what we have isn't sufficient to even determine whether it is $\Theta\bigl((\log n)^{-1}\bigr)$.

However, $(4)$ suffices to show that

$$\begin{split}x_n &= \frac{1}{n}\biggl(\frac{\pi}{2} + \delta_n\biggr) = \frac{\pi}{2n}\biggl(1 + \frac{1 + O\bigl((\log n)^{-1}\bigr)}{n\log n}\biggr)\\ &= \frac{\pi}{2n\bigl(1 - \frac{1 + O((\log n)^{-1})}{n\log n}\bigr)} = \frac{\pi}{2\bigl(n - (\log n)^{-1} + O((\log n)^{-2})\bigr)}.\end{split}\tag{5}$$

Concerning the difference between the empirical best-fit constants and the exact asymptotic values, recall Legendre's constant. The logarithm is a very slowly varying function, to eliminate effects of constants from empirical estimates, one may need very large numbers. However, it may be that $250000$ is large enough, and the difference between the best-fit and the exact values is caused by the larger deviation from the asymptotic behaviour for the smaller $n$. Try a best-fit for e.g. $200000 \leqslant n \leqslant 250000$ to see what that gives.

This is not an answer, just a long comment. When you deal with quantities that converge to $0$ you have to be careful with the sence of the approximation. Notice that if we take any two sequences $x_n,y_n$ such that $$ x_n,y_n \in [\frac{\pi}{2n},\frac{\pi}{2n - 1}]$$ then $|x_n - y_n| \leq \frac{\pi}{2n - 2} - \frac{\pi}{2n} = \frac{2\pi}{(2n - 2)(2n)} = \Theta(\frac{1}{n^2})$ hence $$\frac{|x_n - y_n|}{x_n} = \Theta(\frac1n),$$ so even the relative error goes to $0$. So it may seems that $\frac{\pi}{2n- 1}$ is a good approximation but in fact is as good as any other. So to define a good approximation we should require instead $$\frac{|x_n - y_n|}{x_n} = o(\frac1n).$$ With this definition you can see at least numericaly that $\frac{\pi}{2n}$ is the right approximation as Cave jonhson pointed out.