/** * Configures search parameters. * * @param maxCoefficients The maximum number of coefficients a polynomial will have. * @param searchRadius If {@code >} 0 then roots are searched for inside of -searchRadius ≤ x ≤ searchRadius. If * ≤ 0 then it will find all roots. * @param boundTolerance How accurately roots are found using bisection. Should be close enough that refinement can * accurately converge to the correct root. * @param maxBoundIterations Maximum number of iterations each step can perform when bounding. * @param maxRefineIterations Maximum number of iterations when refining. */ public FindRealRootsSturm(int maxCoefficients, double searchRadius , double boundTolerance, int maxBoundIterations , int maxRefineIterations ) { if( Double.isInfinite(searchRadius) ) searchRadius = -1; sturm = new SturmSequence(maxCoefficients); this.searchRadius = searchRadius; this.boundTolerance = boundTolerance; this.maxBoundIterations = maxBoundIterations; this.maxRefineIterations = maxRefineIterations; roots = new double[maxCoefficients]; }
/** * Configures search parameters. * * @param maxCoefficients The maximum number of coefficients a polynomial will have. * @param searchRadius If {@code >} 0 then roots are searched for inside of -searchRadius ≤ x ≤ searchRadius. If * ≤ 0 then it will find all roots. * @param boundTolerance How accurately roots are found using bisection. Should be close enough that refinement can * accurately converge to the correct root. * @param maxBoundIterations Maximum number of iterations each step can perform when bounding. * @param maxRefineIterations Maximum number of iterations when refining. */ public FindRealRootsSturm(int maxCoefficients, double searchRadius , double boundTolerance, int maxBoundIterations , int maxRefineIterations ) { if( Double.isInfinite(searchRadius) ) searchRadius = -1; sturm = new SturmSequence(maxCoefficients); this.searchRadius = searchRadius; this.boundTolerance = boundTolerance; this.maxBoundIterations = maxBoundIterations; this.maxRefineIterations = maxRefineIterations; roots = new double[maxCoefficients]; }
public static int countRealRoots( Polynomial poly ) { SturmSequence sturm = new SturmSequence(poly.size); sturm.initialize(poly); return sturm.countRealRoots(Double.NEGATIVE_INFINITY,Double.POSITIVE_INFINITY); }
public static int countRealRoots( Polynomial poly ) { SturmSequence sturm = new SturmSequence(poly.size); sturm.initialize(poly); return sturm.countRealRoots(Double.NEGATIVE_INFINITY,Double.POSITIVE_INFINITY); }
public int countRealRoots( Polynomial p , double low , double high ) { SturmSequence alg = new SturmSequence(p.size); alg.initialize(p); return alg.countRealRoots(low,high); }
private void compareSequences( Polynomial p, double value) { List<Double> expected = computeSturm(p,value); SturmSequence alg = new SturmSequence(p.size); alg.initialize(p); alg.computeFunctions(value); assertEquals(expected.size(),alg.sequenceLength); for( int j = 0; j < expected.size(); j++ ) { if( Double.isInfinite(expected.get(j)) ) { assertTrue(expected.get(j) == alg.f[j]); } else { assertEquals(expected.get(j),alg.f[j],Math.abs(alg.f[j])*1e-6); } } }
/** * Check sequence against a hand selected sequence */ @Test public void checkSequence() { Polynomial poly = Polynomial.wrap(-1,0,3,1,5,-3); SturmSequence alg = new SturmSequence(10); alg.initialize(poly); alg.computeFunctions(-3); List<Double> expected = computeSturm(poly,-3); assertEquals(expected.size(),alg.sequenceLength); for( int i = 0; i < expected.size(); i++ ) { assertEquals(expected.get(i),alg.f[i],1e-8); } }
@Test public void checkCount() { SturmSequence alg = new SturmSequence(10); alg.f = new double[]{0,1,1,-1,0,-1,-1,0,1,1}; alg.sequenceLength = alg.f.length; assertEquals(2,alg.countSignChanges()); alg.f = new double[]{1,1,1,-1,0,-1,-1,0,1,1,0,-1,-1,1}; alg.sequenceLength = alg.f.length; assertEquals(4,alg.countSignChanges()); }
/** * Several consistency checks on random polynomials to the number of roots within a random interval */ @Test public void rootCountConsistency_Random() { for( int i = 3; i < 20; i++ ) { Polynomial p = new Polynomial(i); for( int trial = 0; trial < 20; trial++ ) { for( int j = 0; j < p.size; j++ ) { p.c[j] = (rand.nextDouble()-0.5)*2; } SturmSequence alg = new SturmSequence(p.size()); double low = (rand.nextDouble() - 0.5)*200; double high = low + rand.nextDouble()*200; double middle = (low+high)/2.0; alg.initialize(p); int every = alg.countRealRoots(Double.NEGATIVE_INFINITY,Double.POSITIVE_INFINITY); int all = alg.countRealRoots(low,high); int lowN = alg.countRealRoots(low,middle); int highN = alg.countRealRoots(middle,high); assertTrue(all >= 0); assertTrue(lowN >= 0); assertTrue(highN >= 0); assertEquals(all, lowN + highN); assertTrue(all <= every); } } }
/** * Examine a case which was found to cause problems */ @Test public void checkSpecificPoly01() { Polynomial poly = Polynomial.wrap(-41.118263303597175,-120.95384505825373,-417.8477600492497,-634.5308297409192, -347.7885168491812,6.771313016808563,79.70258790927392,31.68212813610444,5.0248961592587875, 0.2879701466217739,0.0); // Compare computed sequence to the standard compareSequences( poly, -500); compareSequences( poly, Double.NEGATIVE_INFINITY); SturmSequence alg = new SturmSequence(poly.size); alg.initialize(poly); int N = alg.countRealRoots(Double.NEGATIVE_INFINITY,Double.POSITIVE_INFINITY); int M = alg.countRealRoots(-500,500); assertTrue(M <= N); }