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.ode.nonstiff.interpolators;
24  
25  import java.util.Arrays;
26  
27  import org.hipparchus.CalculusFieldElement;
28  import org.hipparchus.linear.Array2DRowFieldMatrix;
29  import org.hipparchus.ode.FieldEquationsMapper;
30  import org.hipparchus.ode.FieldODEStateAndDerivative;
31  import org.hipparchus.ode.nonstiff.AdamsBashforthFieldIntegrator;
32  import org.hipparchus.ode.nonstiff.AdamsMoultonFieldIntegrator;
33  import org.hipparchus.ode.sampling.AbstractFieldODEStateInterpolator;
34  import org.hipparchus.util.MathArrays;
35  
36  /**
37   * This class implements an interpolator for Adams integrators using Nordsieck representation.
38   *
39   * <p>This interpolator computes dense output around the current point.
40   * The interpolation equation is based on Taylor series formulas.
41   *
42   * @see AdamsBashforthFieldIntegrator
43   * @see AdamsMoultonFieldIntegrator
44   * @param <T> the type of the field elements
45   */
46  
47  public class AdamsFieldStateInterpolator<T extends CalculusFieldElement<T>> extends AbstractFieldODEStateInterpolator<T> {
48  
49      /** Step size used in the first scaled derivative and Nordsieck vector. */
50      private T scalingH;
51  
52      /** Reference state.
53       * <p>Sometimes, the reference state is the same as globalPreviousState,
54       * sometimes it is the same as globalCurrentState, so we use a separate
55       * field to avoid any confusion.
56       * </p>
57       */
58      private final FieldODEStateAndDerivative<T> reference;
59  
60      /** First scaled derivative. */
61      private final T[] scaled;
62  
63      /** Nordsieck vector. */
64      private final Array2DRowFieldMatrix<T> nordsieck;
65  
66      /** Simple constructor.
67       * @param stepSize step size used in the scaled and Nordsieck arrays
68       * @param reference reference state from which Taylor expansion are estimated
69       * @param scaled first scaled derivative
70       * @param nordsieck Nordsieck vector
71       * @param isForward integration direction indicator
72       * @param globalPreviousState start of the global step
73       * @param globalCurrentState end of the global step
74       * @param equationsMapper mapper for ODE equations primary and secondary components
75       */
76      public AdamsFieldStateInterpolator(final T stepSize, final FieldODEStateAndDerivative<T> reference,
77                                         final T[] scaled, final Array2DRowFieldMatrix<T> nordsieck,
78                                         final boolean isForward,
79                                         final FieldODEStateAndDerivative<T> globalPreviousState,
80                                         final FieldODEStateAndDerivative<T> globalCurrentState,
81                                         final FieldEquationsMapper<T> equationsMapper) {
82          this(stepSize, reference, scaled, nordsieck, isForward, globalPreviousState, globalCurrentState,
83                  globalPreviousState, globalCurrentState, equationsMapper);
84      }
85  
86      /** Simple constructor.
87       * @param stepSize step size used in the scaled and Nordsieck arrays
88       * @param reference reference state from which Taylor expansion are estimated
89       * @param scaled first scaled derivative
90       * @param nordsieck Nordsieck vector
91       * @param isForward integration direction indicator
92       * @param globalPreviousState start of the global step
93       * @param globalCurrentState end of the global step
94       * @param softPreviousState start of the restricted step
95       * @param softCurrentState end of the restricted step
96       * @param equationsMapper mapper for ODE equations primary and secondary components
97       */
98      private AdamsFieldStateInterpolator(final T stepSize, final FieldODEStateAndDerivative<T> reference,
99                                          final T[] scaled, final Array2DRowFieldMatrix<T> nordsieck,
100                                         final boolean isForward,
101                                         final FieldODEStateAndDerivative<T> globalPreviousState,
102                                         final FieldODEStateAndDerivative<T> globalCurrentState,
103                                         final FieldODEStateAndDerivative<T> softPreviousState,
104                                         final FieldODEStateAndDerivative<T> softCurrentState,
105                                         final FieldEquationsMapper<T> equationsMapper) {
106         super(isForward, globalPreviousState, globalCurrentState,
107               softPreviousState, softCurrentState, equationsMapper);
108         this.scalingH  = stepSize;
109         this.reference = reference;
110         this.scaled    = scaled.clone();
111         this.nordsieck = new Array2DRowFieldMatrix<>(nordsieck.getData(), false);
112     }
113 
114     /** Create a new instance.
115      * @param newForward integration direction indicator
116      * @param newGlobalPreviousState start of the global step
117      * @param newGlobalCurrentState end of the global step
118      * @param newSoftPreviousState start of the restricted step
119      * @param newSoftCurrentState end of the restricted step
120      * @param newMapper equations mapper for the all equations
121      * @return a new instance
122      */
123     @Override
124     protected AdamsFieldStateInterpolator<T> create(boolean newForward,
125                                                     FieldODEStateAndDerivative<T> newGlobalPreviousState,
126                                                     FieldODEStateAndDerivative<T> newGlobalCurrentState,
127                                                     FieldODEStateAndDerivative<T> newSoftPreviousState,
128                                                     FieldODEStateAndDerivative<T> newSoftCurrentState,
129                                                     FieldEquationsMapper<T> newMapper) {
130         return new AdamsFieldStateInterpolator<>(scalingH, reference, scaled, nordsieck,
131                                                   newForward,
132                                                   newGlobalPreviousState, newGlobalCurrentState,
133                                                   newSoftPreviousState, newSoftCurrentState,
134                                                   newMapper);
135 
136     }
137 
138     /** Get the first scaled derivative.
139      * @return first scaled derivative
140      */
141     public T[] getScaled() {
142         return scaled.clone();
143     }
144 
145     /** Get the Nordsieck vector.
146      * @return Nordsieck vector
147      */
148     public Array2DRowFieldMatrix<T> getNordsieck() {
149         return new Array2DRowFieldMatrix<>(nordsieck.getData());
150     }
151 
152     /** {@inheritDoc} */
153     @Override
154     protected FieldODEStateAndDerivative<T> computeInterpolatedStateAndDerivatives(final FieldEquationsMapper<T> equationsMapper,
155                                                                                    final T time, final T theta,
156                                                                                    final T thetaH, final T oneMinusThetaH) {
157         return taylor(equationsMapper, reference, time, scalingH, scaled, nordsieck);
158     }
159 
160     /** Estimate state by applying Taylor formula.
161      * @param equationsMapper mapper for ODE equations primary and secondary components
162      * @param reference reference state
163      * @param time time at which state must be estimated
164      * @param stepSize step size used in the scaled and Nordsieck arrays
165      * @param scaled first scaled derivative
166      * @param nordsieck Nordsieck vector
167      * @return estimated state
168      * @param <S> the type of the field elements
169      */
170     public static <S extends CalculusFieldElement<S>> FieldODEStateAndDerivative<S> taylor(final FieldEquationsMapper<S> equationsMapper,
171                                                                                            final FieldODEStateAndDerivative<S> reference,
172                                                                                            final S time, final S stepSize,
173                                                                                            final S[] scaled,
174                                                                                            final Array2DRowFieldMatrix<S> nordsieck) {
175 
176         final S x = time.subtract(reference.getTime());
177         final S normalizedAbscissa = x.divide(stepSize);
178 
179         S[] stateVariation = MathArrays.buildArray(time.getField(), scaled.length);
180         Arrays.fill(stateVariation, time.getField().getZero());
181         S[] estimatedDerivatives = MathArrays.buildArray(time.getField(), scaled.length);
182         Arrays.fill(estimatedDerivatives, time.getField().getZero());
183 
184         // apply Taylor formula from high order to low order,
185         // for the sake of numerical accuracy
186         final S[][] nData = nordsieck.getDataRef();
187         for (int i = nData.length - 1; i >= 0; --i) {
188             final int order = i + 2;
189             final S[] nDataI = nData[i];
190             final S power = normalizedAbscissa.pow(order);
191             for (int j = 0; j < nDataI.length; ++j) {
192                 final S d = nDataI[j].multiply(power);
193                 stateVariation[j]          = stateVariation[j].add(d);
194                 estimatedDerivatives[j] = estimatedDerivatives[j].add(d.multiply(order));
195             }
196         }
197 
198         S[] estimatedState = reference.getCompleteState();
199         for (int j = 0; j < stateVariation.length; ++j) {
200             stateVariation[j] = stateVariation[j].add(scaled[j].multiply(normalizedAbscissa));
201             estimatedState[j] = estimatedState[j].add(stateVariation[j]);
202             estimatedDerivatives[j] =
203                 estimatedDerivatives[j].add(scaled[j].multiply(normalizedAbscissa)).divide(x);
204         }
205 
206         return equationsMapper.mapStateAndDerivative(time, estimatedState, estimatedDerivatives);
207 
208     }
209 
210 }