summaryrefslogtreecommitdiffstats
path: root/scilab/modules/polynomials/tests/unit_tests/poly_roots.dia.ref
blob: ea918869836f979718b1a87fa7a3ccc06c1ec6b5 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
// =============================================================================
// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
// Copyright (C) ????-2008 - INRIA - Michael Baudin
//
//  This file is distributed under the same license as the Scilab package.
// =============================================================================
// <-- JVM NOT MANDATORY -->
// poly_roots.tst --
//   Check the computation of the roots of a polynomial
//   with different kinds of polynomials and different
//   kinds of roots :
//   - real poly,
//   - complex poly,
//   - real roots,
//   - complex roots.
//roots : 3 real roots -> RPOLY
p=-6+11*%s-6*%s^2+%s^3;
myroots=roots(p);
computedroots = gsort(myroots);
expectedroots  = [3. ;2. ; 1.];
if (abs(computedroots-expectedroots)>400*%eps) then bugmes();quit;end
//roots : 3 real roots + polynomials algebra -> RPOLY
p=-6+11*%s-6*%s^2+%s^3;
myroots=roots(p+0);
computedroots = gsort(myroots);
expectedroots  = [3. ;2. ; 1.];
if (abs(computedroots-expectedroots)>400*%eps) then bugmes();quit;end
//roots : 3 complex roots -> Companion matrix used
p=-6-%i*6+(11+%i*5)*%s+(-6-%i)*%s^2+%s^3;
myroots=roots(p);
computedroots = gsort(myroots);
expectedroots  = [3. ;2. ; 1.+%i];
if (abs(computedroots-expectedroots)>400*%eps) then bugmes();quit;end
//roots : 3 complex roots + polynomials algebra -> Companion matrix used
p=-6-%i*6+(11+%i*5)*%s+(-6-%i)*%s^2+%s^3;
myroots=roots(p+0);
computedroots = gsort(myroots);
expectedroots  = [3. ;2. ; 1.+%i];
if (abs(computedroots-expectedroots)>400*%eps) then bugmes();quit;end
// roots : no root at all -> RPOLY
p=1;v=[];
if (roots(p)<>v) then bugmes();quit;end
if (roots(p+0)<>v) then bugmes();quit;end
//roots : 2 complex roots -> RPOLY
p=1+%s+%s^2;
myroots=roots(p);
computedroots = gsort(myroots);
expectedroots  = [-0.5 + sqrt(3.)/2.*%i ; -0.5 - sqrt(3.)/2.*%i];
if (abs(computedroots-expectedroots)>400*%eps) then bugmes();quit;end
//roots : 2 roots equals 0 -> RPOLY
p=%s^2;
myroots=roots(p);
computedroots = gsort(myroots);
expectedroots  = [0. ; 0. ];
if (abs(computedroots-expectedroots)>400*%eps) then bugmes();quit;end
//roots : 3 real roots -> forced to companion matrix with "e" option
p=-6+11*%s-6*%s^2+%s^3;
myroots=roots(p,"e");
computedroots = gsort(myroots);
expectedroots  = [3. ;2. ; 1.];
if (abs(computedroots-expectedroots)>400*%eps) then bugmes();quit;end
//roots : 2 complex roots -> forced to companion matrix with "e" option
p=1+%s+%s^2;
myroots=roots(p,"e");
computedroots = gsort(myroots);
expectedroots  = [-0.5 + sqrt(3.)/2.*%i ; -0.5 - sqrt(3.)/2.*%i];
if (abs(computedroots-expectedroots)>400*%eps) then bugmes();quit;end
// 2 real roots with a zero derivate at the root -> RPOLY
p=(%s-%pi)^2
 p  =
 
                              2  
    9.8696044 - 6.2831853s + s   
myroots=roots(p);
computedroots = gsort(myroots);
expectedroots  = [%pi;%pi];
if (abs(computedroots-expectedroots)>%eps) then bugmes();quit;end
// 2 real roots with a zero derivate at the root -> forced to companion matrix with "e" option
p=(%s-%pi)^2
 p  =
 
                              2  
    9.8696044 - 6.2831853s + s   
myroots=roots(p,"e");
computedroots = gsort(myroots);
expectedroots  = [%pi;%pi];
if (abs(computedroots-expectedroots)>10*%eps) then bugmes();quit;end
//
// Caution !
// The following are difficult root-finding problems
// with expected precision problems.
// See "Principles for testing polynomial
// zerofinding programs"
// Jenkins, Traub
// 1975
// p.28
// "The accuracy which one may expect to achieve in calculating
// zeros is limited by the condition of these zeros. In particular,
// for multiple zeros perturbations of size epsilon in the
// coefficients cause perturbations of size epsilon^(1/m)
// in the zeros."
//
//
// 3 real roots with a zero derivate at the root -> RPOLY
// *** PRECISION PROBLEM : only simple precision computed, instead of double precision ***
p=(%s-%pi)^3
 p  =
 
                                      2   3  
  - 31.006277 + 29.608813s - 9.424778s + s   
myroots=roots(p);
computedroots = gsort(myroots);
expectedroots  = [%pi;%pi;%pi];
if (abs(computedroots-expectedroots)>10^9*%eps) then bugmes();quit;end
// 3 real roots with a zero derivate at the root -> forced to companion matrix with "e" option
// *** PRECISION PROBLEM : not even simple precision computed, instead of double precision expected ***
p=(%s-%pi)^3
 p  =
 
                                      2   3  
  - 31.006277 + 29.608813s - 9.424778s + s   
myroots=roots(p,"e");
computedroots = gsort(myroots);
expectedroots  = [%pi;%pi;%pi];
if (abs(computedroots-expectedroots)>10^11*%eps) then bugmes();quit;end
// 4 real roots with a zero derivate at the root -> RPOLY
// *** PRECISION PROBLEM : only simple precision computed, instead of double precision expected ***
p=(%s-%pi)^4
 p  =
 
                                       2            3   4  
    97.409091 - 124.02511s + 59.217626s - 12.566371s + s   
myroots=roots(p);
computedroots = gsort(myroots);
expectedroots  = [%pi;%pi;%pi;%pi];
if (abs(computedroots-expectedroots)>10^9*%eps) then bugmes();quit;end
// 4 real roots with a zero derivate at the root -> forced to companion matrix with "e" option
// *** PRECISION PROBLEM : not even simple precision computed, instead of double precision expected ***
p=(%s-%pi)^4
 p  =
 
                                       2            3   4  
    97.409091 - 124.02511s + 59.217626s - 12.566371s + s   
myroots=roots(p,"e");
computedroots = gsort(myroots);
expectedroots  = [%pi;%pi;%pi;%pi];
if (abs(computedroots-expectedroots)>10^13*%eps) then bugmes();quit;end
// 10 real roots with a zero derivate at the root -> forced to companion matrix with "e" option
// *** PRECISION PROBLEM : only one one true digit ***
p=(%s-%pi)^10
 p  =
 
                                      2            3            4       
    93648.047 - 298090.99s + 426983.9s - 362435.19s + 201891.73s        
                        5            6            7           8         
            - 77116.961s + 20455.909s - 3720.7532s + 444.1322s          
                         9  10                                          
            - 31.415927s + s                                            
myroots=roots(p);
computedroots = gsort(myroots);
expectedroots  = [%pi;%pi;%pi;%pi;%pi;%pi;%pi;%pi;%pi;%pi];
if (abs(computedroots-expectedroots)>10^10*%eps) then bugmes();quit;end
// 10 real roots with a zero derivate at the root -> forced to companion matrix with "e" option
// *** PRECISION PROBLEM : only one one true digit ***
p=(%s-%pi)^10
 p  =
 
                                      2            3            4       
    93648.047 - 298090.99s + 426983.9s - 362435.19s + 201891.73s        
                        5            6            7           8         
            - 77116.961s + 20455.909s - 3720.7532s + 444.1322s          
                         9  10                                          
            - 31.415927s + s                                            
myroots=roots(p,"e");
computedroots = gsort(myroots);
expectedroots  = [%pi;%pi;%pi;%pi;%pi;%pi;%pi;%pi;%pi;%pi];
if (abs(computedroots-expectedroots)>10^15*%eps) then bugmes();quit;end