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