5
$\begingroup$

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?

$\endgroup$
6
$\begingroup$

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] &)]

Mathematica graphics

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.