Jumps in roots finding solutions

3 次查看（过去 30 天）
Carola Forlini
Carola Forlini 2024-2-7
评论： Matt J 2024-2-8
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 个评论
Dyuman Joshi
Dyuman Joshi 2024-2-7
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 2024-2-7
The jump you see around k=2.3 do not make sense from physics point for the problem I am solving (travelling waves).

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Matt J
Matt J 2024-2-7
I think you'll see continuity if you sort the roots.
  4 个评论
显示 2更早的评论
Carola Forlini
Carola Forlini 2024-2-7
Thank you for the clarification and the help!
Walter Roberson
Walter Roberson 2024-2-7
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.

请先登录，再进行评论。

更多回答（3 个）

Steven Lord
Steven Lord 2024-2-7
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 个评论
Carola Forlini
Carola Forlini 2024-2-7
Yes, I do not see the jump anymore but I still don't really understand why.
Steven Lord
Steven Lord 2024-2-7
Ran in:
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

请先登录，再进行评论。

Walter Roberson
Walter Roberson 2024-2-7
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 个评论
Walter Roberson
Walter Roberson 2024-2-7
Ran in:
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')

请先登录，再进行评论。

John D'Errico
John D'Errico 2024-2-7
编辑：John D'Errico 2024-2-7
Ran in:
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 个评论
Matt J
Matt J 2024-2-8
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.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Polynomials 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by