Way to automatically distinguish branch of roots in parametric sweep?

Specific Domains Modelling & Simulations
1

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?

48%20AM
48%20AM1103×727 37.8 KB

Thanks!

2

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

3

Sorry, dropping it into Numerics

4

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
5

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

2 Likes
6

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

7

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.

8

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):

roots
roots6000×4500 1.3 MB

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]