summaryrefslogtreecommitdiffstats log msg author committer range
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 ``````