roots() function does not give exact roots

33 次查看（过去 30 天）
Mustafa Yasar
Mustafa Yasar 2020-12-15
回答： John D'Errico 2020-12-16
When I write this code:
roots( [1 -3 3 -1] );
I expect to get three 1 roots. However Matlab gives this result:
1.000006571943641, 0.999996714028179, 0.999996714028179.
What should I do to get exact 1's.
  6 个评论
显示 4更早的评论
Bruno Luong
Bruno Luong 2020-12-15
"vandermonde"
Use "\"????
Nah you seem to mix between POLYIT Wlater
The roots are eigen values of companion matrix, and it does not have anything remotely with vandermonde and "\"
Walter Roberson
Walter Roberson 2020-12-15
ah, I must have been up all night again... I need more sleep :(

请先登录，再进行评论。

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

回答（3 个）

John D'Errico
John D'Errico 2020-12-16
Ran in:
What it seems is people have not explained the "cause" of the problem as seen.
format long g
roots( [1 -3 3 -1] )
ans = 3×1
1.00000657194364 + 0i 0.999996714028179 + 5.69145454681503e-06i 0.999996714028179 - 5.69145454681503e-06i
Remember that roots is a numerical rootfinder. It does not understand that this implicit cubic polynomial, described by the integer coefficients of the polynomial, has a nice set of triplicated integer roots.
Anyway, the function roots is a tricky one, in that it smartly converts your polynomial into a matrix, known in mathematics as the companion matrix.
https://en.wikipedia.org/wiki/Companion_matrix#:~:text=In%20this%20sense%2C%20the%20matrix%20C%28p%29%20is%20the,companion%20matrix%20over%20K%20of%20its%20characteristic%20polynomial.
That matrix is simple to construct, and then eig can do its stuff. It assume the polynomial has leading coefficient of the highest order term as 1. Clearly if it is some other value, you can divide all coefficients by that value, but this would not change the roots.
If you assume the polynomial has coefficients
[1, c_(n-1), ..., c_2, c_1, c_0]
representing the polynomial
x^n + c_(n-1)*x^(n-1) + ... + c_2*x^2 + c_1*x + c_0
then we construct the companion matrix as an nxn matrix. For the cubic polynomial, it looks like:
syms c2 c1 c0 x
M = [0 0 -c0;
1 0 -c1;
0 1 -c2]
M = 
Does this have the desired set of roots?
det(x*eye(3) - M)
ans = 
So you should see we can go either way. We can convert a roots problem into an eigenvalue problem. But this is really all that roots does. It calls eig on the companion matrix.
Anyway, where does all of this leave us? Oh, yes. Now we need to discuss root finding on problems with multiple coincident roots. It really does not matter how you will try to solve the problem, multiply replicate roots are difficult. Why? consider the linear polynomial
y(x) = x - x0
this is just a straight line. It passes though the point on the x-axis at location (x0,0) in the (x,y) plane. All is good in the world. Suppose I make an error of dx in the extimated root for y? So consider that I will evaluate the function at the point x0+dx, for some small value of dx? Then we would have
y(x0 + dx) = x0 + dx - x0 = dx
so the difference in y is LINEARLY related to the error in our estimate of the root.
But now consider a quadratic polynomial, with a DOUBLE root at x0? This polynomial would be written as:
z(x) = x^2 - 2*x0*x + x0^2
surely you agree this has a DOUBLE root at x0? But evaluate z at the point x0+dx.
syms x x0 dx
z(x) = x^2 - 2*x0*x + x0^2
z(x) = 
At x==x0, we would expect to get exactly zero.
expand(z(x0))
ans = 
0
Of course, remember that if our computations are done in double precision arithmetic, we might have some fuzz down in the least significant bits.
expand(z(x0 + dx))
ans = 
So here, a tiny perturbation in the value of x from x0, by an amount dx, we see that z(x+x0) is now only dx^2.
Effectively, we see a nice parabolic shape as we move away from x0. But if dx is small, then dx^2 is quadratically small. The curve is actually quite flat down there. For triple roots, things get worse. Now we see the curve gets even more flat. And now our ability to resolve the shape of the curve near zero is even worse.
Consider the cubic polynomial defined by the coefficients:
P3 = [1 -3 3 -1];
We know this has a triple root at x==1.
fplot(@(x) polyval(P3,x),[0.99,1.01])
Do you see the y axis is quite compressed there? While the function looks cubic in nature, see the y axis is multiplied by 1e-6. So within limits of +/- 0.01 around the known root, the function itself varies only by 1.e-6.
fplot(@(x) polyval(P3,x),1 + [-1e-6,1e-6])
Now when I tightened the limits on that plot, see that P3(x) now starts to look itself like almost random fuzz around zero. Effectively, we cannot resolve triplicate roots to closer than roughly nthroot(eps(x0),3). Here, that becomes approximately
nthroot(eps(1),3)
ans =
6.05545445239334e-06
We see that translate into what roots does for this polynomial.
>> roots([1 -3 3 -1])
ans =
1 + 0i
1 + 5.6915e-06i
1 - 5.6915e-06i
Roots simply gets lost in the fuzz down there. And THAT is the ultimate cause of what you saw.
Can you fix it? Not simply done when working in double precision arithmetic. Of course, if you use a symbolic polynomial, then syms will handle it. But roots must fail here to give you more precision than the cube root of eps.
KSSV
KSSV 2020-12-15
Ran in:
p = [1 -3 3 -1] ;
poly = poly2sym(p)
poly = 
r = solve(poly==0)
r = 
% using roots
r = roots(p) ;
real(r)
ans = 3×1
1.0000 1.0000 1.0000
  7 个评论
显示 5更早的评论
Walter Roberson
Walter Roberson 2020-12-16
Ran in:
KSSV: what release and platform are you using?
R2020b on Windows
>> roots([1 -3 3 -1])
ans =
1.00000657194364 + 0i
0.999996714028179 + 5.69145454681503e-06i
0.999996714028179 - 5.69145454681503e-06i
And on Linux:
format long g
roots([1 -3 3 -1])
ans = 3×1
1.00000657194364 + 0i 0.999996714028179 + 5.69145454681503e-06i 0.999996714028179 - 5.69145454681503e-06i
●
KSSV
KSSV 2020-12-16
Walter Roberson: I have used Octave. I have no acess to MATLAB. Agree, this could be due to version. Thought to answer the same to OP but looking at his comments I didn't feel like answering him.

请先登录，再进行评论。

Matt J
Matt J 2020-12-15
编辑：Matt J 2020-12-15
Ran in:
double(roots( sym([1 -3 3 -1] ))) == [1;1;1] %all true
ans = 3x1 logical array
1 1 1

类别

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

标签

产品


版本

R2020b

Community Treasure Hunt

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

Start Hunting!

Translated by