Jumps in roots finding solutions

Hi,

I am using roots to find the solutions of a cubic polynomial equation.

When I plot the solutions (real and imaginary part) I notice some aritifcial jump in the solution which are not possible. Could this be an issue of the roots algorithm?

Thank you

[Dyuman Joshi]

What does "artifical jump" mean in this context?

"Could this be an issue of the roots algorithm?"

I highly doubt so.

[Carola Forlini]

The jump you see around k=2.3 do not make sense from physics point for the problem I am solving (travelling waves).

Answer 1

[Matt J]

I think you'll see continuity if you sort the roots.

[Carola Forlini]

Yes, I tried and I now see continuity.

Why? Am I not changing the solution by sorting the roots?

[Matt J]

No, you are not changing the solution. A set of roots will correspond to the same polynomial coefficients no matter how you order them:

p=poly([1,2,3])

p = 1×4

1 -6 11 -6

p=poly([3,2,1])

p = 1×4

1 -6 11 -6

[Carola Forlini]

Thank you for the clarification and the help!

[Walter Roberson]

Am I not changing the solution by sorting the roots?

Yes? No?

You haven't defined what solution you are going for.

If you have cubics in a parameter, then as the parameter evolves, the roots evolve. There is no inherent ordering of the roots: roots() imposes an (undocumented) ordering on them. It is common in systems that two of the roots to meet. At that point, what was the "upper" branch typically becomes the "lower" branch. If the coefficients are sorted somehow then their identities effectively swap.

The way to avoid this is to use the symbolic toolbox, solve() the equation with 'maxdegree', 3, getting out three solutions, and then plot the three solutions over the parameter.

Answer 2

[Steven Lord]

The documentation of the roots function doesn't state anything about a particular order in which the roots are returned. What happens if you sort the real parts of the roots of your polynomial for each value of k before deciding which one is "root 2" and which one is "root 3"? I'd bet you no longer see that jump.

[Carola Forlini]

Yes, I do not see the jump anymore but I still don't really understand why.

[Steven Lord]

Let's say I create a polynomial with roots 1, 2, and 3.

p1 = poly([1 2 3])

p1 = 1×4

1 -6 11 -6

Now let's say I create a polynomial with roots 2, 3, and 1.

p2 = poly([2 3 1])

p2 = 1×4

1 -6 11 -6

The order doesn't matter; I get the same polynomial either way.

isequal(p1, p2)

ans = logical
   1

Now let's choose one of the polynomials at random. I'll use randi(2) to "flip a coin".

if randi(2) == 1
    p = p1;
else
    p = p2;
end
p

p = 1×4

1 -6 11 -6

Can you tell, just by looking at p, whether I used [1 2 3] or [2 3 1] as the input to poly to create it? No.

So which is the "second root" of p? Is it 2 or 3? Or could it be 1?

Now if I sort the roots of p I could say something about ordering. Note that r is neither [1 2 3] nor [2 3 1]. That is not a bug.

r = roots(p)

r = 3×1

3.0000 2.0000 1.0000

sort(r)

ans = 3×1

1.0000 2.0000 3.0000

Answer 3

[Walter Roberson]

The order of solutions produced by roots() is not documented.

If you need a consistent order then use the Symbolic Toolbox, solve() the equation with 'Maxdegree', 3, and subs() specific numeric values for the free variable.

[Walter Roberson]

syms x t
y = x^3 + (2+t)*x^2 - t^2*x + (1+t)

y = 

sol = solve(y, x, 'maxdegree', 3)

sol = 

tiledlayout('flow')
nexttile();
fplot(real(sol), [-4 2]); title('real')
nexttile();
fplot(imag(sol), [-4 2]); title('imag')

Answer 4

[John D'Errico]

Lets spend some time to understand what happens. Consider the simple polynomial:

y = x^5 + 3*x^4 + (2+t)*x^2 - t^2*x + (1+t)

I just made it up, so I have no clue as to how the roots behave as a function of the parameter t. Remember though, those roots are generated by roots, which actually converts the problem to an eigenvalue problem.

P = @(t) [1,3,0,2+t,t^2,1+t];
T = linspace(-5,5,1001);
R = NaN(numel(T),5);
for i = 1:numel(T)
  R(i,:) = roots(P(T(i))).';
end

Now, plot them in the complex plane. Plot plots the real part on the x axis, and the imaginary part on the y axis.

plot(R,'-')

In this first plot, there are 5 trajectories plotted. The problem is, the blue curve, for example, seems to jump around.

We can better visualize the 5 root trajectories below.

plot(R,'.')

Note that those points are not connected in any specific sequence. They are just plotted as disconnected dots, so our brain can see what should happen. In fact, looking at the first plot, you can see the sequence of those points has them jumping from one trajectory to another.

It looks like there is one root that is always purely real. Another on the left, follows a nearly circular arc. Then there are three other curved arcs, some of which cross each other. Now go back and look at the blue line. Do you see that it appears to hop from one of those root trajectories to another? The problem is that the function roots does not care in which order the roots are generated. It cannot do so. So they just fall out in a pretty arbitrary order.

Can you fix that, by sorting the roots in order of their real part? NO. Clearly that would make as much of a mess as the default order.

The problem comes down to a question of bifurcation. As the coefficients of the polynomial change, at some point, those paths will cross. When the paths cross, there will be some value of t where a pair of roots temporarily becomes a double root, or even possibly a triple root in some cases. At that point, if roots were following that path, roots would need to know which branch to follow for each root. It cannot easily do that. Remember, roots does not see the sequence of polynomials. It sees only one polynomial at a time, and this is the crux of the problem. Can you resolve that?

One idea is to convert the polynomial problem into an eigenvalue problem. My eigenshuffle code can resequence the roots, because it also uses the eigenvectors. Thus, given any polynomial, we can create the companion matrix.

C = NaN(5,5,numel(T));
for i = 1:numel(T)
  C(:,:,i) = compan(P(T(i)));
end

[Vseq,Dseq] = eigenshuffle(C);

plot(Dseq.','-')

Sadly, it has improved, but even so, it got some of those trajectories wrong. Oh well. At least I want you to understand the problem. (Darn, I need to revisit eigenshuffle.)

[Matt J]

Sadly, it has improved, but even so, it got some of those trajectories wrong

By what criteria are they "wrong"? They all appear to be continuous.