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  package org.hipparchus.analysis.solvers;
23  
24  import org.hipparchus.analysis.FunctionUtils;
25  import org.hipparchus.analysis.MonitoredFunction;
26  import org.hipparchus.analysis.QuinticFunction;
27  import org.hipparchus.analysis.UnivariateFunction;
28  import org.hipparchus.analysis.differentiation.DSFactory;
29  import org.hipparchus.analysis.differentiation.UnivariateDifferentiableFunction;
30  import org.hipparchus.analysis.function.Constant;
31  import org.hipparchus.analysis.function.Inverse;
32  import org.hipparchus.analysis.function.Sin;
33  import org.hipparchus.analysis.function.Sqrt;
34  import org.hipparchus.exception.MathIllegalArgumentException;
35  import org.hipparchus.exception.MathIllegalStateException;
36  import org.hipparchus.util.FastMath;
37  import org.junit.Assert;
38  import org.junit.Test;
39  
40  /**
41   * Test case for {@link BrentSolver Brent} solver.
42   * Because Brent-Dekker is guaranteed to converge in less than the default
43   * maximum iteration count due to bisection fallback, it is quite hard to
44   * debug. I include measured iteration counts plus one in order to detect
45   * regressions. On average Brent-Dekker should use 4..5 iterations for the
46   * default absolute accuracy of 10E-8 for sinus and the quintic function around
47   * zero, and 5..10 iterations for the other zeros.
48   *
49   */
50  public final class BrentSolverTest {
51      @Test
52      public void testSinZero() {
53          // The sinus function is behaved well around the root at pi. The second
54          // order derivative is zero, which means linar approximating methods will
55          // still converge quadratically.
56          UnivariateFunction f = new Sin();
57          double result;
58          UnivariateSolver solver = new BrentSolver();
59          // Somewhat benign interval. The function is monotone.
60          result = solver.solve(100, f, 3, 4);
61          // System.out.println(
62          //    "Root: " + result + " Evaluations: " + solver.getEvaluations());
63          Assert.assertEquals(result, FastMath.PI, solver.getAbsoluteAccuracy());
64          Assert.assertTrue(solver.getEvaluations() <= 7);
65          // Larger and somewhat less benign interval. The function is grows first.
66          result = solver.solve(100, f, 1, 4);
67          // System.out.println(
68          //    "Root: " + result + " Evaluations: " + solver.getEvaluations());
69          Assert.assertEquals(result, FastMath.PI, solver.getAbsoluteAccuracy());
70          Assert.assertTrue(solver.getEvaluations() <= 8);
71      }
72  
73      @Test
74      public void testQuinticZero() {
75          // The quintic function has zeros at 0, +-0.5 and +-1.
76          // Around the root of 0 the function is well behaved, with a second derivative
77          // of zero a 0.
78          // The other roots are less well to find, in particular the root at 1, because
79          // the function grows fast for x>1.
80          // The function has extrema (first derivative is zero) at 0.27195613 and 0.82221643,
81          // intervals containing these values are harder for the solvers.
82          UnivariateFunction f = new QuinticFunction();
83          double result;
84          // Brent-Dekker solver.
85          UnivariateSolver solver = new BrentSolver();
86          // Symmetric bracket around 0. Test whether solvers can handle hitting
87          // the root in the first iteration.
88          result = solver.solve(100, f, -0.2, 0.2);
89          //System.out.println(
90          //    "Root: " + result + " Evaluations: " + solver.getEvaluations());
91          Assert.assertEquals(result, 0, solver.getAbsoluteAccuracy());
92          Assert.assertTrue(solver.getEvaluations() <= 3);
93          // 1 iterations on i586 JDK 1.4.1.
94          // Asymmetric bracket around 0, just for fun. Contains extremum.
95          result = solver.solve(100, f, -0.1, 0.3);
96          //System.out.println(
97          //    "Root: " + result + " Evaluations: " + solver.getEvaluations());
98          Assert.assertEquals(result, 0, solver.getAbsoluteAccuracy());
99          // 5 iterations on i586 JDK 1.4.1.
100         Assert.assertTrue(solver.getEvaluations() <= 7);
101         // Large bracket around 0. Contains two extrema.
102         result = solver.solve(100, f, -0.3, 0.45);
103         //System.out.println(
104         //    "Root: " + result + " Evaluations: " + solver.getEvaluations());
105         Assert.assertEquals(result, 0, solver.getAbsoluteAccuracy());
106         // 6 iterations on i586 JDK 1.4.1.
107         Assert.assertTrue(solver.getEvaluations() <= 8);
108         // Benign bracket around 0.5, function is monotonous.
109         result = solver.solve(100, f, 0.3, 0.7);
110         //System.out.println(
111         //    "Root: " + result + " Evaluations: " + solver.getEvaluations());
112         Assert.assertEquals(result, 0.5, solver.getAbsoluteAccuracy());
113         // 6 iterations on i586 JDK 1.4.1.
114         Assert.assertTrue(solver.getEvaluations() <= 9);
115         // Less benign bracket around 0.5, contains one extremum.
116         result = solver.solve(100, f, 0.2, 0.6);
117         // System.out.println(
118         //    "Root: " + result + " Evaluations: " + solver.getEvaluations());
119         Assert.assertEquals(result, 0.5, solver.getAbsoluteAccuracy());
120         Assert.assertTrue(solver.getEvaluations() <= 10);
121         // Large, less benign bracket around 0.5, contains both extrema.
122         result = solver.solve(100, f, 0.05, 0.95);
123         //System.out.println(
124         //    "Root: " + result + " Evaluations: " + solver.getEvaluations());
125         Assert.assertEquals(result, 0.5, solver.getAbsoluteAccuracy());
126         Assert.assertTrue(solver.getEvaluations() <= 11);
127         // Relatively benign bracket around 1, function is monotonous. Fast growth for x>1
128         // is still a problem.
129         result = solver.solve(100, f, 0.85, 1.25);
130         //System.out.println(
131         //    "Root: " + result + " Evaluations: " + solver.getEvaluations());
132         Assert.assertEquals(result, 1.0, solver.getAbsoluteAccuracy());
133         Assert.assertTrue(solver.getEvaluations() <= 11);
134         // Less benign bracket around 1 with extremum.
135         result = solver.solve(100, f, 0.8, 1.2);
136         //System.out.println(
137         //    "Root: " + result + " Evaluations: " + solver.getEvaluations());
138         Assert.assertEquals(result, 1.0, solver.getAbsoluteAccuracy());
139         Assert.assertTrue(solver.getEvaluations() <= 11);
140         // Large bracket around 1. Monotonous.
141         result = solver.solve(100, f, 0.85, 1.75);
142         //System.out.println(
143         //    "Root: " + result + " Evaluations: " + solver.getEvaluations());
144         Assert.assertEquals(result, 1.0, solver.getAbsoluteAccuracy());
145         Assert.assertTrue(solver.getEvaluations() <= 13);
146         // Large bracket around 1. Interval contains extremum.
147         result = solver.solve(100, f, 0.55, 1.45);
148         //System.out.println(
149         //    "Root: " + result + " Evaluations: " + solver.getEvaluations());
150         Assert.assertEquals(result, 1.0, solver.getAbsoluteAccuracy());
151         Assert.assertTrue(solver.getEvaluations() <= 10);
152         // Very large bracket around 1 for testing fast growth behaviour.
153         result = solver.solve(100, f, 0.85, 5);
154         //System.out.println(
155        //     "Root: " + result + " Evaluations: " + solver.getEvaluations());
156         Assert.assertEquals(result, 1.0, solver.getAbsoluteAccuracy());
157         Assert.assertTrue(solver.getEvaluations() <= 15);
158 
159         try {
160             result = solver.solve(5, f, 0.85, 5);
161             Assert.fail("Expected MathIllegalStateException");
162         } catch (MathIllegalStateException e) {
163             // Expected.
164         }
165     }
166 
167     @Test
168     public void testRootEndpoints() {
169         UnivariateFunction f = new Sin();
170         BrentSolver solver = new BrentSolver();
171 
172         // endpoint is root
173         double result = solver.solve(100, f, FastMath.PI, 4);
174         Assert.assertEquals(FastMath.PI, result, solver.getAbsoluteAccuracy());
175 
176         result = solver.solve(100, f, 3, FastMath.PI);
177         Assert.assertEquals(FastMath.PI, result, solver.getAbsoluteAccuracy());
178 
179         result = solver.solve(100, f, FastMath.PI, 4, 3.5);
180         Assert.assertEquals(FastMath.PI, result, solver.getAbsoluteAccuracy());
181 
182         result = solver.solve(100, f, 3, FastMath.PI, 3.07);
183         Assert.assertEquals(FastMath.PI, result, solver.getAbsoluteAccuracy());
184     }
185 
186     @Test
187     public void testBadEndpoints() {
188         UnivariateFunction f = new Sin();
189         BrentSolver solver = new BrentSolver();
190         try {  // bad interval
191             solver.solve(100, f, 1, -1);
192             Assert.fail("Expecting MathIllegalArgumentException - bad interval");
193         } catch (MathIllegalArgumentException ex) {
194             // expected
195         }
196         try {  // no bracket
197             solver.solve(100, f, 1, 1.5);
198             Assert.fail("Expecting MathIllegalArgumentException - non-bracketing");
199         } catch (MathIllegalArgumentException ex) {
200             // expected
201         }
202         try {  // no bracket
203             solver.solve(100, f, 1, 1.5, 1.2);
204             Assert.fail("Expecting MathIllegalArgumentException - non-bracketing");
205         } catch (MathIllegalArgumentException ex) {
206             // expected
207         }
208     }
209 
210     @Test
211     public void testInitialGuess() {
212         MonitoredFunction f = new MonitoredFunction(new QuinticFunction());
213         BrentSolver solver = new BrentSolver();
214         double result;
215 
216         // no guess
217         result = solver.solve(100, f, 0.6, 7.0);
218         Assert.assertEquals(result, 1.0, solver.getAbsoluteAccuracy());
219         int referenceCallsCount = f.getCallsCount();
220         Assert.assertTrue(referenceCallsCount >= 13);
221 
222         // invalid guess (it *is* a root, but outside of the range)
223         try {
224           result = solver.solve(100, f, 0.6, 7.0, 0.0);
225           Assert.fail("a MathIllegalArgumentException was expected");
226         } catch (MathIllegalArgumentException iae) {
227             // expected behaviour
228         }
229 
230         // bad guess
231         f.setCallsCount(0);
232         result = solver.solve(100, f, 0.6, 7.0, 0.61);
233         Assert.assertEquals(result, 1.0, solver.getAbsoluteAccuracy());
234         Assert.assertTrue(f.getCallsCount() > referenceCallsCount);
235 
236         // good guess
237         f.setCallsCount(0);
238         result = solver.solve(100, f, 0.6, 7.0, 0.999999);
239         Assert.assertEquals(result, 1.0, solver.getAbsoluteAccuracy());
240         Assert.assertTrue(f.getCallsCount() < referenceCallsCount);
241 
242         // perfect guess
243         f.setCallsCount(0);
244         result = solver.solve(100, f, 0.6, 7.0, 1.0);
245         Assert.assertEquals(result, 1.0, solver.getAbsoluteAccuracy());
246         Assert.assertEquals(1, solver.getEvaluations());
247         Assert.assertEquals(1, f.getCallsCount());
248     }
249 
250     @Test
251     public void testMath832() {
252         final UnivariateFunction f = new UnivariateFunction() {
253                 private final UnivariateDifferentiableFunction sqrt = new Sqrt();
254                 private final UnivariateDifferentiableFunction inv = new Inverse();
255                 private final UnivariateDifferentiableFunction func
256                     = FunctionUtils.add(FunctionUtils.multiply(new Constant(1e2), sqrt),
257                                         FunctionUtils.multiply(new Constant(1e6), inv),
258                                         FunctionUtils.multiply(new Constant(1e4),
259                                                                FunctionUtils.compose(inv, sqrt)));
260                 private final DSFactory factory = new DSFactory(1, 1);
261 
262                 public double value(double x) {
263                     return func.value(factory.variable(0, x)).getPartialDerivative(1);
264                 }
265 
266             };
267 
268         BrentSolver solver = new BrentSolver();
269         final double result = solver.solve(99, f, 1, 1e30, 1 + 1e-10);
270         Assert.assertEquals(804.93558250, result, 1e-8);
271     }
272 }