Eigenvalues to find roots of polynomials

Eigenvalues to find roots of polynomials

[marius63]

I am trying to use Eigenvalues methods to find the roots of
polynomials... and failing miserably!
My math ability is not stellar. I used the information at

(1)
http://mathworld.wolfram.com/PolynomialRoots.html (equation 12)

and

(2)
http://en.wikipedia.org/wiki/Root_finding
(One possibility is to form the companion matrix of the polynomial.
Since the eigenvalues of this matrix coincide with the roots of the
polynomial, one can use any eigenvalue algorithm to find the roots of
the polynomial.)
http://en.wikipedia.org/wiki/Companion_matrix

I tried both schemes with the test data in (1). The output is first
the matrix, for verification, then the (wrong) roots (they should be
-1, 1, 2) and the result of applying them to the polynomial (and
yielding miserably non-zero values of course). The source code
follows. What am I missing? Thanks a lot.

$ ./find_roots 1 -2 -1 2
 5.00000E-01 1.00000E+00-5.00000E-01
 1.00000E+00 0.00000E+00 0.00000E+00
 0.00000E+00 1.00000E+00 0.00000E+00
-7.6156E-01=> 4.2003E-01
 4.1249E+00=>-1.6015E+01
 6.3667E-01=> 5.9465E-01

$ ./find_roots_wikipedia -2 -1 2
 0.00000000000000E+00,  0.00000000000000E+00, -2.00000000000000E+00,
 1.00000000000000E+00,  0.00000000000000E+00,  1.00000000000000E+00,
 0.00000000000000E+00,  1.00000000000000E+00,  2.00000000000000E+00,
-1.17008648662603E+00 => -3.38499151386795E-02
 6.88892182534018E-01 =>  1.96496630422799E+00
 2.48119430409202E+00 =>  3.00688836109107E+01


--  Trying to compute the roots of a polynomial using the method
--  described in eq. (12) of http://mathworld.wolfram.com/PolynomialRoots.html
--  (consulted July 2011)

with Ada.Command_Line;
with Ada.Numerics.Generic_Real_Arrays;
with Ada.Text_IO;
with Ada.Numerics.Generic_Complex_Types;
with Ada.Numerics.Generic_Complex_Arrays;

procedure Find_Roots is
   type Real is digits 5;
   package Real_Arrays is new Ada.Numerics.Generic_Real_Arrays (Real);
   package Complex_Types is new Ada.Numerics.Generic_Complex_Types
(Real);
   package Complex_Arrays is new Ada.Numerics.Generic_Complex_Arrays
(Real_Arrays, Complex_Types);
   N, M : Natural;
begin
   N := Ada.Command_Line.Argument_Count;
   --  N = number of arguments = degree of the polynomial.
   --  The input consists of all coficients in descending degree,
   --  e.g. for polynomial
   --     x^3 - 2x^2 - x + 2
   --  the input is
   --     1 -2 -1 2
   --  Roots of the above (for testing):
   --     -1, 1, 2

   declare
      C : Complex_Arrays.Complex_Matrix (1 .. N - 1, 1 .. N - 1) :=
(others => (others => (0.0, 0.0)));
      --  Using a complex matrix because real matrices for eigenvalues
must be symmetric.
      --  Below we convert all real values to complex (with null Im
part).
      --  The complex matrix must be Hamiltonian;
      --  I thought a non-symmetric converted to complex would
      --  yield a non-Hamiltonian matrix, but I tried and it
worked :-)... hmm... NOT :-(
      A0 : Real := Real'Value (Ada.Command_Line.Argument (N));
      use Ada.Text_IO;
   begin
      for I in 1 .. N - 1 loop
         C (1, I).Re := - Real'Value (Ada.Command_Line.Argument (N -
I)) / A0;
         if I < N - 1 then
            C (I + 1, I).Re := 1.0;
         end if;
      end loop;

      --  inspect matrix
      for i in c'range(1) loop
         for j in c'range(2) loop
            put (c(i,j).re'img);
         end loop;
         new_line;
      end loop;

      declare
         V : Real_Arrays.Real_Vector := Complex_Arrays.Eigenvalues
(C);
         use Real_Arrays;
         Z : Real;
         use Ada.Text_IO;
      begin
         for I in V'Range loop
            V (I) := 1.0 / V (I);
            Put (Real'Image (V (I)));
            Z := 0.0;
            for J in 1 .. N - 1 loop
               Z := Z + C (1, J).Re * V (I) ** J;
            end loop;
            Put ("=>" & Real'Image (Z));
            New_Line;
         end loop;
      end;
   end;
end;

--  Build with:
--  gnatmake find_roots -largs /Developer/SDKs/MacOSX10.4u.sdk/usr/lib/
libblas.dylib /Developer/SDKs/MacOSX10.4u.sdk/usr/lib/liblapack.dylib


--  Given an Nth degree polynomial A(N)X**N + ... + A(1)X + A(0),
--  the roots can be found by finding the eigenvalues L(I) of the
matrix
--     _                                                           _
--    |   -A(1)/A(0)   -A(2)/A(0)   -A(3)/A(0)   ...   -A(N)/A(0)   |
--    |        1            0            0       ...        0       |
--    |        0            1            0       ...        0       |
--    |       ...          ...           1       ...        0       |
--    |_       0            0            0       ...        0      _|
--
--  and taking R(i) = 1/L(I). This method can be computationally
expensive,
--  but is fairly robust at finding close and multiple roots.
--    (wolphram)


--  This program computes the roots of a polynomial using the method
--  described in http://en.wikipedia.org/wiki/Companion_matrix

--  This method is for monic polynomials, i.e. with coficient of the
term
--  of greater degree equal to 1, so the input is for the remainder
terms,
--  in descending order e.g. for polynomial
--     x^3 - 2x^2 - x + 2
--  the input is
--     -2 -1 2

with Ada.Command_Line;
with Ada.Numerics.Generic_Real_Arrays;
with Ada.Text_IO;
with Ada.Numerics.Generic_Complex_Types;
with Ada.Numerics.Generic_Complex_Arrays;

procedure Find_Roots_Wikipedia is
   type Real is digits 15;
   package Real_Arrays is new Ada.Numerics.Generic_Real_Arrays (Real);
   package Complex_Types is new Ada.Numerics.Generic_Complex_Types
(Real);
   package Complex_Arrays is new Ada.Numerics.Generic_Complex_Arrays
(Real_Arrays, Complex_Types);
   N : Natural;
begin
   N := Ada.Command_Line.Argument_Count;
   declare
      C : Complex_Arrays.Complex_Matrix (1 .. N, 1 .. N) := (others =>
(others => (0.0, 0.0)));
      --  Using a complex matrix because real matrices for eigenvalues
must be symmetric.
      --  Below we convert all real values to complex (with null Im
part).
      --  (Still, the complex matrix must be Hamiltonian...)
      use Ada.Text_IO;
   begin
      for I in 1 .. N loop
         C (I, N).Re := - Real'Value (Ada.Command_Line.Argument (N - I
+ 1));
         if I > 1 then
            C (I, I - 1).Re := 1.0;
         end if;
      end loop;

      --  inspect matrix
      for i in c'range(1) loop
         for j in c'range(2) loop
            put (c(i,j).re'img & ", ");
         end loop;
         new_line;
      end loop;

      declare
         V : Real_Arrays.Real_Vector := Complex_Arrays.Eigenvalues
(C);
         use Real_Arrays;
         Z : Real;
         use Ada.Text_IO;
      begin
         for I in V'Range loop
            Put (Real'Image (V (I)));
            Z := V (I) ** N;
            for J in 1 .. N loop
               Z := Z + C (N, J).Re * V (I) ** (J - 1);
            end loop;
            Put (" => " & Real'Image (Z));
            New_Line;
         end loop;
      end;
   end;
end;

--  Build with:
--  gnatmake find_roots_wikipedia -largs /Developer/SDKs/
MacOSX10.4u.sdk/usr/lib/libblas.dylib /Developer/SDKs/MacOSX10.4u.sdk/
usr/lib/liblapack.dylib

Re: Eigenvalues to find roots of polynomials

[Dmitry A. Kazakov]

On Mon, 25 Jul 2011 10:36:14 -0700 (PDT), marius63 wrote:

> --    |   -A(1)/A(0)   -A(2)/A(0)   -A(3)/A(0)   ...   -A(N)/A(0)   |
> --    |        1            0            0       ...        0       |
> --    |        0            1            0       ...        0       |
> --    |       ...          ...           1       ...        0       |
> --    |_       0            0            0       ...        0      _|

The matrix passed to Eigenvalues (ARM G.3.2) must be Hermitian. The above
matrix is not Hermitian.

P.S. A matrix M is Hermitian if M (I,J) if the complex conjugate of M
(J,I)), i.e. Re (M(I,J)) = Re (M(J,I)) and Im (M(I,J)) = -Im (M(J,I)) 

-- 
Regards,
Dmitry A. Kazakov
http://www.dmitry-kazakov.de

Re: Eigenvalues to find roots of polynomials

[Simon Wright]

marius63 <amado.alves@gmail.com> writes:
[...]

See http://sourceforge.net/projects/gnat-math-extn/ for work I did
supporting asymmetric real and non-hermitian complex matrices.

> --  Build with:
> --  gnatmake find_roots -largs /Developer/SDKs/MacOSX10.4u.sdk/usr/lib/
> libblas.dylib /Developer/SDKs/MacOSX10.4u.sdk/usr/lib/liblapack.dylib

You shouldn't need to specify the libraries - which GNAT are you using?
On this Mac with GNAT GPL 2011 (built for Snow Leopard),

  $ gnatmake find_roots
  gcc -c find_roots.adb
  find_roots.adb:17:07: warning: variable "M" is never read and never assigned
  gnatbind -x find_roots.ali
  gnatlink find_roots.ali
  $ 

Even if your Ada can't manage it itself, I'd have thought
  -largs -llapack -lblas
would do.

Re: Eigenvalues to find roots of polynomials

[John B. Matthews]

In article 
<3ec148c9-a033-40a1-bc32-1aa52867478a@e35g2000yqc.googlegroups.com>,
 marius63 <amado.alves@gmail.com> wrote:

> (2) <http://en.wikipedia.org/wiki/Root_finding>

Later in the same article, this alternative:

<http://en.wikipedia.org/wiki/Durand–Kerner_method>

is implemented here:

<http://home.roadrunner.com/~jbmatthews/misc/groots.html>

-- 
John B. Matthews
trashgod at gmail dot com
<http://sites.google.com/site/drjohnbmatthews>

Re: Eigenvalues to find roots of polynomials

[Ada novice]

On Jul 25, 9:22 pm, Simon Wright <si...@pushface.org> wrote:
> marius63 <amado.al...@gmail.com> writes:
>
> [...]
>
> Seehttp://sourceforge.net/projects/gnat-math-extn/for work I did
> supporting asymmetric real and non-hermitian complex matrices.

Thanks Simon for making this possible. I have been away from Ada for a
few months. Last week I looked at your codes and tested them for a few
cases. The codes worked fine :) as expected. In structural dynamics,
one can have three matrices to deal with. These three matrices are
combined with each other to form then the matrix A and B that you have
in your codes and hence it is a generalized eigenproblem. But the
eigenvalues are obtained in conjugate pairs as we are dealing with
double-sized matrices. I am waiting for a book to arrive here and then
I will formulate and use your Ada codes. For such a problem, I have my
codes in Matlab to compare with.

A very useful theorem that I learned several years back is the
Gershgorin Theorem. It lets you know the bounds for your eigenvalues.
More information here:

http://en.wikipedia.org/wiki/Gershgorin_circle_theorem


Cheers
YC

Re: Eigenvalues to find roots of polynomials

[Simon Wright]

Ada novice <yogeshwarsing@gmx.com> writes:

> I have been away from Ada for a few months. Last week I looked at your
> codes and tested them for a few cases.

I've been trying to reply to your e-mails, there's something wrong with
the communications channel!

Re: Eigenvalues to find roots of polynomials

[marius63]

> See http://sourceforge.net/projects/gnat-math-extn/ for work I did
> supporting asymmetric real and non-hermitian complex matrices. (Simon)

Thanks. Cannot reset my sourceforge password though. How you guys did
it?

> > --  gnatmake find_roots -largs /Developer/SDKs/MacOSX10.4u.sdk/usr/lib/
> > libblas.dylib /Developer/SDKs/MacOSX10.4u.sdk/usr/lib/liblapack.dylib
>
> You shouldn't need to specify the libraries - which GNAT are you using?

I know - old Mac, old GNAT...

Re: Eigenvalues to find roots of polynomials

[Simon Wright]

marius63 <amado.alves@gmail.com> writes:

> Thanks. Cannot reset my sourceforge password though. How you guys did
> it?

From the e-mail that SF sent:

"To access the site again, you'll need to go through the email
recovery process and choose a shiny new password:

https://sourceforge.net/account/registration/recover.php

If you need help with this, feel free to e-mail us:

sfnet_ops@geek.net"

Re: Eigenvalues to find roots of polynomials

[marius63]

> > --    |   -A(1)/A(0)   -A(2)/A(0)   -A(3)/A(0)   ...   -A(N)/A(0)   |
> > --    |        1            0            0       ...        0       |
> > --    |        0            1            0       ...        0       |
> > --    |       ...          ...           1       ...        0       |
> > --    |_       0            0            0       ...        0      _|
>
> The matrix passed to Eigenvalues (ARM G.3.2) must be Hermitian. The above
> matrix is not Hermitian.
>
> P.S. A matrix M is Hermitian if M (I,J) if the complex conjugate of M
> (J,I)), i.e. Re (M(I,J)) = Re (M(J,I)) and Im (M(I,J)) = -Im (M(J,I))
> (Dmitry A. Kazakov)

I suspected that. I couldn't tell an hermitian matrix from a
hamiltonian gargle blaster so when the program did not raise
Argument_Error I assumed the matrix was kosher enough. It seems not.
(So, Ada Eigenvalues is not a very good way to find roots, is it?)
Thanks.

Re: Eigenvalues to find roots of polynomials

[marius63]

On Jul 26, 3:04 pm, Simon Wright <si...@pushface.org> wrote:
> marius63 <amado.al...@gmail.com> writes:
> > Thanks. Cannot reset my sourceforge password though. How you guys did
> > it?
>
> From the e-mail that SF sent:
>
> "To access the site again, you'll need to go through the email
> recovery process and choose a shiny new password:
>
> https://sourceforge.net/account/registration/recover.php

I know this page. Looks like something out of the sirius cybernetic
corporation. It offers one(!) option: recover by email. And that's
all. No email address field, no submit button. That's why I asked how
you did it. Voice activation?

Re: Eigenvalues to find roots of polynomials

[marius63]

And the winner is...

> <http://home.roadrunner.com/~jbmatthews/misc/groots.html>

John B. Matthews Ada...Generic_Roots works like a charm.
It passed the simple test case, and solved my quintic problem x^5 -
4178x + 4177

$ ./find_roots_mathews 1 0 0 0 -4178 4177
Index =>  0, Position =>  6, Value => 4177
Index =>  1, Position =>  5, Value => -4178
Index =>  2, Position =>  4, Value => 0
Index =>  3, Position =>  3, Value => 0
Index =>  4, Position =>  2, Value => 0
Index =>  5, Position =>  1, Value => 1
Root 1 => (Re=> 1.00000000000000E+00, Im=> 0.00000000000000E+00)
Root 2 => (Re=>-2.47582467311450E-01, Im=>-8.05881611194044E+00)
Root 3 => (Re=> 7.76752669636581E+00, Im=>-1.40129846432482E-45)
Root 4 => (Re=>-2.47582467311450E-01, Im=> 8.05881611194044E+00)
Root 5 => (Re=>-8.27236176174291E+00, Im=> 0.00000000000000E+00)

I know from previous mathematical derivation (the polynomial comes
from extracting a variable in a summation), that the solution is a non-
negative scalar different from one, so Root 3 is my number, with the
non-zero but very small value of the imaginary part certainly only a
residue of computation done inside Matthews's magic machine. Thanks a
very great lot. Source code follows. Sorry for the mispelling,
"Mathews".

--  Trying to compute the roots of a polynomial using John Mathews
procedure.

with Ada.Command_Line;
with Ada.Numerics.Generic_Real_Arrays;
with Ada.Text_IO;
with Ada.Numerics.Generic_Complex_Types;
with Ada.Numerics.Generic_Complex_Arrays;
with Ada.Numerics.Generic_Complex_Arrays.Generic_Roots;

procedure Find_Roots_Mathews is
   type Real is digits 15;
   package Real_Arrays is new Ada.Numerics.Generic_Real_Arrays (Real);
   package Complex_Types is new Ada.Numerics.Generic_Complex_Types
(Real);
   package Complex_Arrays is new Ada.Numerics.Generic_Complex_Arrays
(Real_Arrays, Complex_Types);
   procedure Complex_Roots is new Complex_Arrays.Generic_Roots;
   M, N : Natural;
begin
   M := Ada.Command_Line.Argument_Count;
   N := M - 1;
   --  N = number of arguments = degree of the polynomial + 1.
   --  The input consists of all coficients in descending degree,
   --  e.g. for polynomial
   --     x^3 - 2x^2 - x + 2
   --  the input is
   --     1 -2 -1 2
   --  Roots of the above (for testing):
   --     -1, 1, 2

   declare
      P : Complex_Arrays.Complex_Vector (0 .. N);
      R : Complex_Arrays.Complex_Vector (1 .. N);
      use Ada.Text_IO;
   begin
      for I in 0 .. N loop
         Put_Line ("Index => " & I'img &
                   ", Position => " & Integer'Image(M - I) &
                   ", Value => " & Ada.Command_Line.Argument (M - I));
         P (I) := (Re => Real'Value (Ada.Command_Line.Argument (M -
I)), Im => 0.0);
      end loop;
      Complex_Roots (P, R);
      for I in R'Range loop
         Put_Line ("Root" & I'img &
                   " => (Re=>" & R(I).Re'img &
                   ", Im=>" & R(I).Im'img & ")");
      end loop;
   end;
end;

Re: Eigenvalues to find roots of polynomials

[Simon Wright]

marius63 <amado.alves@gmail.com> writes:

> On Jul 26, 3:04 pm, Simon Wright <si...@pushface.org> wrote:
>> marius63 <amado.al...@gmail.com> writes:
>> > Thanks. Cannot reset my sourceforge password though. How you guys did
>> > it?
>>
>> From the e-mail that SF sent:
>>
>> "To access the site again, you'll need to go through the email
>> recovery process and choose a shiny new password:
>>
>> https://sourceforge.net/account/registration/recover.php
>
> I know this page. Looks like something out of the sirius cybernetic
> corporation. It offers one(!) option: recover by email. And that's
> all. No email address field, no submit button. That's why I asked how
> you did it. Voice activation?

The 'Email Recovery' radio button needs to be selected, then a text box
for you to enter your email magically appears (what HCI design!  a
winner!)

Re: Eigenvalues to find roots of polynomials

[Georg Bauhaus]

On 26.07.11 16:47, marius63 wrote:
> On Jul 26, 3:04 pm, Simon Wright <si...@pushface.org> wrote:

>> https://sourceforge.net/account/registration/recover.php
> 
> I know this page. Looks like something out of the sirius cybernetic
> corporation. It offers one(!) option: recover by email. And that's
> all. No email address field, no submit button. That's why I asked how
> you did it. Voice activation?

No, Netscape Lisp, a.k.a. Javascript is needed. The radio button is
UI polymorphic, and has a handler that will take you to the next
step.

Re: Eigenvalues to find roots of polynomials

[marius63]

On Jul 26, 7:04 pm, Georg Bauhaus <rm.dash-bauh...@futureapps.de>
wrote:
> On 26.07.11 16:47, marius63 wrote:
>
> > On Jul 26, 3:04 pm, Simon Wright <si...@pushface.org> wrote:
> >>https://sourceforge.net/account/registration/recover.php
>
> > I know this page. Looks like something out of the sirius cybernetic
> > corporation. It offers one(!) option: recover by email. And that's
> > all. No email address field, no submit button. That's why I asked how
> > you did it. Voice activation?
>
> No, Netscape Lisp, a.k.a. Javascript is needed. The radio button is
> UI polymorphic, and has a handler that will take you to the next
> step.

Yeah, well, it's buggy. The relevant code fragment is

<input type="radio" id="method" name="method" value="email"
onClick="select_method('#form_email');" />

and the Javascript debugger says "reference error: cant find variable
select_method"

Maybe this system, Safari Version 4.1.3 (4533.19.4) on Mac OS X 10.4,
is too old.

/*
Chrome is not available for 10.4
Someday will upgrade Mac OS.
Not now.
Too many critical stuff installed.
*/

Will try from a different computer soon.
Thanks.

Re: Eigenvalues to find roots of polynomials

[John B. Matthews]

In article 
<8d424bae-53cb-4418-8ec1-1fbf3ec260e3@y16g2000yqk.googlegroups.com>,
 marius63 <amado.alves@gmail.com> wrote:

> And the winner is...
> 
> > <http://home.roadrunner.com/~jbmatthews/misc/groots.html>
> 
> John B. Matthews Ada...Generic_Roots works like a charm.
> It passed the simple test case, and solved my quintic problem
> x^5 - 4178x + 4177
> 
> $ ./find_roots_mathews 1 0 0 0 -4178 4177
> Index =>  0, Position =>  6, Value => 4177
> Index =>  1, Position =>  5, Value => -4178
> Index =>  2, Position =>  4, Value => 0
> Index =>  3, Position =>  3, Value => 0
> Index =>  4, Position =>  2, Value => 0
> Index =>  5, Position =>  1, Value => 1
> Root 1 => (Re=> 1.00000000000000E+00, Im=> 0.00000000000000E+00)
> Root 2 => (Re=>-2.47582467311450E-01, Im=>-8.05881611194044E+00)
> Root 3 => (Re=> 7.76752669636581E+00, Im=>-1.40129846432482E-45)
> Root 4 => (Re=>-2.47582467311450E-01, Im=> 8.05881611194044E+00)
> Root 5 => (Re=>-8.27236176174291E+00, Im=> 0.00000000000000E+00)
> 
> I know from previous mathematical derivation (the polynomial comes 
> from extracting a variable in a summation), that the solution is a 
> non- negative scalar different from one, so Root 3 is my number, with 
> the non-zero but very small value of the imaginary part certainly 
> only a residue of computation done inside Matthews's magic machine. 
> Thanks a very great lot. Source code follows. Sorry for the 
> mispelling, "Mathews".

No problem; I'm glad you found it useful! You might also like the test 
program, croot.adb. For example, adding this line:

Show_Result((4177.0, -4178.0, 0.0, 0.0, 0.0, 1.0));

produces the following table:

Poly:  1.00x^5 + 0.00x^4 + 0.00x^3 + 0.00x^2 - 4178.00x^1 + 4177.00
Real: -8.27236176174291E+00
Comp: -2.47582467311450E-01 +- 8.05881611194044E+00i
Real:  1.00000000000000E+00
Real:  7.76752669636581E+00
Largest error: = 9.09494701772928E-12

The program sorts the results to find conjugate pairs and checks values 
against zero using the same epsilon found in the implementation. It then 
evaluates the polynomial at each root, recording the largest deviation 
from zero.

[Example omitted.]
-- 
John B. Matthews
trashgod at gmail dot com
<http://sites.google.com/site/drjohnbmatthews>

Re: Eigenvalues to find roots of polynomials

[Ada novice]

On Jul 26, 11:44 am, Simon Wright <si...@pushface.org> wrote:

> I've been trying to reply to your e-mails, there's something wrong with
> the communications channel!

Hi Simon, I have just sent an email to you and I hope we can continue
from there.

Cheers
YC

Re: Eigenvalues to find roots of polynomials

[marius63]

> > > <http://home.roadrunner.com/~jbmatthews/misc/groots.html>
...
> You might also like the test program, croot.adb.

I do, I've read this program, it's wonderful, and I'll use it someday,
but now I needed to avoid Makefiles, and found writing my own test
program easier for that. Thanks.