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 org.hipparchus.CalculusFieldElement;
26  import org.hipparchus.analysis.CalculusFieldUnivariateFunction;
27  import org.hipparchus.analysis.FieldUnivariateFunction;
28  import org.hipparchus.dfp.Dfp;
29  import org.hipparchus.dfp.DfpField;
30  import org.hipparchus.dfp.DfpMath;
31  import org.hipparchus.exception.MathIllegalArgumentException;
32  import org.hipparchus.exception.MathRuntimeException;
33  import org.hipparchus.util.FastMath;
34  import org.junit.jupiter.api.BeforeEach;
35  import org.junit.jupiter.api.Test;
36  
37  import static org.junit.jupiter.api.Assertions.assertEquals;
38  import static org.junit.jupiter.api.Assertions.assertThrows;
39  import static org.junit.jupiter.api.Assertions.assertTrue;
40  
41  /**
42   * Test case for {@link FieldBracketingNthOrderBrentSolver bracketing n<sup>th</sup> order Brent} solver.
43   *
44   */
45  final class FieldBracketingNthOrderBrentSolverTest {
46  
47      @Test
48      void testInsufficientOrder3() {
49          assertThrows(MathIllegalArgumentException.class, () -> {
50              new FieldBracketingNthOrderBrentSolver<Dfp>(relativeAccuracy,
51                                                          absoluteAccuracy,
52                                                          functionValueAccuracy,
53                                                          1);
54          });
55      }
56  
57      @Test
58      void testConstructorOK() {
59          FieldBracketingNthOrderBrentSolver<Dfp> solver =
60                  new FieldBracketingNthOrderBrentSolver<Dfp>(relativeAccuracy, absoluteAccuracy,
61                                                              functionValueAccuracy, 2);
62          assertEquals(2, solver.getMaximalOrder());
63      }
64  
65      @Test
66      void testConvergenceOnFunctionAccuracy() {
67          FieldBracketingNthOrderBrentSolver<Dfp> solver =
68                  new FieldBracketingNthOrderBrentSolver<Dfp>(relativeAccuracy, absoluteAccuracy,
69                                                              field.newDfp(1.0e-20), 20);
70          FieldUnivariateFunction f = new FieldUnivariateFunction() {
71              public <T extends CalculusFieldElement<T>> T value(T x) {
72                  T one     = x.getField().getOne();
73                  T oneHalf = one.divide(2);
74                  T xMo     = x.subtract(one);
75                  T xMh     = x.subtract(oneHalf);
76                  T xPh     = x.add(oneHalf);
77                  T xPo     = x.add(one);
78                  return xMo.multiply(xMh).multiply(x).multiply(xPh).multiply(xPo);
79              }
80          };
81  
82          Dfp result = solver.solve(20, f.toCalculusFieldUnivariateFunction(field), field.newDfp(0.2), field.newDfp(0.9),
83                                    field.newDfp(0.4), AllowedSolution.BELOW_SIDE);
84          assertTrue(f.value(result).abs().lessThan(solver.getFunctionValueAccuracy()));
85          assertTrue(f.value(result).negativeOrNull());
86          assertTrue(result.subtract(field.newDfp(0.5)).subtract(solver.getAbsoluteAccuracy()).positiveOrNull());
87          result = solver.solve(20, f.toCalculusFieldUnivariateFunction(field), field.newDfp(-0.9), field.newDfp(-0.2),
88                                field.newDfp(-0.4), AllowedSolution.ABOVE_SIDE);
89          assertTrue(f.value(result).abs().lessThan(solver.getFunctionValueAccuracy()));
90          assertTrue(f.value(result).positiveOrNull());
91          assertTrue(result.add(field.newDfp(0.5)).subtract(solver.getAbsoluteAccuracy()).negativeOrNull());
92      }
93  
94      @Test
95      void testToleranceLessThanUlp() {
96          // function that is never zero
97          Dfp zero = field.getZero();
98          Dfp one = field.getOne();
99          CalculusFieldUnivariateFunction<Dfp> f = (x) -> x.getReal() <= 2.1 ? one.negate(): one;
100         // tolerance less than 1 ulp(x)
101         FieldBracketingNthOrderBrentSolver<Dfp> solver =
102                 new FieldBracketingNthOrderBrentSolver<>(zero, field.newDfp(1e-55), zero, 5);
103 
104         // make sure it doesn't throw a maxIterations exception
105         Dfp result = solver.solve(200, f, zero, zero.add(5.0), AllowedSolution.LEFT_SIDE);
106         double difference = field.newDfp(2.1).subtract(result).abs().getReal();
107         assertTrue( difference < FastMath.ulp(2.1), "difference: " + difference);
108     }
109 
110     @Test
111     void testNeta() {
112 
113         // the following test functions come from Beny Neta's paper:
114         // "Several New Methods for solving Equations"
115         // intern J. Computer Math Vol 23 pp 265-282
116         // available here: http://www.math.nps.navy.mil/~bneta/SeveralNewMethods.PDF
117         for (AllowedSolution allowed : AllowedSolution.values()) {
118             check(new CalculusFieldUnivariateFunction<Dfp>() {
119                 public Dfp value(Dfp x) {
120                     return DfpMath.sin(x).subtract(x.divide(2));
121                 }
122             }, 200, -2.0, 2.0, allowed);
123 
124             check(new CalculusFieldUnivariateFunction<Dfp>() {
125                 public Dfp value(Dfp x) {
126                     return DfpMath.pow(x, 5).add(x).subtract(field.newDfp(10000));
127                 }
128             }, 200, -5.0, 10.0, allowed);
129 
130             check(new CalculusFieldUnivariateFunction<Dfp>() {
131                 public Dfp value(Dfp x) {
132                     return x.sqrt().subtract(field.getOne().divide(x)).subtract(field.newDfp(3));
133                 }
134             }, 200, 0.001, 10.0, allowed);
135 
136             check(new CalculusFieldUnivariateFunction<Dfp>() {
137                 public Dfp value(Dfp x) {
138                     return DfpMath.exp(x).add(x).subtract(field.newDfp(20));
139                 }
140             }, 200, -5.0, 5.0, allowed);
141 
142             check(new CalculusFieldUnivariateFunction<Dfp>() {
143                 public Dfp value(Dfp x) {
144                     return DfpMath.log(x).add(x.sqrt()).subtract(field.newDfp(5));
145                 }
146             }, 200, 0.001, 10.0, allowed);
147 
148             check(new CalculusFieldUnivariateFunction<Dfp>() {
149                 public Dfp value(Dfp x) {
150                     return x.subtract(field.getOne()).multiply(x).multiply(x).subtract(field.getOne());
151                 }
152             }, 200, -0.5, 1.5, allowed);
153         }
154 
155     }
156 
157     private void check(CalculusFieldUnivariateFunction<Dfp> f, int maxEval, double min, double max,
158                        AllowedSolution allowedSolution) {
159         FieldBracketingNthOrderBrentSolver<Dfp> solver =
160                 new FieldBracketingNthOrderBrentSolver<Dfp>(relativeAccuracy, absoluteAccuracy,
161                                                      functionValueAccuracy, 20);
162         Dfp xResult = solver.solve(maxEval, f, field.newDfp(min), field.newDfp(max),
163                                    allowedSolution);
164         Dfp yResult = f.value(xResult);
165         switch (allowedSolution) {
166         case ANY_SIDE :
167             assertTrue(yResult.abs().lessThan(functionValueAccuracy.multiply(2)));
168             break;
169         case LEFT_SIDE : {
170             boolean increasing = f.value(xResult).add(absoluteAccuracy).greaterThan(yResult);
171             assertTrue(increasing ? yResult.negativeOrNull() : yResult.positiveOrNull());
172             break;
173         }
174         case RIGHT_SIDE : {
175             boolean increasing = f.value(xResult).add(absoluteAccuracy).greaterThan(yResult);
176             assertTrue(increasing ? yResult.positiveOrNull() : yResult.negativeOrNull());
177             break;
178         }
179         case BELOW_SIDE :
180             assertTrue(yResult.negativeOrNull());
181             break;
182         case ABOVE_SIDE :
183             assertTrue(yResult.positiveOrNull());
184             break;
185         default :
186             // this should never happen
187             throw new MathRuntimeException(null);
188         }
189     }
190 
191     @BeforeEach
192     void setUp() {
193         field                 = new DfpField(50);
194         absoluteAccuracy      = field.newDfp(1.0e-45);
195         relativeAccuracy      = field.newDfp(1.0e-45);
196         functionValueAccuracy = field.newDfp(1.0e-45);
197     }
198 
199     private DfpField field;
200     private Dfp      absoluteAccuracy;
201     private Dfp      relativeAccuracy;
202     private Dfp      functionValueAccuracy;
203 
204 }