Repository: commons-math Updated Branches: refs/heads/field-ode 8949c0b99 -> 7ec2a6342
Intermediate level implementations of fixed-step Rung-Kutta methods. Project: http://git-wip-us.apache.org/repos/asf/commons-math/repo Commit: http://git-wip-us.apache.org/repos/asf/commons-math/commit/7ec2a634 Tree: http://git-wip-us.apache.org/repos/asf/commons-math/tree/7ec2a634 Diff: http://git-wip-us.apache.org/repos/asf/commons-math/diff/7ec2a634 Branch: refs/heads/field-ode Commit: 7ec2a6342bb9c78680d58e6f35c579e52f730a2c Parents: 8949c0b Author: Luc Maisonobe <l...@apache.org> Authored: Sun Nov 15 10:37:04 2015 +0100 Committer: Luc Maisonobe <l...@apache.org> Committed: Sun Nov 15 10:37:04 2015 +0100 ---------------------------------------------------------------------- .../ode/nonstiff/RungeKuttaFieldIntegrator.java | 271 +++++++++++++++++++ .../RungeKuttaFieldStepInterpolator.java | 144 ++++++++++ 2 files changed, 415 insertions(+) ---------------------------------------------------------------------- http://git-wip-us.apache.org/repos/asf/commons-math/blob/7ec2a634/src/main/java/org/apache/commons/math3/ode/nonstiff/RungeKuttaFieldIntegrator.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/ode/nonstiff/RungeKuttaFieldIntegrator.java b/src/main/java/org/apache/commons/math3/ode/nonstiff/RungeKuttaFieldIntegrator.java new file mode 100644 index 0000000..577687b --- /dev/null +++ b/src/main/java/org/apache/commons/math3/ode/nonstiff/RungeKuttaFieldIntegrator.java @@ -0,0 +1,271 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +package org.apache.commons.math3.ode.nonstiff; + + +import org.apache.commons.math3.Field; +import org.apache.commons.math3.RealFieldElement; +import org.apache.commons.math3.exception.DimensionMismatchException; +import org.apache.commons.math3.exception.MaxCountExceededException; +import org.apache.commons.math3.exception.NoBracketingException; +import org.apache.commons.math3.exception.NumberIsTooSmallException; +import org.apache.commons.math3.ode.AbstractFieldIntegrator; +import org.apache.commons.math3.ode.FieldExpandableODE; +import org.apache.commons.math3.ode.FieldFirstOrderDifferentialEquations; +import org.apache.commons.math3.ode.FieldODEState; +import org.apache.commons.math3.ode.FieldODEStateAndDerivative; +import org.apache.commons.math3.util.MathArrays; + +/** + * This class implements the common part of all fixed step Runge-Kutta + * integrators for Ordinary Differential Equations. + * + * <p>These methods are explicit Runge-Kutta methods, their Butcher + * arrays are as follows : + * <pre> + * 0 | + * c2 | a21 + * c3 | a31 a32 + * ... | ... + * cs | as1 as2 ... ass-1 + * |-------------------------- + * | b1 b2 ... bs-1 bs + * </pre> + * </p> + * + * @see EulerFieldIntegrator + * @see ClassicalRungeKuttaFieldIntegrator + * @see GillFieldIntegrator + * @see MidpointFieldIntegrator + * @param <T> the type of the field elements + * @since 3.6 + */ + +public abstract class RungeKuttaFieldIntegrator<T extends RealFieldElement<T>> + extends AbstractFieldIntegrator<T> { + + /** Time steps from Butcher array (without the first zero). */ + private final double[] c; + + /** Internal weights from Butcher array (without the first empty row). */ + private final double[][] a; + + /** External weights for the high order method from Butcher array. */ + private final double[] b; + + /** Prototype of the step interpolator. */ + private final RungeKuttaFieldStepInterpolator<T> prototype; + + /** Integration step. */ + private final T step; + + /** Simple constructor. + * Build a Runge-Kutta integrator with the given + * step. The default step handler does nothing. + * @param field field to which the time and state vector elements belong + * @param name name of the method + * @param c time steps from Butcher array (without the first zero) + * @param a internal weights from Butcher array (without the first empty row) + * @param b propagation weights for the high order method from Butcher array + * @param prototype prototype of the step interpolator to use + * @param step integration step + */ + protected RungeKuttaFieldIntegrator(final Field<T> field, final String name, + final double[] c, final double[][] a, final double[] b, + final RungeKuttaFieldStepInterpolator<T> prototype, + final T step) { + super(field, name); + this.c = c; + this.a = a; + this.b = b; + this.prototype = prototype; + this.step = step.abs(); + } + + /** {@inheritDoc} */ + @Override + public FieldODEStateAndDerivative<T> integrate(final FieldExpandableODE<T> equations, + final FieldODEState<T> initialState, final T finalTime) + throws NumberIsTooSmallException, DimensionMismatchException, + MaxCountExceededException, NoBracketingException { + + sanityChecks(initialState, finalTime); + final T t0 = initialState.getTime(); + final T[] y0 = equations.getMapper().mapState(initialState); + stepStart = initIntegration(equations, t0, y0, finalTime); + final boolean forward = finalTime.subtract(initialState.getTime()).getReal() > 0; + + // create some internal working arrays + final int stages = c.length + 1; + T[] y = y0; + final T[][] yDotK = MathArrays.buildArray(getField(), stages, -1); + final T[] yTmp = MathArrays.buildArray(getField(), y0.length); + + // set up an interpolator sharing the integrator arrays + final RungeKuttaFieldStepInterpolator<T> interpolator = (RungeKuttaFieldStepInterpolator<T>) prototype.copy(); + interpolator.reinitialize(this, y0, yDotK, forward, equations.getMapper()); + interpolator.storeState(stepStart); + + // set up integration control objects + if (forward) { + if (stepStart.getTime().add(step).subtract(finalTime).getReal() >= 0) { + stepSize = finalTime.subtract(stepStart.getTime()); + } else { + stepSize = step; + } + } else { + if (stepStart.getTime().subtract(step).subtract(finalTime).getReal() <= 0) { + stepSize = finalTime.subtract(stepStart.getTime()); + } else { + stepSize = step.negate(); + } + } + + // main integration loop + isLastStep = false; + do { + + interpolator.shift(); + + // first stage + yDotK[0] = stepStart.getDerivative(); + + // next stages + for (int k = 1; k < stages; ++k) { + + for (int j = 0; j < y0.length; ++j) { + T sum = yDotK[0][j].multiply(a[k-1][0]); + for (int l = 1; l < k; ++l) { + sum = sum.add(yDotK[l][j].multiply(a[k-1][l])); + } + yTmp[j] = y[j].add(stepSize.multiply(sum)); + } + + yDotK[k] = computeDerivatives(stepStart.getTime().add(stepSize.multiply(c[k-1])), yTmp); + + } + + // estimate the state at the end of the step + for (int j = 0; j < y0.length; ++j) { + T sum = yDotK[0][j].multiply(b[0]); + for (int l = 1; l < stages; ++l) { + sum = sum.add(yDotK[l][j].multiply(b[l])); + } + yTmp[j] = y[j].add(stepSize.multiply(sum)); + } + final T stepEnd = stepStart.getTime().add(stepSize); + final T[] yDotTmp = computeDerivatives(stepEnd, yTmp); + final FieldODEStateAndDerivative<T> stateTmp = new FieldODEStateAndDerivative<T>(stepEnd, yTmp, yDotTmp); + + // discrete events handling + interpolator.storeState(stateTmp); + System.arraycopy(yTmp, 0, y, 0, y0.length); + System.arraycopy(yDotK[stages - 1], 0, yDotTmp, 0, y0.length); + stepStart = acceptStep(interpolator, finalTime); + + if (!isLastStep) { + + // prepare next step + interpolator.storeState(stepStart); + + // stepsize control for next step + final T nextT = stepStart.getTime().add(stepSize); + final boolean nextIsLast = forward ? + (nextT.subtract(finalTime).getReal() >= 0) : + (nextT.subtract(finalTime).getReal() <= 0); + if (nextIsLast) { + stepSize = finalTime.subtract(stepStart.getTime()); + } + } + + } while (!isLastStep); + + final FieldODEStateAndDerivative<T> finalState = stepStart; + stepStart = null; + stepSize = null; + return finalState; + + } + + /** Fast computation of a single step of ODE integration. + * <p>This method is intended for the limited use case of + * very fast computation of only one step without using any of the + * rich features of general integrators that may take some time + * to set up (i.e. no step handlers, no events handlers, no additional + * states, no interpolators, no error control, no evaluations count, + * no sanity checks ...). It handles the strict minimum of computation, + * so it can be embedded in outer loops.</p> + * <p> + * This method is <em>not</em> used at all by the {@link #integrate(FieldExpandableODE, + * FieldODEState, RealFieldElement)} method. It also completely ignores the step set at + * construction time, and uses only a single step to go from {@code t0} to {@code t}. + * </p> + * <p> + * As this method does not use any of the state-dependent features of the integrator, + * it should be reasonably thread-safe <em>if and only if</em> the provided differential + * equations are themselves thread-safe. + * </p> + * @param equations differential equations to integrate + * @param t0 initial time + * @param y0 initial value of the state vector at t0 + * @param t target time for the integration + * (can be set to a value smaller than {@code t0} for backward integration) + * @return state vector at {@code t} + */ + public T[] singleStep(final FieldFirstOrderDifferentialEquations<T> equations, + final T t0, final T[] y0, final T t) { + + // create some internal working arrays + final T[] y = y0.clone(); + final int stages = c.length + 1; + final T[][] yDotK = MathArrays.buildArray(getField(), stages, -1); + final T[] yTmp = y0.clone(); + + // first stage + final T h = t.subtract(t0); + yDotK[0] = equations.computeDerivatives(t0, y); + + // next stages + for (int k = 1; k < stages; ++k) { + + for (int j = 0; j < y0.length; ++j) { + T sum = yDotK[0][j].multiply(a[k-1][0]); + for (int l = 1; l < k; ++l) { + sum = sum.add(yDotK[l][j].multiply(a[k-1][l])); + } + yTmp[j] = y[j].add(h.multiply(sum)); + } + + yDotK[k] = equations.computeDerivatives(t0.add(h.multiply(c[k-1])), yTmp); + + } + + // estimate the state at the end of the step + for (int j = 0; j < y0.length; ++j) { + T sum = yDotK[0][j].multiply(b[0]); + for (int l = 1; l < stages; ++l) { + sum = sum.add(yDotK[l][j].multiply(b[l])); + } + y[j] = y[j].add(h.multiply(sum)); + } + + return y; + + } + +} http://git-wip-us.apache.org/repos/asf/commons-math/blob/7ec2a634/src/main/java/org/apache/commons/math3/ode/nonstiff/RungeKuttaFieldStepInterpolator.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/ode/nonstiff/RungeKuttaFieldStepInterpolator.java b/src/main/java/org/apache/commons/math3/ode/nonstiff/RungeKuttaFieldStepInterpolator.java new file mode 100644 index 0000000..29fbecb --- /dev/null +++ b/src/main/java/org/apache/commons/math3/ode/nonstiff/RungeKuttaFieldStepInterpolator.java @@ -0,0 +1,144 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +package org.apache.commons.math3.ode.nonstiff; + +import org.apache.commons.math3.RealFieldElement; +import org.apache.commons.math3.ode.AbstractFieldIntegrator; +import org.apache.commons.math3.ode.FieldEquationsMapper; +import org.apache.commons.math3.ode.sampling.AbstractFieldStepInterpolator; +import org.apache.commons.math3.util.MathArrays; + +/** This class represents an interpolator over the last step during an + * ODE integration for Runge-Kutta and embedded Runge-Kutta integrators. + * + * @see RungeKuttaFieldIntegrator + * @see EmbeddedRungeKuttaFieldIntegrator + * + * @param <T> the type of the field elements + * @since 3.6 + */ + +abstract class RungeKuttaFieldStepInterpolator<T extends RealFieldElement<T>> + extends AbstractFieldStepInterpolator<T> { + + /** Previous state. */ + protected T[] previousState; + + /** Slopes at the intermediate points */ + protected T[][] yDotK; + + /** Reference to the integrator. */ + protected AbstractFieldIntegrator<T> integrator; + + /** Simple constructor. + * This constructor builds an instance that is not usable yet, the + * {@link #reinitialize} method should be called before using the + * instance in order to initialize the internal arrays. This + * constructor is used only in order to delay the initialization in + * some cases. The {@link RungeKuttaIntegrator} and {@link + * EmbeddedRungeKuttaIntegrator} classes use the prototyping design + * pattern to create the step interpolators by cloning an + * uninitialized model and latter initializing the copy. + */ + protected RungeKuttaFieldStepInterpolator() { + previousState = null; + yDotK = null; + integrator = null; + } + + /** Copy constructor. + + * <p>The copied interpolator should have been finalized before the + * copy, otherwise the copy will not be able to perform correctly any + * interpolation and will throw a {@link NullPointerException} + * later. Since we don't want this constructor to throw the + * exceptions finalization may involve and since we don't want this + * method to modify the state of the copied interpolator, + * finalization is <strong>not</strong> done automatically, it + * remains under user control.</p> + + * <p>The copy is a deep copy: its arrays are separated from the + * original arrays of the instance.</p> + + * @param interpolator interpolator to copy from. + + */ + RungeKuttaFieldStepInterpolator(final RungeKuttaFieldStepInterpolator<T> interpolator) { + + super(interpolator); + + if (interpolator.currentState != null) { + + previousState = interpolator.previousState.clone(); + + yDotK = MathArrays.buildArray(interpolator.integrator.getField(), + interpolator.yDotK.length, interpolator.yDotK[0].length); + for (int k = 0; k < interpolator.yDotK.length; ++k) { + System.arraycopy(interpolator.yDotK[k], 0, yDotK[k], 0, interpolator.yDotK[k].length); + } + + } else { + previousState = null; + yDotK = null; + } + + // we cannot keep any reference to the equations in the copy + // the interpolator should have been finalized before + integrator = null; + + } + + /** Reinitialize the instance + * <p>Some Runge-Kutta integrators need fewer functions evaluations + * than their counterpart step interpolators. So the interpolator + * should perform the last evaluations they need by themselves. The + * {@link RungeKuttaFieldIntegrator RungeKuttaFieldIntegrator} and {@link + * EmbeddedRungeKuttaFieldIntegrator EmbeddedRungeKuttaFieldIntegrator} + * abstract classes call this method in order to let the step + * interpolator perform the evaluations it needs. These evaluations + * will be performed during the call to <code>doFinalize</code> if + * any, i.e. only if the step handler either calls the {@link + * AbstractFieldStepInterpolator#finalizeStep finalizeStep} method or the + * {@link AbstractFieldStepInterpolator#getInterpolatedState + * getInterpolatedState} method (for an interpolator which needs a + * finalization) or if it clones the step interpolator.</p> + * @param rkIntegrator integrator being used + * @param y reference to the integrator array holding the state at + * the end of the step + * @param yDotArray reference to the integrator array holding all the + * intermediate slopes + * @param forward integration direction indicator + * @param mapper equations mapper for the all equations + */ + public void reinitialize(final AbstractFieldIntegrator<T> rkIntegrator, + final T[] y, final T[][] yDotArray, final boolean forward, + final FieldEquationsMapper<T> mapper) { + reinitialize(y, forward, mapper); + this.previousState = null; + this.yDotK = yDotArray; + this.integrator = rkIntegrator; + } + + /** {@inheritDoc} */ + @Override + public void shift() { + previousState = currentState.clone(); + super.shift(); + } + +}
