Which is Better for Polynomial Equations: roots or solve?

Which one is better in terms of accuracy?

Answer 1

[John D'Errico]

Roots versus solve can be VERY different in the result. The answer is you need to understand the differences. When is increased precision meaningless?

For example, consider this simple problem, to find the three solutions of the polynomial:

(x-1)^3 = 0

Yes, I know, the solution is 1, a root with multiplicity of 3. But let me try it now:

P = [1 -3 3 -1];
format long g
roots(P)

ans = 
           1.00000657194364 +                     0i
          0.999996714028179 +  5.69145454681503e-06i
          0.999996714028179 -  5.69145454681503e-06i

Do you see the roots found are accurate only to about 5-6 significant digits? In fact, we expect a result that is no better than the cube root of eps for roots with multiplicty 3. But now try solve.

syms x
Psym = x^3 - 3*x^2 + 3*x - 1;
solve(Psym)

ans = 

The very same polynomial, but now essentially perfect results.

So in the case above, one would naturally prefer solve. As I said, they can be VERY different.

Next, consider a slightly different polynomial.

Psym = prod(x - (1:5))

Psym = 

Yes, I know you can find those roots yourself.

Psym = expand(Psym)

Psym = 

solve(Psym)

ans = 

And again, MATLAB finds the roots exactly. Yes, I'm cheating bit in some of these examples. They are really pretty easy.

Pcoef = double(flip(coeffs(Psym)))

Pcoef = 1×6

1 -15 85 -225 274 -120

roots(Pcoef)

ans = 5×1

5.00000000000014 3.99999999999978 3.00000000000012 1.99999999999997 1

And here roots find the solution, but because it is now a numerical method, working in double precision, there is some floating point trash down in the weeds. Roots was pretty good, but not perfect.

At the same time, I'd bet that roots will be WAY faster than solve.

timeit(@() roots(Pcoef))

ans = 
       2.0198380952381e-05

timeit(@() solve(Psym))

ans = 
               0.070501726

As expected, roots is several orders of magnitude faster than solve. This is a common tradeoff. In fact, on some problems, solve just never terminates, but numerical methods like roots are blazingly fast. Again, understanding what problem you are solving and the methods involved is crucial. Do you really need those extra digits? Do they mean anything to you?

Next, suppose you have a problem that is not so unusual in terms of what I see.

y = 4.363*x^3 + 25.42*x^2 - 148.4*x - 956.5

y = 

So each coefficient is given to 4 significant digits. Are those the EXACT VALUES OF THOSE NUMBERS? Are you kidding me? However those coefficients were derived, if we were able o compute them to higher precison, we would have more digits on the end for EACH coefficient. So now we might do this:

xsol = solve(y,'maxdegree',3)

xsol = 

vpa(xsol)

ans = 

And this is a significant problem with mathematical tools, in the hands of a novice. It is far too easy to look at that result, and imagine those are the EXACT solutions, down to 40 significant digits. Sorry, but that is a pile of complete and utter buffalo chips.

If you know those coefficients only to 4 significant digits, hoping to obtain 40 significant digits on a result is a meaningless waste of time.

I suppose it would be nice if MATLAB offered an interval arithmetic tool, that would allow me to compute those roots based on the knowledge that 4.363 might truly be written as the interval [4.3625,4.3635]. So the set of values for that coefficient that would have been rounded to 4.363.

Anyway, hoping to get high precision using solve on such a polynomial is silly in the extreme. Is roots better than solve there? WHO CARES? Just use roots, because again, it would be way faster, and those extra digits are meaningless.

(Do I need to expand on this, showing what the distribution of the set of roots, assuming that each of those coefficients was in fact uncertain, so now we would have a distribution for each root?)

P = vpa(expand(prod(x - (1:3))*.4567434325))

P = 

I have constructed P so that it has 3 roots, at [1 2 3]. Extract the coefficients, but I'll round them to 4 significant digits.

C = round(double(flip(coeffs(P))),4)

C = 1×4

0.4567 -2.7405 5.0242 -2.7405

N = 1e6;
noise = (rand(N,4)*2 - 1)*1e-4;
Csim = C + noise;
R = zeros(N,3);
for i = 1:N
  R(i,:) = sort(roots(Csim(i,:)));
end
histogram(R(:,1),100)

histogram(R(:,2),100)

histogram(R(:,3),100)

So the 3 roots in fact are not truly known. Hoping for an exact root from rounded coefficients is a waste of time.

In the end, what matters is using common sense, having an understanding of the various methods, and the problems they will sometimes suffer. What matters is understanding the problem you want to solve. If you NEED a symbolic result for the next step, then roots may be the wrong tool.

Answer 2

[Kareem Elgindy]

An outstanding answer!