|
Commons Math example source code file (EmbeddedRungeKuttaIntegrator.java)
The Commons Math EmbeddedRungeKuttaIntegrator.java source code/* * 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.math.ode.nonstiff; import org.apache.commons.math.ode.DerivativeException; import org.apache.commons.math.ode.FirstOrderDifferentialEquations; import org.apache.commons.math.ode.IntegratorException; import org.apache.commons.math.ode.events.CombinedEventsManager; import org.apache.commons.math.ode.sampling.AbstractStepInterpolator; import org.apache.commons.math.ode.sampling.DummyStepInterpolator; import org.apache.commons.math.ode.sampling.StepHandler; /** * This class implements the common part of all embedded Runge-Kutta * integrators for Ordinary Differential Equations. * * <p>These methods are embedded explicit Runge-Kutta methods with two * sets of coefficients allowing to estimate the error, their Butcher * arrays are as follows : * <pre> * 0 | * c2 | a21 * c3 | a31 a32 * ... | ... * cs | as1 as2 ... ass-1 * |-------------------------- * | b1 b2 ... bs-1 bs * | b'1 b'2 ... b's-1 b's * </pre> * </p> * * <p>In fact, we rather use the array defined by ej = bj - b'j to * compute directly the error rather than computing two estimates and * then comparing them.</p> * * <p>Some methods are qualified as fsal (first same as last) * methods. This means the last evaluation of the derivatives in one * step is the same as the first in the next step. Then, this * evaluation can be reused from one step to the next one and the cost * of such a method is really s-1 evaluations despite the method still * has s stages. This behaviour is true only for successful steps, if * the step is rejected after the error estimation phase, no * evaluation is saved. For an <i>fsal method, we have cs = 1 and * asi = bi for all i.</p> * * @version $Revision: 927202 $ $Date: 2010-03-24 18:11:51 -0400 (Wed, 24 Mar 2010) $ * @since 1.2 */ public abstract class EmbeddedRungeKuttaIntegrator extends AdaptiveStepsizeIntegrator { /** Indicator for <i>fsal methods. */ private final boolean fsal; /** 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; /** Stepsize control exponent. */ private final double exp; /** Safety factor for stepsize control. */ private double safety; /** Minimal reduction factor for stepsize control. */ private double minReduction; /** Maximal growth factor for stepsize control. */ private double maxGrowth; /** Build a Runge-Kutta integrator with the given Butcher array. * @param name name of the method * @param fsal indicate that the method is an <i>fsal * @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 minStep minimal step (must be positive even for backward * integration), the last step can be smaller than this * @param maxStep maximal step (must be positive even for backward * integration) * @param scalAbsoluteTolerance allowed absolute error * @param scalRelativeTolerance allowed relative error */ protected EmbeddedRungeKuttaIntegrator(final String name, final boolean fsal, final double[] c, final double[][] a, final double[] b, final RungeKuttaStepInterpolator prototype, final double minStep, final double maxStep, final double scalAbsoluteTolerance, final double scalRelativeTolerance) { super(name, minStep, maxStep, scalAbsoluteTolerance, scalRelativeTolerance); this.fsal = fsal; this.c = c; this.a = a; this.b = b; this.prototype = prototype; exp = -1.0 / getOrder(); // set the default values of the algorithm control parameters setSafety(0.9); setMinReduction(0.2); setMaxGrowth(10.0); } /** Build a Runge-Kutta integrator with the given Butcher array. * @param name name of the method * @param fsal indicate that the method is an <i>fsal * @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 minStep minimal step (must be positive even for backward * integration), the last step can be smaller than this * @param maxStep maximal step (must be positive even for backward * integration) * @param vecAbsoluteTolerance allowed absolute error * @param vecRelativeTolerance allowed relative error */ protected EmbeddedRungeKuttaIntegrator(final String name, final boolean fsal, final double[] c, final double[][] a, final double[] b, final RungeKuttaStepInterpolator prototype, final double minStep, final double maxStep, final double[] vecAbsoluteTolerance, final double[] vecRelativeTolerance) { super(name, minStep, maxStep, vecAbsoluteTolerance, vecRelativeTolerance); this.fsal = fsal; this.c = c; this.a = a; this.b = b; this.prototype = prototype; exp = -1.0 / getOrder(); // set the default values of the algorithm control parameters setSafety(0.9); setMinReduction(0.2); setMaxGrowth(10.0); } /** Get the order of the method. * @return order of the method */ public abstract int getOrder(); /** Get the safety factor for stepsize control. * @return safety factor */ public double getSafety() { return safety; } /** Set the safety factor for stepsize control. * @param safety safety factor */ public void setSafety(final double safety) { this.safety = safety; } /** {@inheritDoc} */ @Override public double integrate(final FirstOrderDifferentialEquations equations, final double t0, final double[] y0, final double t, final double[] y) throws DerivativeException, IntegratorException { sanityChecks(equations, t0, y0, t, y); setEquations(equations); resetEvaluations(); final boolean forward = t > t0; // create some internal working arrays final int stages = c.length + 1; if (y != y0) { System.arraycopy(y0, 0, y, 0, y0.length); } final double[][] yDotK = new double[stages][y0.length]; final double[] yTmp = new double[y0.length]; // set up an interpolator sharing the integrator arrays AbstractStepInterpolator interpolator; if (requiresDenseOutput() || (! eventsHandlersManager.isEmpty())) { final RungeKuttaStepInterpolator rki = (RungeKuttaStepInterpolator) prototype.copy(); rki.reinitialize(this, yTmp, yDotK, forward); interpolator = rki; } else { interpolator = new DummyStepInterpolator(yTmp, yDotK[stages - 1], forward); } interpolator.storeTime(t0); // set up integration control objects stepStart = t0; double hNew = 0; boolean firstTime = true; for (StepHandler handler : stepHandlers) { handler.reset(); } CombinedEventsManager manager = addEndTimeChecker(t0, t, eventsHandlersManager); boolean lastStep = false; // main integration loop while (!lastStep) { interpolator.shift(); double error = 0; for (boolean loop = true; loop;) { if (firstTime || !fsal) { // first stage computeDerivatives(stepStart, y, yDotK[0]); } if (firstTime) { final double[] scale = new double[y0.length]; if (vecAbsoluteTolerance == null) { for (int i = 0; i < scale.length; ++i) { scale[i] = scalAbsoluteTolerance + scalRelativeTolerance * Math.abs(y[i]); } } else { for (int i = 0; i < scale.length; ++i) { scale[i] = vecAbsoluteTolerance[i] + vecRelativeTolerance[i] * Math.abs(y[i]); } } hNew = initializeStep(equations, forward, getOrder(), scale, stepStart, y, yDotK[0], yTmp, yDotK[1]); firstTime = false; } stepSize = hNew; // 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; } // estimate the error at the end of the step error = estimateError(yDotK, y, yTmp, stepSize); if (error <= 1.0) { // discrete events handling interpolator.storeTime(stepStart + stepSize); if (manager.evaluateStep(interpolator)) { final double dt = manager.getEventTime() - stepStart; if (Math.abs(dt) <= Math.ulp(stepStart)) { // we cannot simply truncate the step, reject the current computation // and let the loop compute another state with the truncated step. // it is so small (much probably exactly 0 due to limited accuracy) // that the code above would fail handling it. // So we set up an artificial 0 size step by copying states interpolator.storeTime(stepStart); System.arraycopy(y, 0, yTmp, 0, y0.length); hNew = 0; stepSize = 0; loop = false; } else { // reject the step to match exactly the next switch time hNew = dt; } } else { // accept the step loop = false; } } else { // reject the step and attempt to reduce error by stepsize control final double factor = Math.min(maxGrowth, Math.max(minReduction, safety * Math.pow(error, exp))); hNew = filterStep(stepSize * factor, forward, false); } } // the step has been accepted final double nextStep = stepStart + stepSize; System.arraycopy(yTmp, 0, y, 0, y0.length); manager.stepAccepted(nextStep, y); lastStep = manager.stop(); // provide the step data to the step handler interpolator.storeTime(nextStep); for (StepHandler handler : stepHandlers) { handler.handleStep(interpolator, lastStep); } stepStart = nextStep; if (fsal) { // save the last evaluation for the next step System.arraycopy(yDotK[stages - 1], 0, yDotK[0], 0, y0.length); } if (manager.reset(stepStart, y) && ! lastStep) { // some event handler has triggered changes that // invalidate the derivatives, we need to recompute them computeDerivatives(stepStart, y, yDotK[0]); } if (! lastStep) { // in some rare cases we may get here with stepSize = 0, for example // when an event occurs at integration start, reducing the first step // to zero; we have to reset the step to some safe non zero value stepSize = filterStep(stepSize, forward, true); // stepsize control for next step final double factor = Math.min(maxGrowth, Math.max(minReduction, safety * Math.pow(error, exp))); final double scaledH = stepSize * factor; final double nextT = stepStart + scaledH; final boolean nextIsLast = forward ? (nextT >= t) : (nextT <= t); hNew = filterStep(scaledH, forward, nextIsLast); } } final double stopTime = stepStart; resetInternalState(); return stopTime; } /** Get the minimal reduction factor for stepsize control. * @return minimal reduction factor */ public double getMinReduction() { return minReduction; } /** Set the minimal reduction factor for stepsize control. * @param minReduction minimal reduction factor */ public void setMinReduction(final double minReduction) { this.minReduction = minReduction; } /** Get the maximal growth factor for stepsize control. * @return maximal growth factor */ public double getMaxGrowth() { return maxGrowth; } /** Set the maximal growth factor for stepsize control. * @param maxGrowth maximal growth factor */ public void setMaxGrowth(final double maxGrowth) { this.maxGrowth = maxGrowth; } /** Compute the error ratio. * @param yDotK derivatives computed during the first stages * @param y0 estimate of the step at the start of the step * @param y1 estimate of the step at the end of the step * @param h current step * @return error ratio, greater than 1 if step should be rejected */ protected abstract double estimateError(double[][] yDotK, double[] y0, double[] y1, double h); } Other Commons Math examples (source code examples)Here is a short list of links related to this Commons Math EmbeddedRungeKuttaIntegrator.java source code file: |
... this post is sponsored by my books ... | |
#1 New Release! |
FP Best Seller |
Copyright 1998-2024 Alvin Alexander, alvinalexander.com
All Rights Reserved.
A percentage of advertising revenue from
pages under the /java/jwarehouse
URI on this website is
paid back to open source projects.