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.math4.legacy.ode.nonstiff; 019 020 021import org.apache.commons.math4.legacy.core.Field; 022import org.apache.commons.math4.legacy.core.RealFieldElement; 023import org.apache.commons.math4.legacy.exception.DimensionMismatchException; 024import org.apache.commons.math4.legacy.exception.MaxCountExceededException; 025import org.apache.commons.math4.legacy.exception.NoBracketingException; 026import org.apache.commons.math4.legacy.exception.NumberIsTooSmallException; 027import org.apache.commons.math4.legacy.ode.AbstractFieldIntegrator; 028import org.apache.commons.math4.legacy.ode.FieldEquationsMapper; 029import org.apache.commons.math4.legacy.ode.FieldExpandableODE; 030import org.apache.commons.math4.legacy.ode.FirstOrderFieldDifferentialEquations; 031import org.apache.commons.math4.legacy.ode.FieldODEState; 032import org.apache.commons.math4.legacy.ode.FieldODEStateAndDerivative; 033import org.apache.commons.math4.legacy.core.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 * 051 * @see EulerFieldIntegrator 052 * @see ClassicalRungeKuttaFieldIntegrator 053 * @see GillFieldIntegrator 054 * @see MidpointFieldIntegrator 055 * @param <T> the type of the field elements 056 * @since 3.6 057 */ 058 059public abstract class RungeKuttaFieldIntegrator<T extends RealFieldElement<T>> 060 extends AbstractFieldIntegrator<T> 061 implements FieldButcherArrayProvider<T> { 062 063 /** Time steps from Butcher array (without the first zero). */ 064 private final T[] c; 065 066 /** Internal weights from Butcher array (without the first empty row). */ 067 private final T[][] a; 068 069 /** External weights for the high order method from Butcher array. */ 070 private final T[] b; 071 072 /** Integration step. */ 073 private final T step; 074 075 /** Simple constructor. 076 * Build a Runge-Kutta integrator with the given 077 * step. The default step handler does nothing. 078 * @param field field to which the time and state vector elements belong 079 * @param name name of the method 080 * @param step integration step 081 */ 082 protected RungeKuttaFieldIntegrator(final Field<T> field, final String name, final T step) { 083 super(field, name); 084 this.c = getC(); 085 this.a = getA(); 086 this.b = getB(); 087 this.step = step.abs(); 088 } 089 090 /** Create a fraction. 091 * @param p numerator 092 * @param q denominator 093 * @return p/q computed in the instance field 094 */ 095 protected T fraction(final int p, final int q) { 096 return getField().getZero().add(p).divide(q); 097 } 098 099 /** Create an interpolator. 100 * @param forward integration direction indicator 101 * @param yDotK slopes at the intermediate points 102 * @param globalPreviousState start of the global step 103 * @param globalCurrentState end of the global step 104 * @param mapper equations mapper for the all equations 105 * @return external weights for the high order method from Butcher array 106 */ 107 protected abstract RungeKuttaFieldStepInterpolator<T> createInterpolator(boolean forward, T[][] yDotK, 108 FieldODEStateAndDerivative<T> globalPreviousState, 109 FieldODEStateAndDerivative<T> globalCurrentState, 110 FieldEquationsMapper<T> mapper); 111 112 /** {@inheritDoc} */ 113 @Override 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 // estimate the state at the end of the step 169 for (int j = 0; j < y0.length; ++j) { 170 T sum = yDotK[0][j].multiply(b[0]); 171 for (int l = 1; l < stages; ++l) { 172 sum = sum.add(yDotK[l][j].multiply(b[l])); 173 } 174 yTmp[j] = y[j].add(getStepSize().multiply(sum)); 175 } 176 final T stepEnd = getStepStart().getTime().add(getStepSize()); 177 final T[] yDotTmp = computeDerivatives(stepEnd, yTmp); 178 final FieldODEStateAndDerivative<T> stateTmp = new FieldODEStateAndDerivative<>(stepEnd, yTmp, yDotTmp); 179 180 // discrete events handling 181 System.arraycopy(yTmp, 0, y, 0, y0.length); 182 setStepStart(acceptStep(createInterpolator(forward, yDotK, getStepStart(), stateTmp, equations.getMapper()), 183 finalTime)); 184 185 if (!isLastStep()) { 186 187 // stepsize control for next step 188 final T nextT = getStepStart().getTime().add(getStepSize()); 189 final boolean nextIsLast = forward ? 190 (nextT.subtract(finalTime).getReal() >= 0) : 191 (nextT.subtract(finalTime).getReal() <= 0); 192 if (nextIsLast) { 193 setStepSize(finalTime.subtract(getStepStart().getTime())); 194 } 195 } 196 } while (!isLastStep()); 197 198 final FieldODEStateAndDerivative<T> finalState = getStepStart(); 199 setStepStart(null); 200 setStepSize(null); 201 return finalState; 202 } 203 204 /** Fast computation of a single step of ODE integration. 205 * <p>This method is intended for the limited use case of 206 * very fast computation of only one step without using any of the 207 * rich features of general integrators that may take some time 208 * to set up (i.e. no step handlers, no events handlers, no additional 209 * states, no interpolators, no error control, no evaluations count, 210 * no sanity checks ...). It handles the strict minimum of computation, 211 * so it can be embedded in outer loops.</p> 212 * <p> 213 * This method is <em>not</em> used at all by the {@link #integrate(FieldExpandableODE, 214 * FieldODEState, RealFieldElement)} method. It also completely ignores the step set at 215 * construction time, and uses only a single step to go from {@code t0} to {@code t}. 216 * </p> 217 * <p> 218 * As this method does not use any of the state-dependent features of the integrator, 219 * it should be reasonably thread-safe <em>if and only if</em> the provided differential 220 * equations are themselves thread-safe. 221 * </p> 222 * @param equations differential equations to integrate 223 * @param t0 initial time 224 * @param y0 initial value of the state vector at t0 225 * @param t target time for the integration 226 * (can be set to a value smaller than {@code t0} for backward integration) 227 * @return state vector at {@code t} 228 */ 229 public T[] singleStep(final FirstOrderFieldDifferentialEquations<T> equations, 230 final T t0, final T[] y0, final T t) { 231 232 // create some internal working arrays 233 final T[] y = y0.clone(); 234 final int stages = c.length + 1; 235 final T[][] yDotK = MathArrays.buildArray(getField(), stages, -1); 236 final T[] yTmp = y0.clone(); 237 238 // first stage 239 final T h = t.subtract(t0); 240 yDotK[0] = equations.computeDerivatives(t0, y); 241 242 // next stages 243 for (int k = 1; k < stages; ++k) { 244 245 for (int j = 0; j < y0.length; ++j) { 246 T sum = yDotK[0][j].multiply(a[k-1][0]); 247 for (int l = 1; l < k; ++l) { 248 sum = sum.add(yDotK[l][j].multiply(a[k-1][l])); 249 } 250 yTmp[j] = y[j].add(h.multiply(sum)); 251 } 252 253 yDotK[k] = equations.computeDerivatives(t0.add(h.multiply(c[k-1])), yTmp); 254 } 255 256 // estimate the state at the end of the step 257 for (int j = 0; j < y0.length; ++j) { 258 T sum = yDotK[0][j].multiply(b[0]); 259 for (int l = 1; l < stages; ++l) { 260 sum = sum.add(yDotK[l][j].multiply(b[l])); 261 } 262 y[j] = y[j].add(h.multiply(sum)); 263 } 264 265 return y; 266 } 267}