[R] solving cubic/quartic equations non-iteratively -- comparisons

[Larry Hotchkiss]

Hi,

I'm responding to a post about finding roots of a cubic or quartic equation non-iteratively. One obviously could create functions using the explicit algebraic solutions. One post on the subject noted that the square-roots in those solutions also require iteration, and one post claimed iterative solutions are more accurate than the explicit solutions.

This post, however, is about comparative accuracy of (1) the R solve functionused in the included post, (2) the R polyroot function, (3) the Matlab roots function, (4) the SAS IML polyroot function. I tried the posted polynomial:

  -8 + 14*x - 7*x^2 + x^3 = 0

and a repeating-roots example:

  8 - 36*x + 54*x^2 - 27*x^3 = 0

I used Mathematica solutions as the reference:

  (* Posted example *)
  Roots[-8 + 14 x - 7 x^2 + x^3 == 0, x]
  x == 1 || x == 2 || x == 4

  (* Repeating-roots example  *)
  Roots[8 - 36 x + 54 x^2 - 27 x^3 == 0, x]
  
  x == 2/3 || x == 2/3 || x == 2/3

The results indicate that the R polyroot function is the most accurate for these examples. The R solve function is quite inaccurate for the repeating-roots example. It appears to be single-precision arithmetic on the real and imaginary parts. The same appears to be true for the Matlab function roots. SAS IML polyroot function appears to use double-precision calculations but was not nearly as accurate as the R polyroot function for the repeating-roots example.

The syntax and output for these examples are as follows:

 # ------------------------------------------------------------------------- #
 Mathematica:

  (* Posted example *)
  Roots[-8 + 14 x - 7 x^2 + x^3 == 0, x]
  x == 1 || x == 2 || x == 4

  (* Repeating-roots example  *)
  Roots[8 - 36 x + 54 x^2 - 27 x^3 == 0, x]
  
  x == 2/3 || x == 2/3 || x == 2/3

 # ------------------------------------------------------------------------- #
 R:

>   options(digits=16)
>   library(polynom) 

 # Posted example
>   p <- polynomial(c(-8,14,-7,1))
>   solve(p)
[1] 0.999999999999999 2.000000000000002 3.999999999999997
> 
>   polyroot(c(-8,14,-7,1))
[1] 1-0i 2+0i 4-0i
> 
 # Repeating-roots example
>   lp <- polynomial(c(8,-36, 54, -27))
>   solve(lp)
[1] 0.6666648437558125-0.000003157332198i 0.6666648437558125+0.000003157332198i
[3] 0.6666703124883749+0.000000000000000i
>   polyroot(lp)
[1] 6.666666666666670e-01+2e-16i 6.666666666666666e-01-8e-17i
[3] 6.666666666666664e-01-1e-16i
 
 # ------------------------------------------------------------------------- #
 Matlab:

>>  format long
>>  format compact

% Note: Matlab order of polynomial is reverse of R
% Posted example
>>  p = [-8,14,-7,1]; p = p(4:(-1):1)
p = 
    1    -7    14    -8
>>  roots(p)
ans =
   4.000000000000000
   2.000000000000004
   0.999999999999999

% Repeating-roots example
>>  lp = [8,-36, 54, -27]; lp = lp(4:(-1):1)
lp =
   -27    54   -36     8
>>  roots(lp)
ans =
  0.666669450738337 + 0.000004822149713i
  0.666669450738337 - 0.000004822149713i
  0.666661098523329                     
 # ------------------------------------------------------------------------- #
 SAS proc IML:
 proc IML;
  * Note: SAS polyroot also uses high-to-low order of the polynomial, reverse
          of R fuctions. *;
  * Posted example *;
   p = {-8 14 -7 1}; p = p[4:1];
   rts = polyroot(p);
   print p rts[format=19.16];

  * Repeating-roots example *;
   lp = {8 -36 54 -27}; lp = lp[4:1];
   rts = polyroot(lp);
   print lp rts[format=19.16];
  quit;
 
  * Output;
       p                 rts
        1  1.0000000000001500  0.0000000000000000
       -7  1.9999999999997700  0.0000000000000000
       14  4.0000000000000700  0.0000000000000000
       -8                                        


       lp                 rts
      -27  0.6666666666666700  0.0000000000000000
       54  0.6666666526177100  0.0000000000000000
      -36  0.6666666807156100  0.0000000000000000
        8 
 # ------------------------------------------------------------------------- #

Larry Hotchkiss





--------------------------------------------------------------------------------------
Message: 7
Date: Wed, 6 Jan 2010 13:03:14 +0100
From: "Kasper Kristensen" <kkr at aqua.dtu.dk>
Subject: Re: [R] solving cubic/quartic equations non-iteratively
To: <s02mjtj at math.ku.dk>
Cc: r-help at r-project.org
Message-ID:
<88431FCDD055A44BA5EEE1EC14E7E151143E05 at lu-mail-san.dfu.local>
Content-Type: text/plain; charset="iso-8859-1"


Try,

library(polynom)
p <- polynomial(c(-8,14,-7,1))
solve(p)


Kasper