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