Just a tangent-ing series

Hints:

It is quite obvious that $x_k$ are the positive roots of the equation $\tan x= x$

Using some coordinate geometry we can find that the locus of points of contact is $x^2-y^2=x^2y^2$

Using this information we see that we need to find $\sum_{k=1}^{\infty} \frac {x_k^2+1}{x_k^2(x_k^2+2)}=\frac 12\left( \sum_{k=1}^{\infty} \frac {1}{x_k^2} +\sum_{k=1}^{\infty} \frac {1}{x_k^2 +2}\right)$

Now writing power series of $\frac {\sin x-x\cos x}{x^3}$ about $0$ will give a polynomial with positive roots of $\tan x=x$ . Using Vieta's formula we see that $\sum_{k=1}^{\infty} \frac {1}{x_k^2}=\frac {\text{ -Coefficient of x²}}{\text {Constant term}}=\frac {1}{10}$

Also changing $x\rightarrow \sqrt {x-2}$ in the power series and using the same method above gives $\sum_{k=1}^{\infty} \frac {1}{x_k^2+2}=\frac {\text{ -Coefficient of x}}{\text {Constant term}}=\frac {5\sqrt 2\sinh(\sqrt 2)-6\cosh(\sqrt 2)}{4(2\cosh(\sqrt 2)-\sqrt 2\sinh(\sqrt 2))}$

You could also also check out another method used by robjohn on My question here

For the tangent to the curve $y=f(x)$ at $(t,f(t))$ to pass through $O$ , we need

$\begin{aligned} f'(t)=\frac{f(t)}{t} \end{aligned}$

In this case, this becomes $t=\tan t$ . The $x_k$ are the positive roots of this equation. It didn't seem likely that I would get anywhere with this analytically, so I proceeded numerically from here.

My idea was to use Newton's method to solve $t-\tan t=0$ in the intervals $(k\pi,(k+1)\pi)$ to work out the $x_k$ , which would then give $a_k$ . However, the gradients involved in the $\tan$ formulation are very large near the roots, so this iteration scheme is not numerically stable.

I therefore set up Newton's method for the (equivalent) equation $t\cos t - \sin t=0$ .

This gives as a scheme for iteration

$\begin{aligned} t_0 &= k\pi+\frac{\pi}{2} \\ t_{n+1} &= t_n + \frac{t_n\cos {t_n} - \sin {t_n}}{t_n\sin {t_n}} = t_n - \frac{1}{t_n}+\cot{t_n} \end{aligned}$

which converges rapidly to $x_k$ .

The first few roots are below:

We can then calculate the corresponding values of $a_k$ . If we let $S_N=\sum_{k=1}^N \frac{1}{a_k^2}$ (so the answer we want is $S_{\infty}$ ), we have the following partial sums:

We could stop here, but the partial sums converge quite slowly in this case. It's worth noting that we can approximate the error term $S_\infty-S_N$ as well: note that for large $k$ , the tangents occur nearer and nearer to the peaks and troughs of the sine curve, so that $x_k\approx k\pi+\frac{\pi}{2}$ and $y_k\approx\pm1$ and $a_k^2\approx 1+(k\pi+\frac{\pi}{2})^2$

If we define

$\begin{aligned} T_N=\sum_{k=0}^N \frac{1}{1+(k\pi+\frac{\pi}{2})^2} \end{aligned}$

(note the changed lower limit of the sum - to be explained in a moment) then $S_\infty-S_N \approx T_\infty-T_N$

The reason for the change of lower limit is that with this definition, we can use the (remarkable) identity $T_\infty=\frac{\tanh {1}}{2}$ (I can't claim credit for this - thank you, Wolfram|Alpha!)

So $S^*_N = S_N+\frac{\tanh {1}}{2}-T_N$ should be an even better approximation to $S_\infty$ . Indeed, we get

which clearly converges rapidly, giving the answer $\boxed{9737}$ .