[Scilab-Dev] Problems with function roots

[Michaël Baudin]

Dear Pedro,

The following script is a simple way of checking
that, as Serge wrote, the companion matrix gives the correct result.
It explicitly uses the companion matrix to compute the
roots.
There is a small step to reverse the coefficients of the
polynomial ; indeed, "roots" expects the coefficients in
decreasing degree order, while "poly" expects the coefficients
in increasing degree order.

v= [1.12119799 0 3.512D+13 32 3.275D+27 0 1.117D+41 4.952D+27 1.722D+54 
0 1.224D+67 0 3.262D+79 ];
r1 = roots(v) // may fail
r2 = roots(v,"e")
dv = size(v,"*")
p = poly(v(dv:-1:1),"x","coeff")
A = companion(p)
r3 = spec(A)
plot(real(r3),imag(r3),"bo")

The following script measures the numerical difficulty of this
eigenvalue problem.

[R,D] = spec(A);
cond(R)

The condition number of the matrix of right
eigenvectors is extremely large (~10^72, which is much larger than 10^16).
This problem is numerically difficult so we are
warned that the eigenvalues that we get are extremely
sensitive on the entries of the companion matrix.

As an example of this, consider an example by Wilkinson,
presented by Moler in NCM, section 10.4 "Characteristic Polynomial".
The exact roots of the following polynomial are 1, 2, ..., 20, but
the computed roots are varying slightly from this.
The condition number of the matrix of right eigenvectors is ~10^20.

A = diag(1:20)
p = prod((1:20)' - x)
roots(p)
[R,D] = spec(companion(p));
cond(R)

Best regards,

Michaël

PS
We can get a more inaccurate result by computing the
determinant of the matrix A, as in the following script.
This time, the eigenvalues are completely wrong, (one of the
computed eigenvalues is ~ -23.07, i.e. even the sign is wrong !)
because the coefficients of the polynomial produced by the det function
were not computed sufficiently accurately, because of the
limited precision of doubles and the numerical difficulty of
this particular problem.

x = poly(0,"x")
p = det(A - x * eye())
roots(p)

This is an example where the exact mathematical
statement ("the eigenvalues are the roots of the matrix A - lambda I")
leads to a completely wrong numerical computation.

Serge Steer a écrit :
> Le 17/07/2010 23:28, Pedro Ledoux a écrit :
>>
>> Hello
>>
>> I'm Pedro Ledoux, I'm a electrical engineering student. At university 
>> I try to do as many as possible using free softwares including 
>> Scilab. During the project of a passive high pass filter  I had to 
>> get the roots of a 12 order polynomial function. So I typed in Scilab 
>> environment:
>>
>>  v= [1.12119799 0 3.512D+13 32 3.275D+27 0 1.117D+41 4.952D+27 
>> 1.722D+54 0 1.224D+67 0 3.262D+79 ];
>> roots(v)
>>
>> Unfortunately Scilab failed. The message :" !--error 24 Convergence 
>> problems,," has appear. I didn't tryed it in Matlab but GNU octave 
>> was able to calculate those roots. By now I'm not a great programmer 
>> yet and the only way that I have to contribute to Scilab is using in 
>> and reporting bugs.
>>
>> If almost developers are French I can speak French also if it is more 
>> comfortable.
>>
>>
>>   
> You can use roots(v,"e"). In this the algorithm computes the 
> eigenvalues of the companion matrix. It is a slower algorithm, but it 
> never fail.
>
> Serge Steer
>


-- 
Michaël Baudin
Ingénieur de développement
michael.baudin at scilab.org
-------------------------
Consortium Scilab - Digiteo
Domaine de Voluceau - Rocquencourt
B.P. 105 - 78153 Le Chesnay Cedex
Tel. : 01 39 63 56 87 - Fax : 01 39 63 55 94