# Finding smooth root sets to a transcendental equation

I'm trying to solve for the roots of an equation like

$f(a, n) = \frac{n^2}{n+1} e^{i a n }$

for some list of complex $f$ values. For each $f$ value I know the (real) value of $a$, and I'd like to find an $n$ (complex) that satisfies. The catch is, I would like output list of $n$ values to have "smooth" real and imaginary parts. That is, plotting $n(a)$ I would see no "jumps" from one branch of solutions to another.

I've seen a number of related posts around here already about continuous phase unwrapping, collecting transcendental equations' roots, and plotting branching functions. But here I'd like to ask more specifically about finding the roots in a smooth, possibly iterative way.

My attempt at a solution is to search for a root starting at previous root, making it more difficult for FindRoot[] to jump away:

 smoothFindRoot[
list_,
list2_,
function_,
initialApproximation_] := First@Last@Reap@
For[i = 1, i <= Length@list, i++,
If[i == 1,
rootApproximation = initialApproximation,
rootApproximation = Round[z, .1];
];
Sow[z = x /.
FindRoot[Evaluate[list[[i]]] == function[list2[[i]]],
{x, rootApproximation}]]]


But doing something iteratively like this is new to me in Mathematica, and I feel there must be a better way. For one thing, it defines global variables. Can one use Module[] or Fold[] or some other functions to make this better?

• Commented Sep 11, 2016 at 21:03

You can use NDSolve to trace a solution by differentiating the equation:

eqn = f == n^2/(n + 1) Exp[I a n];

{sol} = NDSolve[{D[eqn /. n -> n[a], a], n[1] == 1}, n, {a, -10, 10},
InterpolationOrder -> All];


Or by also using the "Projection" method, which when coupled with InterpolationOrder -> All, gives a better approximation. The code below gives a single-precision approximation. Uncomment the options to get a double-precision approximation.

Block[{f = n^2/(n + 1) Exp[I a n] /. {n -> 1, a -> 1}},
{sol} =
NDSolve[{D[eqn /. n -> n[a], a], n[1] == 1}, n, {a, -10, 10},
Method -> {"Projection",
"Invariants" -> {eqn /. n -> n[a] /. Equal -> Subtract}},
InterpolationOrder -> All(*, AccuracyGoal -> 24,
PrecisionGoal -> 24, WorkingPrecision -> 32*)]
];


Check the residual of the equation:

Block[{f = E^I/2},
Table[eqn /. n -> n[a] /. Equal -> Subtract /. sol, {a, -10, 10, 0.1}]
] // Abs // Max
(*  2.73044*10^-9  *)


Visualize the trace of the root:

ParametricPlot[ReIm[n[a]] /. sol, {a, -10, 10},
Epilog -> {PointSize[Medium],
Table[Tooltip[Point[ReIm[n[a]] /. sol], a], {a, -10, 10}]},
ColorFunction -> (ColorData["Rainbow"][#3] &)]