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.math3.ode.nonstiff; 019 020 021import org.apache.commons.math3.Field; 022import org.apache.commons.math3.RealFieldElement; 023import org.apache.commons.math3.exception.DimensionMismatchException; 024import org.apache.commons.math3.exception.MaxCountExceededException; 025import org.apache.commons.math3.exception.NoBracketingException; 026import org.apache.commons.math3.exception.NumberIsTooSmallException; 027import org.apache.commons.math3.ode.AbstractFieldIntegrator; 028import org.apache.commons.math3.ode.FieldEquationsMapper; 029import org.apache.commons.math3.ode.FieldExpandableODE; 030import org.apache.commons.math3.ode.FirstOrderFieldDifferentialEquations; 031import org.apache.commons.math3.ode.FieldODEState; 032import org.apache.commons.math3.ode.FieldODEStateAndDerivative; 033import org.apache.commons.math3.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 * </p> 051 * 052 * @see EulerFieldIntegrator 053 * @see ClassicalRungeKuttaFieldIntegrator 054 * @see GillFieldIntegrator 055 * @see MidpointFieldIntegrator 056 * @param <T> the type of the field elements 057 * @since 3.6 058 */ 059 060public abstract class RungeKuttaFieldIntegrator<T extends RealFieldElement<T>> 061 extends AbstractFieldIntegrator<T> 062 implements FieldButcherArrayProvider<T> { 063 064 /** Time steps from Butcher array (without the first zero). */ 065 private final T[] c; 066 067 /** Internal weights from Butcher array (without the first empty row). */ 068 private final T[][] a; 069 070 /** External weights for the high order method from Butcher array. */ 071 private final T[] b; 072 073 /** Integration step. */ 074 private final T step; 075 076 /** Simple constructor. 077 * Build a Runge-Kutta integrator with the given 078 * step. The default step handler does nothing. 079 * @param field field to which the time and state vector elements belong 080 * @param name name of the method 081 * @param step integration step 082 */ 083 protected RungeKuttaFieldIntegrator(final Field<T> field, final String name, final T step) { 084 super(field, name); 085 this.c = getC(); 086 this.a = getA(); 087 this.b = getB(); 088 this.step = step.abs(); 089 } 090 091 /** Create a fraction. 092 * @param p numerator 093 * @param q denominator 094 * @return p/q computed in the instance field 095 */ 096 protected T fraction(final int p, final int q) { 097 return getField().getZero().add(p).divide(q); 098 } 099 100 /** Create an interpolator. 101 * @param forward integration direction indicator 102 * @param yDotK slopes at the intermediate points 103 * @param globalPreviousState start of the global step 104 * @param globalCurrentState end of the global step 105 * @param mapper equations mapper for the all equations 106 * @return external weights for the high order method from Butcher array 107 */ 108 protected abstract RungeKuttaFieldStepInterpolator<T> createInterpolator(boolean forward, T[][] yDotK, 109 final FieldODEStateAndDerivative<T> globalPreviousState, 110 final FieldODEStateAndDerivative<T> globalCurrentState, 111 FieldEquationsMapper<T> mapper); 112 113 /** {@inheritDoc} */ 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<T>(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}