roots() function does not give exact roots

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.

[KALYAN ACHARJYA]

"What should I do to get exact 1's."

>> real(roots( [1 -3 3 -1] ))

[KALYAN ACHARJYA]

Mustafa Yasar Comments moved here

shortG did not solve my problem. Still I am getting the same numbers

[KALYAN ACHARJYA]

Yes, please wait, other experts will answer the question. I read the docs, still, the exact answer (satisfactory) answer is unknown to me.

[Walter Roberson]

I see what you mean, my tests agree with your tests.

roots() constructs a vandermonde companion matrix and uses \ which is a process that risks numeric problems.

[Bruno Luong]

"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]

ah, I must have been up all night again... I need more sleep :(

Answer 1

[John D'Errico]

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.

Answer 2

[KSSV]

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

[Walter Roberson]

roots(sym([1 -3 3 1]))

[Mustafa Yasar]

r = roots(p) ;

real(r) ??????????*

Are you for real?????

Click to these 1.0000's and see the result. THIS DOES NOT SOLVE MY PROBLEM!

Please do not write your answer if you do not know matlab at all.

[KSSV]

The real(r), gives all one. Attached below image.

Please, you try to understand the basics..first make your basics clear. I think you don't have strong good basics.

[Mustafa Yasar]

Impossible to explain something to you.

numerator_roots = roots([1 -3 3 -1]);

numerator_roots = real(numerator_roots);

[KSSV]

Yes agree...you cannot explain things clearly unless you have good idea of the problem.

[Walter Roberson]

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]

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.

Answer 3

[Matt J]

double(roots( sym([1 -3 3 -1] ))) == [1;1;1]  %all true

ans = 3x1 logical array
   1
   1
   1