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