summaryrefslogtreecommitdiffstats
path: root/scilab/modules/polynomials/tests/unit_tests/poly_roots.tst
blob: da35f6682eecf42e288877bb8d458c6a1c3f8989 (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
// =============================================================================
// 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 pause,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 pause,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 pause,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 pause,end
// roots : no root at all -> RPOLY
p=1;v=[];
if (roots(p)<>v) then pause,end
if (roots(p+0)<>v) then pause,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 pause,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 pause,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 pause,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 pause,end
// 2 real roots with a zero derivate at the root -> RPOLY
p=(%s-%pi)^2
myroots=roots(p);
computedroots = gsort(myroots);
expectedroots  = [%pi;%pi];
if (abs(computedroots-expectedroots)>%eps) then pause,end
// 2 real roots with a zero derivate at the root -> forced to companion matrix with "e" option
p=(%s-%pi)^2
myroots=roots(p,"e");
computedroots = gsort(myroots);
expectedroots  = [%pi;%pi];
if (abs(computedroots-expectedroots)>10*%eps) then pause,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
myroots=roots(p);
computedroots = gsort(myroots);
expectedroots  = [%pi;%pi;%pi];
if (abs(computedroots-expectedroots)>10^9*%eps) then pause,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
myroots=roots(p,"e");
computedroots = gsort(myroots);
expectedroots  = [%pi;%pi;%pi];
if (abs(computedroots-expectedroots)>10^11*%eps) then pause,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
myroots=roots(p);
computedroots = gsort(myroots);
expectedroots  = [%pi;%pi;%pi;%pi];
if (abs(computedroots-expectedroots)>10^9*%eps) then pause,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
myroots=roots(p,"e");
computedroots = gsort(myroots);
expectedroots  = [%pi;%pi;%pi;%pi];
if (abs(computedroots-expectedroots)>10^13*%eps) then pause,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
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 pause,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
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 pause,end