Way to automatically distinguish branch of roots in parametric sweep?

Let’s say you have some nonlinear function that you plug into a root solver to get all the roots

// where the maximum number of real roots possible <= N (i.e. N = 3)

Then you sweep this function over M other values to get an array of M root lists (each of size 0 to N).

Is there a way to automatically distinguish the families of roots (i.e. assuming there are at most N families)


Through pictures,

  • can you create separate lines for the two families below
  • so you don’t have a saw-tooth looking function?

Thanks!

Is this a julia question or should it go into one of the domain categories?

Sorry, dropping it into Numerics

If I understand correctly, you are working on a parametric rootfinding problem.

In general, given some regularity assumptions, once you have the roots for a particular parameter you can use homotopy continuation methods to find the “same” root for other parameters. I am not aware of a Julia package, but writing one is on my summer todo list.

3 Likes

You are asking about “numerical continuation”. The standard software is Auto, which desperately needs rewriting in Julia…

2 Likes

To clarify, the name of the software is “Auto”, principal author Eusebius Doedel.

Once you have one solution, change the parameter a bit and look for the solution that is close. In that way you trace that “branch” of solutions. Then use the other one and trace its branch.

Thanks for the help! I solved my simple example with this code:

// where cur_array has dimensions: length(x_list) by max_roots_count


Before explaining it quick, here are the pretty pictures (where markers denote root id):


Now the quick explanation.

  • we find roots using find_zeros in Roots.jl (you’d have to use my branch that allows the function to go complex in some regions)
  • we piece together all the points that display the max root count (obviously a big assumption!)
  • we fit the curves using a quadratic polynomial // i.e. a parabola :stuck_out_tongue:
  • each new point is added to the closest polynomial with some pretty loose discretion

Obviously this is too simple, but it is a step in the right direction

// with no optimization (premature or otherwise)

edit: one way to optimize it for sure would be tracking roots and only doing find_zeros on mystery areas – i.e. [a, min_root) and (max_root, b]