All Downloads are FREE. Search and download functionalities are using the official Maven repository.

org.apache.commons.math3.ode.nonstiff.RungeKuttaIntegrator Maven / Gradle / Ivy

There is a newer version: 2.12.15
Show newest version
/*
 * 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.
 *
 * 

These methods are explicit Runge-Kutta methods, their Butcher * arrays are as follows : *

 *    0  |
 *   c2  | a21
 *   c3  | a31  a32
 *   ... |        ...
 *   cs  | as1  as2  ...  ass-1
 *       |--------------------------
 *       |  b1   b2  ...   bs-1  bs
 * 
*

* * @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. *

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.

*

* This method is not 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}. *

*

* As this method does not use any of the state-dependent features of the integrator, * it should be reasonably thread-safe if and only if the provided differential * equations are themselves thread-safe. *

* @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; } }




© 2015 - 2024 Weber Informatics LLC | Privacy Policy