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