# Jumps in roots finding solutions

7 views (last 30 days)

Show older comments

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

##### 2 Comments

Dyuman Joshi
on 7 Feb 2024

What does "artifical jump" mean in this context?

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

I highly doubt so.

### Accepted Answer

Matt J
on 7 Feb 2024

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

##### 4 Comments

Walter Roberson
on 7 Feb 2024

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.

### More Answers (3)

Steven Lord
on 7 Feb 2024

##### 2 Comments

Steven Lord
on 7 Feb 2024

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

p1 = poly([1 2 3])

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

p2 = poly([2 3 1])

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

isequal(p1, p2)

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

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)

sort(r)

Walter Roberson
on 7 Feb 2024

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.

##### 1 Comment

Walter Roberson
on 7 Feb 2024

syms x t

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

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

tiledlayout('flow')

nexttile();

fplot(real(sol), [-4 2]); title('real')

nexttile();

fplot(imag(sol), [-4 2]); title('imag')

John D'Errico
on 7 Feb 2024

Edited: John D'Errico
on 7 Feb 2024

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

##### 1 Comment

Matt J
on 8 Feb 2024

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.

### See Also

### Categories

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!