I am trying to solve a bunch of equations for the zeros of the derivative of an analytic function, and I would like to know if there exist methods that exploit this structure to provide better performance than the standard algorithms.
At the moment I am using Mathematica's FindRoot function, which I understand relies on Newton's method. (I iterate over quasirandom seeds as described here.) My problem, however, has additional, structure, of the form
$$f'(z)=0$$
for an analytic $f$, so maybe there is a better way to do this. So: are there methods that exploit this structure to provide better performance?
To be a bit more explicit and provide some context in case it's useful: I am writing a Mathematica package to solve the saddle-point equations for high-order harmonic generation (as explained e.g. here), which are of the form
$$
\left\{\begin{aligned}
(p(t,u)-A(u))^2+\gamma^2 & =0 \\
(p(t,u)-A(t))^2+\gamma^2 & = \omega,
\end{aligned}\right.
\tag{1}
$$
where $p(t,u)=\frac{1}{i\varepsilon+t-u}\int_u^tA(\tau)\mathrm d\tau$ for $\varepsilon$ a small, positive constant, $\gamma>0$ is fixed, and $\omega>0$ is a parameter. Here $A(t)$ is the vector potential of a laser field and can be a three-dimensional vector, but a simple example is $A(t)=A_0\sin(t)$. These equations can be seen as looking for the zeros of both partial derivatives of
$$
S(t,u)=\gamma^2(t-u)-\omega t+\int_u^t(p(t,u)-A(\tau))^2\mathrm d\tau.
$$
$S$ can be assumed to be known explicitly but in some cases it is relatively awkward (with LeafCounts in the >6000 range).
Some notes on the behaviour of this system:
- I am interested in a specific box in the complex plane for $t$ and $u$. In some regimes some roots can wander off the top and bottom of this box (i.e. to higher $|\mathrm{Im}(t)|$ than prescribed), and in those cases the roots would mostly be ignored anyway. It would be nice to have them as it simplifies the accounting, but they're not crucial.
- Some roots which are mostly not of interest can also wander off the sides of the box, and those would definitely get ignored anyway. However, I have no guarantees on the number of roots in my chosen box, and it takes some fiddling after I've got the entire curves w.r.t. $\omega$ to design a function $g(t_s,u_s)$ that will take a root $(t_s,u_s)$ and tell me whether it's a keeper or not.
- The regime $A_0^2\gg\gamma^2$ is of particular interest. In this case the roots will mostly have real $t$ roots, and as $\omega$ varies they will approach each other, have avoided crossings, and veer off into nonzero $\mathrm{Im}(t)$. When $A_0^2\gg\gamma^2$ these avoided crossings can be very tight and happen very quickly with respect to variations in $\omega$, but I'm OK with having to deal with those $\omega$ regions separately.
- At the moment I have an explicit expression for $S(t,u)$ which I differentiate symbolically with Mathematica, and this gives derivatives with
LeafCountroughly 6 to 8 times those of $S$. They're not stuff I'd like to simplify by hand. - I would appreciate methods which integrate well with Mathematica (either already implemented, existent as third-party tools, or relatively easy to implement), but if none are available I'm interested in any methods in this area.
- Sometimes I do have access to reasonable guesses for the roots (via e.g. a related, simpler $A(t)$, or using the results from the previous $\omega$) but I would rather avoid this or at least have methods which work even without such guesses.
- If a simpler model is desired, setting $\gamma=0$ can yield good estimates for $t$ (but not necessarily for $u$).
