Jumps in roots finding solutions

Question

Carola Forlini 2024-2-7

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2079541-jumps-in-roots-finding-solutions

评论： 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 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 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).

请先登录，再进行评论。

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

Answer 1

Matt J 2024-2-7

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2079541-jumps-in-roots-finding-solutions#answer_1404641

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

4 个评论
显示 2更早的评论隐藏 2更早的评论

Carola Forlini 2024-2-7

Thank you for the clarification and the help!

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.

请先登录，再进行评论。

Answer 2

Steven Lord 2024-2-7

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2079541-jumps-in-roots-finding-solutions#answer_1404646

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 2024-2-7

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

Steven Lord 2024-2-7

在 MATLAB Online 中打开

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 2024-2-7

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2079541-jumps-in-roots-finding-solutions#answer_1404656

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 个评论
显示 -1更早的评论隐藏 -1更早的评论

Walter Roberson 2024-2-7

在 MATLAB Online 中打开

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 2024-2-7

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2079541-jumps-in-roots-finding-solutions#answer_1404726

编辑：John D'Errico 2024-2-7

在 MATLAB Online 中打开

eigenshuffle.m

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