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