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