View Javadoc
1   /*
2    * Licensed to the Hipparchus project 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 Hipparchus project 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  package org.hipparchus.ode.nonstiff;
19  
20  
21  import org.hipparchus.analysis.UnivariateFunction;
22  import org.hipparchus.analysis.polynomials.PolynomialFunction;
23  import org.hipparchus.linear.Array2DRowRealMatrix;
24  import org.junit.jupiter.api.Test;
25  
26  import static org.junit.jupiter.api.Assertions.assertEquals;
27  
28  class AdamsNordsieckTransformerTest {
29  
30      @Test
31      void testPolynomialExtraDerivative() {
32          checkNordsieckStart(new PolynomialFunction(new double[] { 6, 5, 4, 3, 2, 1 }),
33                              5, 0.0, 0.125, 3.2e-16);
34      }
35  
36      @Test
37      void testPolynomialRegular() {
38          checkNordsieckStart(new PolynomialFunction(new double[] { 6, 5, 4, 3, 2, 1 }),
39                              4, 0.0, 0.125, 3.1e-16);
40      }
41  
42      @Test
43      void testPolynomialMissingLastDerivative() {
44          // this test intentionally uses not enough start points,
45          // the Nordsieck vector is therefore not expected to match the exact scaled derivatives
46          checkNordsieckStart(new PolynomialFunction(new double[] { 6, 5, 4, 3, 2, 1 }),
47                              3, 0.0, 0.125, 1.6e-4);
48      }
49  
50      @Test
51      void testTransformExact() {
52          // a 5 steps transformer handles a degree 5 polynomial exactly
53          // the Nordsieck vector holds the full information about the function
54          // transforming the vector from t0 to t0+h or recomputing it from scratch
55          // at t0+h yields the same result
56          checkTransform(new PolynomialFunction(new double[] { 6, 5, 4, 3, 2, 1 }), 5, 2.567e-15);
57      }
58  
59      @Test
60      void testTransformInexact() {
61          // a 4 steps transformer cannot handle a degree 5 polynomial exactly
62          // the Nordsieck vector lacks some high degree information about the function
63          // transforming the vector from t0 to t0+h or recomputing it from scratch
64          // at t0+h yields different results
65          checkTransform(new PolynomialFunction(new double[] { 6, 5, 4, 3, 2, 1 }), 4, 5.658e-4);
66      }
67  
68      private void checkNordsieckStart(final PolynomialFunction polynomial, final int nbSteps, final double t0,
69                                       final double h, final double epsilon) {
70  
71          final AdamsNordsieckTransformer transformer = AdamsNordsieckTransformer.getInstance(nbSteps);
72          PolynomialFunction derivative = polynomial.polynomialDerivative();
73          final Array2DRowRealMatrix nordsieck = start(transformer, nbSteps, t0, h, polynomial, derivative);
74  
75          assertEquals(nbSteps - 1, nordsieck.getRowDimension());
76          double coeff = h;
77          for (int i = 0; i < nordsieck.getRowDimension(); ++i) {
78              coeff *= h / (i + 2);
79              derivative = derivative.polynomialDerivative();
80              assertEquals(derivative.value(t0) * coeff, nordsieck.getEntry(i, 0), epsilon);
81          }
82  
83      }
84  
85      private void checkTransform(final PolynomialFunction polynomial, final int nbSteps, final double expectedError) {
86  
87          final AdamsNordsieckTransformer transformer = AdamsNordsieckTransformer.getInstance(nbSteps);
88          final PolynomialFunction derivative = polynomial.polynomialDerivative();
89  
90          final double t0 = 0.0;
91          final double h  = 0.125;
92          final Array2DRowRealMatrix n0 = start(transformer, nbSteps, t0, h, polynomial, derivative);
93          final Array2DRowRealMatrix n1 = transformer.updateHighOrderDerivativesPhase1(n0);
94          transformer.updateHighOrderDerivativesPhase2(new double[] { h * derivative.value(t0)     },
95                                                       new double[] { h * derivative.value(t0 + h) },
96                                                       n1);
97          final Array2DRowRealMatrix n2 = start(transformer, nbSteps, t0 + h, h, polynomial, derivative);
98  
99          assertEquals(expectedError, n2.subtract(n1).getNorm1(), expectedError * 0.001);
100 
101     }
102 
103     private Array2DRowRealMatrix start(final AdamsNordsieckTransformer transformer, final int nbSteps,
104                                        final double t0, final double h,
105                                        final UnivariateFunction f0, final UnivariateFunction f1) {
106 
107         final int        nbStartPoints = (nbSteps + 3) / 2;
108         final double[]   t             = new double[nbStartPoints];
109         final double[][] y             = new double[nbStartPoints][1];
110         final double[][] yDot          = new double[nbStartPoints][1];
111         for (int i = 0; i < nbStartPoints; ++i) {
112             t[i]       = t0 + i * h;
113             y[i][0]    = f0.value(t[i]);
114             yDot[i][0] = f1.value(t[i]);
115         }
116 
117         return transformer.initializeHighOrderDerivatives(h, t, y, yDot);
118 
119     }
120 
121 }