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 }