001/*
002 * Licensed to the Apache Software Foundation (ASF) under one or more
003 * contributor license agreements.  See the NOTICE file distributed with
004 * this work for additional information regarding copyright ownership.
005 * The ASF licenses this file to You under the Apache License, Version 2.0
006 * (the "License"); you may not use this file except in compliance with
007 * the License.  You may obtain a copy of the License at
008 *
009 *      http://www.apache.org/licenses/LICENSE-2.0
010 *
011 * Unless required by applicable law or agreed to in writing, software
012 * distributed under the License is distributed on an "AS IS" BASIS,
013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014 * See the License for the specific language governing permissions and
015 * limitations under the License.
016 */
017
018package org.apache.commons.math4.ode.nonstiff;
019
020
021import org.apache.commons.math4.Field;
022import org.apache.commons.math4.RealFieldElement;
023import org.apache.commons.math4.exception.DimensionMismatchException;
024import org.apache.commons.math4.exception.MaxCountExceededException;
025import org.apache.commons.math4.exception.NoBracketingException;
026import org.apache.commons.math4.exception.NumberIsTooSmallException;
027import org.apache.commons.math4.ode.AbstractFieldIntegrator;
028import org.apache.commons.math4.ode.FieldEquationsMapper;
029import org.apache.commons.math4.ode.FieldExpandableODE;
030import org.apache.commons.math4.ode.FirstOrderFieldDifferentialEquations;
031import org.apache.commons.math4.ode.FieldODEState;
032import org.apache.commons.math4.ode.FieldODEStateAndDerivative;
033import org.apache.commons.math4.util.MathArrays;
034
035/**
036 * This class implements the common part of all fixed step Runge-Kutta
037 * integrators for Ordinary Differential Equations.
038 *
039 * <p>These methods are explicit Runge-Kutta methods, their Butcher
040 * arrays are as follows :
041 * <pre>
042 *    0  |
043 *   c2  | a21
044 *   c3  | a31  a32
045 *   ... |        ...
046 *   cs  | as1  as2  ...  ass-1
047 *       |--------------------------
048 *       |  b1   b2  ...   bs-1  bs
049 * </pre>
050 *
051 * @see EulerFieldIntegrator
052 * @see ClassicalRungeKuttaFieldIntegrator
053 * @see GillFieldIntegrator
054 * @see MidpointFieldIntegrator
055 * @param <T> the type of the field elements
056 * @since 3.6
057 */
058
059public abstract class RungeKuttaFieldIntegrator<T extends RealFieldElement<T>>
060    extends AbstractFieldIntegrator<T>
061    implements FieldButcherArrayProvider<T> {
062
063    /** Time steps from Butcher array (without the first zero). */
064    private final T[] c;
065
066    /** Internal weights from Butcher array (without the first empty row). */
067    private final T[][] a;
068
069    /** External weights for the high order method from Butcher array. */
070    private final T[] b;
071
072    /** Integration step. */
073    private final T step;
074
075    /** Simple constructor.
076     * Build a Runge-Kutta integrator with the given
077     * step. The default step handler does nothing.
078     * @param field field to which the time and state vector elements belong
079     * @param name name of the method
080     * @param step integration step
081     */
082    protected RungeKuttaFieldIntegrator(final Field<T> field, final String name, final T step) {
083        super(field, name);
084        this.c    = getC();
085        this.a    = getA();
086        this.b    = getB();
087        this.step = step.abs();
088    }
089
090    /** Create a fraction.
091     * @param p numerator
092     * @param q denominator
093     * @return p/q computed in the instance field
094     */
095    protected T fraction(final int p, final int q) {
096        return getField().getZero().add(p).divide(q);
097    }
098
099    /** Create an interpolator.
100     * @param forward integration direction indicator
101     * @param yDotK slopes at the intermediate points
102     * @param globalPreviousState start of the global step
103     * @param globalCurrentState end of the global step
104     * @param mapper equations mapper for the all equations
105     * @return external weights for the high order method from Butcher array
106     */
107    protected abstract RungeKuttaFieldStepInterpolator<T> createInterpolator(boolean forward, T[][] yDotK,
108                                                                             final FieldODEStateAndDerivative<T> globalPreviousState,
109                                                                             final FieldODEStateAndDerivative<T> globalCurrentState,
110                                                                             FieldEquationsMapper<T> mapper);
111
112    /** {@inheritDoc} */
113    @Override
114    public FieldODEStateAndDerivative<T> integrate(final FieldExpandableODE<T> equations,
115                                                   final FieldODEState<T> initialState, final T finalTime)
116        throws NumberIsTooSmallException, DimensionMismatchException,
117        MaxCountExceededException, NoBracketingException {
118
119        sanityChecks(initialState, finalTime);
120        final T   t0 = initialState.getTime();
121        final T[] y0 = equations.getMapper().mapState(initialState);
122        setStepStart(initIntegration(equations, t0, y0, finalTime));
123        final boolean forward = finalTime.subtract(initialState.getTime()).getReal() > 0;
124
125        // create some internal working arrays
126        final int   stages = c.length + 1;
127        T[]         y      = y0;
128        final T[][] yDotK  = MathArrays.buildArray(getField(), stages, -1);
129        final T[]   yTmp   = MathArrays.buildArray(getField(), y0.length);
130
131        // set up integration control objects
132        if (forward) {
133            if (getStepStart().getTime().add(step).subtract(finalTime).getReal() >= 0) {
134                setStepSize(finalTime.subtract(getStepStart().getTime()));
135            } else {
136                setStepSize(step);
137            }
138        } else {
139            if (getStepStart().getTime().subtract(step).subtract(finalTime).getReal() <= 0) {
140                setStepSize(finalTime.subtract(getStepStart().getTime()));
141            } else {
142                setStepSize(step.negate());
143            }
144        }
145
146        // main integration loop
147        setIsLastStep(false);
148        do {
149
150            // first stage
151            y        = equations.getMapper().mapState(getStepStart());
152            yDotK[0] = equations.getMapper().mapDerivative(getStepStart());
153
154            // next stages
155            for (int k = 1; k < stages; ++k) {
156
157                for (int j = 0; j < y0.length; ++j) {
158                    T sum = yDotK[0][j].multiply(a[k-1][0]);
159                    for (int l = 1; l < k; ++l) {
160                        sum = sum.add(yDotK[l][j].multiply(a[k-1][l]));
161                    }
162                    yTmp[j] = y[j].add(getStepSize().multiply(sum));
163                }
164
165                yDotK[k] = computeDerivatives(getStepStart().getTime().add(getStepSize().multiply(c[k-1])), yTmp);
166
167            }
168
169            // estimate the state at the end of the step
170            for (int j = 0; j < y0.length; ++j) {
171                T sum = yDotK[0][j].multiply(b[0]);
172                for (int l = 1; l < stages; ++l) {
173                    sum = sum.add(yDotK[l][j].multiply(b[l]));
174                }
175                yTmp[j] = y[j].add(getStepSize().multiply(sum));
176            }
177            final T stepEnd   = getStepStart().getTime().add(getStepSize());
178            final T[] yDotTmp = computeDerivatives(stepEnd, yTmp);
179            final FieldODEStateAndDerivative<T> stateTmp = new FieldODEStateAndDerivative<>(stepEnd, yTmp, yDotTmp);
180
181            // discrete events handling
182            System.arraycopy(yTmp, 0, y, 0, y0.length);
183            setStepStart(acceptStep(createInterpolator(forward, yDotK, getStepStart(), stateTmp, equations.getMapper()),
184                                    finalTime));
185
186            if (!isLastStep()) {
187
188                // stepsize control for next step
189                final T  nextT      = getStepStart().getTime().add(getStepSize());
190                final boolean nextIsLast = forward ?
191                                           (nextT.subtract(finalTime).getReal() >= 0) :
192                                           (nextT.subtract(finalTime).getReal() <= 0);
193                if (nextIsLast) {
194                    setStepSize(finalTime.subtract(getStepStart().getTime()));
195                }
196            }
197
198        } while (!isLastStep());
199
200        final FieldODEStateAndDerivative<T> finalState = getStepStart();
201        setStepStart(null);
202        setStepSize(null);
203        return finalState;
204
205    }
206
207    /** Fast computation of a single step of ODE integration.
208     * <p>This method is intended for the limited use case of
209     * very fast computation of only one step without using any of the
210     * rich features of general integrators that may take some time
211     * to set up (i.e. no step handlers, no events handlers, no additional
212     * states, no interpolators, no error control, no evaluations count,
213     * no sanity checks ...). It handles the strict minimum of computation,
214     * so it can be embedded in outer loops.</p>
215     * <p>
216     * This method is <em>not</em> used at all by the {@link #integrate(FieldExpandableODE,
217     * FieldODEState, RealFieldElement)} method. It also completely ignores the step set at
218     * construction time, and uses only a single step to go from {@code t0} to {@code t}.
219     * </p>
220     * <p>
221     * As this method does not use any of the state-dependent features of the integrator,
222     * it should be reasonably thread-safe <em>if and only if</em> the provided differential
223     * equations are themselves thread-safe.
224     * </p>
225     * @param equations differential equations to integrate
226     * @param t0 initial time
227     * @param y0 initial value of the state vector at t0
228     * @param t target time for the integration
229     * (can be set to a value smaller than {@code t0} for backward integration)
230     * @return state vector at {@code t}
231     */
232    public T[] singleStep(final FirstOrderFieldDifferentialEquations<T> equations,
233                          final T t0, final T[] y0, final T t) {
234
235        // create some internal working arrays
236        final T[] y       = y0.clone();
237        final int stages  = c.length + 1;
238        final T[][] yDotK = MathArrays.buildArray(getField(), stages, -1);
239        final T[] yTmp    = y0.clone();
240
241        // first stage
242        final T h = t.subtract(t0);
243        yDotK[0] = equations.computeDerivatives(t0, y);
244
245        // next stages
246        for (int k = 1; k < stages; ++k) {
247
248            for (int j = 0; j < y0.length; ++j) {
249                T sum = yDotK[0][j].multiply(a[k-1][0]);
250                for (int l = 1; l < k; ++l) {
251                    sum = sum.add(yDotK[l][j].multiply(a[k-1][l]));
252                }
253                yTmp[j] = y[j].add(h.multiply(sum));
254            }
255
256            yDotK[k] = equations.computeDerivatives(t0.add(h.multiply(c[k-1])), yTmp);
257
258        }
259
260        // estimate the state at the end of the step
261        for (int j = 0; j < y0.length; ++j) {
262            T sum = yDotK[0][j].multiply(b[0]);
263            for (int l = 1; l < stages; ++l) {
264                sum = sum.add(yDotK[l][j].multiply(b[l]));
265            }
266            y[j] = y[j].add(h.multiply(sum));
267        }
268
269        return y;
270
271    }
272
273}