Algorithm for calculating real, positive roots of transcendental equation involving tangens

Question

Crank ("The mathematics of diffusion", 2nd editon, 1975, p.57) describes a diffusion modelling algorithm which relies on the non-zero positive roots of

$$ \tan{q_n} = -\alpha\cdot q_n$$

Typical values I use for $\alpha : 0.1, 1, 10, 15, 25, 50,... $

To achieve the necessary accuracy for small diffusion coefficients (e. g. $10^{-16} $) I need to calculate several hundred roots. I currently use Mathematica to calculate the roots, but this is a quite painful procedure. Therefore I would like to write a C++ program which will help me calculating the roots. For that I would need an algorithm to solve above equation. Any help with that would be appreciated.

Answer 1

You are looking for the zero's of function $$f(x)=\tan(x)+\alpha x$$ They will be closer and closer to $(2n+1)\frac \pi 2$. To remove the unpleasant discontinuities, consider instead $$g(x)=\sin(x)+\alpha x \cos(x)$$ and use Taylor series around $x=(2n+1)\frac \pi 2$ $$g(x)=\frac 12\sum_{k=0}^\infty(-1)^n\frac{ 2 (\alpha k+1) \cos \left(\frac{\pi k}{2}\right)-\pi \alpha (2 n+1) \sin \left(\frac{\pi k}{2}\right)}{ k!}\Big[x- (2 n+1)\frac \pi 2\Big]^k$$ Truncate to some order and use series reversion to obtain fot the $n^{\text{th}}$ root $$x_{(n)}= t+\sum_{m=0}^p (-1)^{m}\frac {c_m} {(\alpha t)^{2m+1}}\qquad \text{where} \qquad t=(2n+1)\frac \pi 2$$

The first coefficients $c_m$ are $$\left( \begin{array}{cc} m & c_m \\ 0 & 1 \\ 1 & \alpha +\frac{1}{3} \\ 2 & 2 \alpha ^2+\frac{4 \alpha }{3}+\frac{1}{5} \\ 3 & 5 \alpha ^3+5 \alpha ^2+\frac{23 \alpha }{15}+\frac{1}{7} \\ 4 & 14 \alpha ^4+\frac{56 \alpha ^3}{3}+\frac{392 \alpha ^2}{45}+\frac{176 \alpha }{105}+\frac{1}{9} \\ 5 & 42 \alpha ^5+70 \alpha ^4+44 \alpha ^3+\frac{818 \alpha ^2}{63}+\frac{563 \alpha }{315}+\frac{1}{11} \\ 6 & 132 \alpha ^6+264 \alpha ^5+209 \alpha ^4+\frac{15796 \alpha ^3}{189}+\frac{3102 \alpha ^2}{175}+\frac{6508 \alpha }{3465}+\frac{61}{1024} \\ 7 & 429 \alpha ^7+1001 \alpha ^6+\frac{43043 \alpha ^5}{45}+\frac{65351 \alpha ^4}{135}+\frac{282841 \alpha ^3}{2025}+\frac{1534247 \alpha ^2}{57600}+\frac{141583 \alpha }{23040}+\frac{18811}{15360} \end{array} \right)$$

Trying for $\alpha=10$ with $p=7$

$$\left( \begin{array}{ccc} 0 &\color{red}{1.63199452}646140205048274685183 & 1.63199452721480006371688983888 \\ 1 & \color{red}{4.733511802356786}19428997112737 & 4.73351180235678620084907851480 \\ 2 & \color{red}{7.86669277156157419748}480845369 & 7.86669277156157419748592987705 \\ 3 & \color{red}{11.004661096033518512264}5971500 & 11.0046610960335185122646008382 \\ 4 & \color{red}{14.1442368407184333320452158}130 & 14.1442368407184333320452158645 \\ 5 & \color{red}{17.28454504550461827956001048}46 & 17.2845450455046182795600104863 \\ 6 & \color{red}{20.425248110576069042708349747}2 & 20.4252481105760690427083497473 \\ 7 & \color{red}{23.5661882440696287891706461263} & 23.5661882440696287891706461263 \end{array} \right)$$

Edit

If you are not looking for too much accuracy, rewrite

$$x_{(n)}= t+\sum_{m=0}^p (-1)^{m}\frac {c_m} {(\alpha t)^{2m+1}}=t+\frac 1 {\alpha t}\sum_{m=0}^p (-1)^{m} c_m \,y^m$$ with $y=\frac 1{(\alpha t)^2} $ and use the simplest Padé approximant $$\sum_{m=0}^p (-1)^{m} c_m \,y^m=\frac{c_0 c_1+\left(c_0 c_2-c_1^2\right) y } {c_1+c_2y }+O\left(y^3\right)$$

For $\alpha=10$ as above, the solutions will be $$\{1.631998,4.733512,7.866693,11.00466,14.14424,17.28455,20.42525,23.56619\}$$