In general, are there any clever tricks to help find the roots of a finite Fourier series? Presumably there aren't analytic methods, but can we use the fact that our function is a finite Fourier series to facilitate numerical methods?
The hardest part of numerical root finding is bracketing the root, so is there a way we can do that easily? I know for a general function it's a hard problem.
The short answer to your question is "yes", you can exploit the structure of the Fourier-series to reduce your rootfinding problem to a matrix eigenvalue problem, which can then be solved efficiently using QR decomposition. None of this is my own, but I'm paraphrasing a paper by John P. Boyd at University of Michigan from 2006, entitled "Computing the zeros, maxima and inflection points of Chebyshev, Legendre and Fourier series: solving transcendental equations by spectral interpolation and polynomial rootfinding". Let a finite Fourier-series with $2N+1$ terms be given by
$$ f_N(t) = \sum_{j=0}^N a_j\cos(jt) + \sum_{j=1}^N b_j\sin(jt) $$
Define the following $2N+1$ coefficients:
$$ h_k = \left\{\begin{array}{ll} a_{N-k}+i b_{N-k}, & k=0,1,\ldots,(N-1)\\ 2 a_0 & k=N \\ a_{k-N} - i b_{k-N} & k=N+1,N+2,\ldots, 2N \end{array}\right. $$
And define the following $2N\times2N$ matrix $\bf B$ with entries $B_{jk}$ at indicies $j,k$ as $$ B_{jk} = \left\{ \begin{array}{ll} \delta_{j,k-1}, & j=1,2,\ldots,(2N-1)\\ - \frac{h_{k-1}}{a_N - i b_{N}} & j=2N \end{array}\right. $$
The delta above is the Kronecker delta that is one when its two arguments are the same, zero otherwise. Let the matrix $\bf B$ have eigenvalues given by $z_k$, then the roots of $f_N(t)$ are given by
$$ t_k = -i \log(z_k) $$
In general $z_k$ will be complex or negative, so we mean the complex logarithm (this nicely gives that each root is also periodic, the expected behavior of a periodic function)
$$ \log(z) = \log(|z|) + i (\arg(z)+2\pi m) \;\;\text{for integer }m $$
For a final formula for the periodic roots
$$ t_k = \arg(z_k)+2\pi m-i \log(|z_k|),\;\;k=1,2,\ldots,2N,\;\;m\text{ an integer} $$
So the numerical task of Fourier series rootfinding is reduced to basically solving for eigenvalues of the matrix $\bf B$.
