Jumps in roots finding solutions

7 views (last 30 days)
Carola Forlini
Carola Forlini on 7 Feb 2024
Commented: Matt J on 8 Feb 2024
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
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.
Carola Forlini
Carola Forlini on 7 Feb 2024
The jump you see around k=2.3 do not make sense from physics point for the problem I am solving (travelling waves).

Sign in to comment.

Accepted Answer

Matt J
Matt J on 7 Feb 2024
I think you'll see continuity if you sort the roots.
  4 Comments
Carola Forlini
Carola Forlini on 7 Feb 2024
Thank you for the clarification and the help!
Walter Roberson
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.

Sign in to comment.

More Answers (3)

Steven Lord
Steven Lord on 7 Feb 2024
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.
  2 Comments
Carola Forlini
Carola Forlini on 7 Feb 2024
Yes, I do not see the jump anymore but I still don't really understand why.
Steven Lord
Steven Lord on 7 Feb 2024
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

Sign in to comment.


Walter Roberson
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
Walter Roberson on 7 Feb 2024
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')

Sign in to comment.


John D'Errico
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
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.

Sign in to comment.

Categories

Find more on Data Preprocessing in Help Center and File Exchange

Community Treasure Hunt

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

Start Hunting!