View Javadoc
1   /*
2    * Licensed to the Apache Software Foundation (ASF) under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * The ASF licenses this file to You under the Apache License, Version 2.0
6    * (the "License"); you may not use this file except in compliance with
7    * the License.  You may obtain a copy of the License at
8    *
9    *      https://www.apache.org/licenses/LICENSE-2.0
10   *
11   * Unless required by applicable law or agreed to in writing, software
12   * distributed under the License is distributed on an "AS IS" BASIS,
13   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14   * See the License for the specific language governing permissions and
15   * limitations under the License.
16   */
17  
18  /*
19   * This is not the original file distributed by the Apache Software Foundation
20   * It has been modified by the Hipparchus project
21   */
22  
23  package org.hipparchus.analysis.solvers;
24  
25  import java.util.ArrayList;
26  import java.util.List;
27  
28  import org.hipparchus.analysis.QuinticFunction;
29  import org.hipparchus.analysis.UnivariateFunction;
30  import org.hipparchus.analysis.differentiation.DSFactory;
31  import org.hipparchus.analysis.differentiation.Derivative;
32  import org.hipparchus.analysis.differentiation.UnivariateDifferentiableFunction;
33  import org.hipparchus.exception.MathIllegalArgumentException;
34  import org.hipparchus.exception.MathIllegalStateException;
35  import org.hipparchus.util.FastMath;
36  import org.junit.Assert;
37  import org.junit.Test;
38  
39  /**
40   * Test case for {@link BracketingNthOrderBrentSolver bracketing n<sup>th</sup> order Brent} solver.
41   *
42   */
43  public final class BracketingNthOrderBrentSolverTest extends BaseSecantSolverAbstractTest {
44      /** {@inheritDoc} */
45      @Override
46      protected UnivariateSolver getSolver() {
47          return new BracketingNthOrderBrentSolver();
48      }
49  
50      /** {@inheritDoc} */
51      @Override
52      protected int[] getQuinticEvalCounts() {
53          return new int[] {1, 3, 8, 1, 9, 4, 8, 1, 12, 1, 16};
54      }
55  
56      @Test(expected=MathIllegalArgumentException.class)
57      public void testInsufficientOrder1() {
58          new BracketingNthOrderBrentSolver(1.0e-10, 1);
59      }
60  
61      @Test(expected=MathIllegalArgumentException.class)
62      public void testInsufficientOrder2() {
63          new BracketingNthOrderBrentSolver(1.0e-10, 1.0e-10, 1);
64      }
65  
66      @Test(expected=MathIllegalArgumentException.class)
67      public void testInsufficientOrder3() {
68          new BracketingNthOrderBrentSolver(1.0e-10, 1.0e-10, 1.0e-10, 1);
69      }
70  
71      @Test
72      public void testConstructorsOK() {
73          Assert.assertEquals(2, new BracketingNthOrderBrentSolver(1.0e-10, 2).getMaximalOrder());
74          Assert.assertEquals(2, new BracketingNthOrderBrentSolver(1.0e-10, 1.0e-10, 2).getMaximalOrder());
75          Assert.assertEquals(2, new BracketingNthOrderBrentSolver(1.0e-10, 1.0e-10, 1.0e-10, 2).getMaximalOrder());
76      }
77  
78      @Test
79      public void testConvergenceOnFunctionAccuracy() {
80          BracketingNthOrderBrentSolver solver =
81                  new BracketingNthOrderBrentSolver(1.0e-12, 1.0e-10, 0.001, 3);
82          QuinticFunction f = new QuinticFunction();
83          double result = solver.solve(20, f, 0.2, 0.9, 0.4, AllowedSolution.BELOW_SIDE);
84          Assert.assertEquals(0, f.value(result), solver.getFunctionValueAccuracy());
85          Assert.assertTrue(f.value(result) <= 0);
86          Assert.assertTrue(result - 0.5 > solver.getAbsoluteAccuracy());
87          result = solver.solve(20, f, -0.9, -0.2,  -0.4, AllowedSolution.ABOVE_SIDE);
88          Assert.assertEquals(0, f.value(result), solver.getFunctionValueAccuracy());
89          Assert.assertTrue(f.value(result) >= 0);
90          Assert.assertTrue(result + 0.5 < -solver.getAbsoluteAccuracy());
91      }
92  
93      @Test
94      public void testIssue716() {
95          BracketingNthOrderBrentSolver solver =
96                  new BracketingNthOrderBrentSolver(1.0e-12, 1.0e-10, 1.0e-22, 5);
97          UnivariateFunction sharpTurn = new UnivariateFunction() {
98              public double value(double x) {
99                  return (2 * x + 1) / (1.0e9 * (x + 1));
100             }
101         };
102         double result = solver.solve(100, sharpTurn, -0.9999999, 30, 15, AllowedSolution.RIGHT_SIDE);
103         Assert.assertEquals(0, sharpTurn.value(result), solver.getFunctionValueAccuracy());
104         Assert.assertTrue(sharpTurn.value(result) >= 0);
105         Assert.assertEquals(-0.5, result, 1.0e-10);
106     }
107 
108     @Test
109     public void testToleranceLessThanUlp() {
110         // function that is never zero
111         UnivariateFunction f = (x) -> x < 2.1 ? -1 : 1;
112         // tolerance less than 1 ulp(x)
113         UnivariateSolver solver = new BracketingNthOrderBrentSolver(0, 1e-18, 0, 5);
114 
115         // make sure it doesn't throw a maxIterations exception
116         double result = solver.solve(100, f, 0.0, 5.0);
117         Assert.assertEquals(2.1, result, 0.0);
118     }
119 
120     @Test
121     public void testFasterThanNewton() {
122         // the following test functions come from Beny Neta's paper:
123         // "Several New Methods for solving Equations"
124         // intern J. Computer Math Vol 23 pp 265-282
125         // available here: http://www.math.nps.navy.mil/~bneta/SeveralNewMethods.PDF
126         // the reference roots have been computed by the Dfp solver to more than
127         // 80 digits and checked with emacs (only the first 20 digits are reproduced here)
128         compare(new TestFunction(0.0, -2, 2) {
129             @Override
130             public <T extends Derivative<T>> T value(T x) {
131                 return x.sin().subtract(x.multiply(0.5));
132             }
133         });
134         compare(new TestFunction(6.3087771299726890947, -5, 10) {
135             @Override
136             public <T extends Derivative<T>> T value(T x) {
137                 return x.pow(5).add(x).subtract(10000);
138             }
139         });
140         compare(new TestFunction(9.6335955628326951924, 0.001, 10) {
141             @Override
142             public <T extends Derivative<T>> T value(T x) {
143                 return x.sqrt().subtract(x.reciprocal()).subtract(3);
144             }
145         });
146         compare(new TestFunction(2.8424389537844470678, -5, 5) {
147             @Override
148             public <T extends Derivative<T>> T value(T x) {
149                 return x.exp().add(x).subtract(20);
150             }
151         });
152         compare(new TestFunction(8.3094326942315717953, 0.001, 10) {
153             @Override
154             public <T extends Derivative<T>> T value(T x) {
155                 return x.log().add(x.sqrt()).subtract(5);
156             }
157         });
158         compare(new TestFunction(1.4655712318767680266, -0.5, 1.5) {
159             @Override
160             public <T extends Derivative<T>> T value(T x) {
161                 return x.subtract(1).multiply(x).multiply(x).subtract(1);
162             }
163         });
164 
165     }
166 
167     @Test
168     public void testSolverStopIteratingOnceSolutionIsFound() {
169         final double absoluteAccuracy = 0.1;
170         final double valueAccuracy = 1.;
171         BracketingNthOrderBrentSolver solver = new BracketingNthOrderBrentSolver(1e-14, absoluteAccuracy, valueAccuracy, 5);
172         FunctionHipparcus function = new FunctionHipparcus();
173         solver.solve(100, function, -100, 100);
174         Assert.assertEquals(1, function.values.stream().filter(value -> FastMath.abs(value) < valueAccuracy).count());
175         Assert.assertEquals(7, function.values.size());
176     }
177 
178     private static class FunctionHipparcus implements UnivariateFunction {
179 
180         List<Double> values = new ArrayList<>();
181 
182         @Override
183         public double value(double x) {
184             double value = -0.01 * x * x + 20;
185             values.add(value);
186             return value;
187         }
188     }
189 
190 
191     private void compare(TestFunction f) {
192         compare(f, f.getRoot(), f.getMin(), f.getMax());
193     }
194 
195     private void compare(final UnivariateDifferentiableFunction f,
196                          double root, double min, double max) {
197         NewtonRaphsonSolver newton = new NewtonRaphsonSolver(1.0e-12);
198         BracketingNthOrderBrentSolver bracketing =
199                 new BracketingNthOrderBrentSolver(1.0e-12, 1.0e-12, 1.0e-18, 5);
200         double resultN;
201         try {
202             resultN = newton.solve(100, f, min, max);
203         } catch (MathIllegalStateException tmee) {
204             resultN = Double.NaN;
205         }
206         double resultB;
207         try {
208             resultB = bracketing.solve(100, f, min, max);
209         } catch (MathIllegalStateException tmee) {
210             resultB = Double.NaN;
211         }
212         Assert.assertEquals(root, resultN, newton.getAbsoluteAccuracy());
213         Assert.assertEquals(root, resultB, bracketing.getAbsoluteAccuracy());
214 
215         // bracketing solver evaluates only function value, we set the weight to 1
216         final int weightedBracketingEvaluations = bracketing.getEvaluations();
217 
218         // Newton-Raphson solver evaluates both function value and derivative, we set the weight to 2
219         final int weightedNewtonEvaluations = 2 * newton.getEvaluations();
220 
221         Assert.assertTrue(weightedBracketingEvaluations < weightedNewtonEvaluations);
222 
223     }
224 
225     private static abstract class TestFunction implements UnivariateDifferentiableFunction {
226 
227         private final double root;
228         private final double min;
229         private final double max;
230 
231         protected TestFunction(final double root, final double min, final double max) {
232             this.root = root;
233             this.min  = min;
234             this.max  = max;
235         }
236 
237         public double getRoot() {
238             return root;
239         }
240 
241         public double getMin() {
242             return min;
243         }
244 
245         public double getMax() {
246             return max;
247         }
248 
249         public double value(final double x) {
250             return value(new DSFactory(0, 0).constant(x)).getValue();
251         }
252 
253         public abstract <T extends Derivative<T>> T value(final T t);
254 
255     }
256 
257 }