diff options
author | Karl Shaffer <karlshaffer@google.com> | 2023-08-10 22:35:48 +0000 |
---|---|---|
committer | Automerger Merge Worker <android-build-automerger-merge-worker@system.gserviceaccount.com> | 2023-08-10 22:35:48 +0000 |
commit | 5484895ffd3d0c8337d159667cafc127c459f677 (patch) | |
tree | ace24ba4307d4978ee3134f7da671a77ad172da0 /src/main/java/org/apache/commons/math3/ode/nonstiff/RungeKuttaIntegrator.java | |
parent | bbf9548f049f99fd8e5a593baae983532dd983f4 (diff) | |
parent | b3715644fba79ef08acd9a2e157d078865281767 (diff) | |
download | apache-commons-math-5484895ffd3d0c8337d159667cafc127c459f677.tar.gz |
Check-in commons-math 3.6.1 am: 1354beaf45 am: 0018f64b87 am: b3715644fb
Original change: https://android-review.googlesource.com/c/platform/external/apache-commons-math/+/2702413
Change-Id: I5ad9b2a0822d668b5b6a62933c6d4c1f0b802001
Signed-off-by: Automerger Merge Worker <android-build-automerger-merge-worker@system.gserviceaccount.com>
Diffstat (limited to 'src/main/java/org/apache/commons/math3/ode/nonstiff/RungeKuttaIntegrator.java')
-rw-r--r-- | src/main/java/org/apache/commons/math3/ode/nonstiff/RungeKuttaIntegrator.java | 269 |
1 files changed, 269 insertions, 0 deletions
diff --git a/src/main/java/org/apache/commons/math3/ode/nonstiff/RungeKuttaIntegrator.java b/src/main/java/org/apache/commons/math3/ode/nonstiff/RungeKuttaIntegrator.java new file mode 100644 index 0000000..5f7d5d8 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/ode/nonstiff/RungeKuttaIntegrator.java @@ -0,0 +1,269 @@ +/* + * 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.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.AbstractIntegrator; +import org.apache.commons.math3.ode.ExpandableStatefulODE; +import org.apache.commons.math3.ode.FirstOrderDifferentialEquations; +import org.apache.commons.math3.util.FastMath; + +/** + * 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 EulerIntegrator + * @see ClassicalRungeKuttaIntegrator + * @see GillIntegrator + * @see MidpointIntegrator + * @since 1.2 + */ + +public abstract class RungeKuttaIntegrator extends AbstractIntegrator { + + /** 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 RungeKuttaStepInterpolator prototype; + + /** Integration step. */ + private final double step; + + /** Simple constructor. + * Build a Runge-Kutta integrator with the given + * step. The default step handler does nothing. + * @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 RungeKuttaIntegrator(final String name, + final double[] c, final double[][] a, final double[] b, + final RungeKuttaStepInterpolator prototype, + final double step) { + super(name); + this.c = c; + this.a = a; + this.b = b; + this.prototype = prototype; + this.step = FastMath.abs(step); + } + + /** {@inheritDoc} */ + @Override + public void integrate(final ExpandableStatefulODE equations, final double t) + throws NumberIsTooSmallException, DimensionMismatchException, + MaxCountExceededException, NoBracketingException { + + sanityChecks(equations, t); + setEquations(equations); + final boolean forward = t > equations.getTime(); + + // create some internal working arrays + final double[] y0 = equations.getCompleteState(); + final double[] y = y0.clone(); + final int stages = c.length + 1; + final double[][] yDotK = new double[stages][]; + for (int i = 0; i < stages; ++i) { + yDotK [i] = new double[y0.length]; + } + final double[] yTmp = y0.clone(); + final double[] yDotTmp = new double[y0.length]; + + // set up an interpolator sharing the integrator arrays + final RungeKuttaStepInterpolator interpolator = (RungeKuttaStepInterpolator) prototype.copy(); + interpolator.reinitialize(this, yTmp, yDotK, forward, + equations.getPrimaryMapper(), equations.getSecondaryMappers()); + interpolator.storeTime(equations.getTime()); + + // set up integration control objects + stepStart = equations.getTime(); + if (forward) { + if (stepStart + step >= t) { + stepSize = t - stepStart; + } else { + stepSize = step; + } + } else { + if (stepStart - step <= t) { + stepSize = t - stepStart; + } else { + stepSize = -step; + } + } + initIntegration(equations.getTime(), y0, t); + + // main integration loop + isLastStep = false; + do { + + interpolator.shift(); + + // first stage + computeDerivatives(stepStart, y, yDotK[0]); + + // next stages + for (int k = 1; k < stages; ++k) { + + for (int j = 0; j < y0.length; ++j) { + double sum = a[k-1][0] * yDotK[0][j]; + for (int l = 1; l < k; ++l) { + sum += a[k-1][l] * yDotK[l][j]; + } + yTmp[j] = y[j] + stepSize * sum; + } + + computeDerivatives(stepStart + c[k-1] * stepSize, yTmp, yDotK[k]); + + } + + // estimate the state at the end of the step + for (int j = 0; j < y0.length; ++j) { + double sum = b[0] * yDotK[0][j]; + for (int l = 1; l < stages; ++l) { + sum += b[l] * yDotK[l][j]; + } + yTmp[j] = y[j] + stepSize * sum; + } + + // discrete events handling + interpolator.storeTime(stepStart + stepSize); + System.arraycopy(yTmp, 0, y, 0, y0.length); + System.arraycopy(yDotK[stages - 1], 0, yDotTmp, 0, y0.length); + stepStart = acceptStep(interpolator, y, yDotTmp, t); + + if (!isLastStep) { + + // prepare next step + interpolator.storeTime(stepStart); + + // stepsize control for next step + final double nextT = stepStart + stepSize; + final boolean nextIsLast = forward ? (nextT >= t) : (nextT <= t); + if (nextIsLast) { + stepSize = t - stepStart; + } + } + + } while (!isLastStep); + + // dispatch results + equations.setTime(stepStart); + equations.setCompleteState(y); + + stepStart = Double.NaN; + stepSize = Double.NaN; + + } + + /** 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(ExpandableStatefulODE, double)} + * 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 double[] singleStep(final FirstOrderDifferentialEquations equations, + final double t0, final double[] y0, final double t) { + + // create some internal working arrays + final double[] y = y0.clone(); + final int stages = c.length + 1; + final double[][] yDotK = new double[stages][]; + for (int i = 0; i < stages; ++i) { + yDotK [i] = new double[y0.length]; + } + final double[] yTmp = y0.clone(); + + // first stage + final double h = t - t0; + equations.computeDerivatives(t0, y, yDotK[0]); + + // next stages + for (int k = 1; k < stages; ++k) { + + for (int j = 0; j < y0.length; ++j) { + double sum = a[k-1][0] * yDotK[0][j]; + for (int l = 1; l < k; ++l) { + sum += a[k-1][l] * yDotK[l][j]; + } + yTmp[j] = y[j] + h * sum; + } + + equations.computeDerivatives(t0 + c[k-1] * h, yTmp, yDotK[k]); + + } + + // estimate the state at the end of the step + for (int j = 0; j < y0.length; ++j) { + double sum = b[0] * yDotK[0][j]; + for (int l = 1; l < stages; ++l) { + sum += b[l] * yDotK[l][j]; + } + y[j] += h * sum; + } + + return y; + + } + +} |