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 020import org.apache.commons.math4.legacy.core.Field; 021import org.apache.commons.math4.legacy.core.RealFieldElement; 022import org.apache.commons.math4.legacy.ode.FieldEquationsMapper; 023import org.apache.commons.math4.legacy.ode.FieldODEStateAndDerivative; 024import org.apache.commons.math4.legacy.core.MathArrays; 025 026 027/** 028 * This class implements the 5(4) Dormand-Prince integrator for Ordinary 029 * Differential Equations. 030 031 * <p>This integrator is an embedded Runge-Kutta integrator 032 * of order 5(4) used in local extrapolation mode (i.e. the solution 033 * is computed using the high order formula) with stepsize control 034 * (and automatic step initialization) and continuous output. This 035 * method uses 7 functions evaluations per step. However, since this 036 * is an <i>fsal</i>, the last evaluation of one step is the same as 037 * the first evaluation of the next step and hence can be avoided. So 038 * the cost is really 6 functions evaluations per step.</p> 039 * 040 * <p>This method has been published (whithout the continuous output 041 * that was added by Shampine in 1986) in the following article : 042 * <pre> 043 * A family of embedded Runge-Kutta formulae 044 * J. R. Dormand and P. J. Prince 045 * Journal of Computational and Applied Mathematics 046 * volume 6, no 1, 1980, pp. 19-26 047 * </pre> 048 * 049 * @param <T> the type of the field elements 050 * @since 3.6 051 */ 052 053public class DormandPrince54FieldIntegrator<T extends RealFieldElement<T>> 054 extends EmbeddedRungeKuttaFieldIntegrator<T> { 055 056 /** Integrator method name. */ 057 private static final String METHOD_NAME = "Dormand-Prince 5(4)"; 058 059 /** Error array, element 1. */ 060 private final T e1; 061 062 // element 2 is zero, so it is neither stored nor used 063 064 /** Error array, element 3. */ 065 private final T e3; 066 067 /** Error array, element 4. */ 068 private final T e4; 069 070 /** Error array, element 5. */ 071 private final T e5; 072 073 /** Error array, element 6. */ 074 private final T e6; 075 076 /** Error array, element 7. */ 077 private final T e7; 078 079 /** Simple constructor. 080 * Build a fifth order Dormand-Prince integrator with the given step bounds 081 * @param field field to which the time and state vector elements belong 082 * @param minStep minimal step (sign is irrelevant, regardless of 083 * integration direction, forward or backward), the last step can 084 * be smaller than this 085 * @param maxStep maximal step (sign is irrelevant, regardless of 086 * integration direction, forward or backward), the last step can 087 * be smaller than this 088 * @param scalAbsoluteTolerance allowed absolute error 089 * @param scalRelativeTolerance allowed relative error 090 */ 091 public DormandPrince54FieldIntegrator(final Field<T> field, 092 final double minStep, final double maxStep, 093 final double scalAbsoluteTolerance, 094 final double scalRelativeTolerance) { 095 super(field, METHOD_NAME, 6, 096 minStep, maxStep, scalAbsoluteTolerance, scalRelativeTolerance); 097 e1 = fraction( 71, 57600); 098 e3 = fraction( -71, 16695); 099 e4 = fraction( 71, 1920); 100 e5 = fraction(-17253, 339200); 101 e6 = fraction( 22, 525); 102 e7 = fraction( -1, 40); 103 } 104 105 /** Simple constructor. 106 * Build a fifth order Dormand-Prince integrator with the given step bounds 107 * @param field field to which the time and state vector elements belong 108 * @param minStep minimal step (sign is irrelevant, regardless of 109 * integration direction, forward or backward), the last step can 110 * be smaller than this 111 * @param maxStep maximal step (sign is irrelevant, regardless of 112 * integration direction, forward or backward), the last step can 113 * be smaller than this 114 * @param vecAbsoluteTolerance allowed absolute error 115 * @param vecRelativeTolerance allowed relative error 116 */ 117 public DormandPrince54FieldIntegrator(final Field<T> field, 118 final double minStep, final double maxStep, 119 final double[] vecAbsoluteTolerance, 120 final double[] vecRelativeTolerance) { 121 super(field, METHOD_NAME, 6, 122 minStep, maxStep, vecAbsoluteTolerance, vecRelativeTolerance); 123 e1 = fraction( 71, 57600); 124 e3 = fraction( -71, 16695); 125 e4 = fraction( 71, 1920); 126 e5 = fraction(-17253, 339200); 127 e6 = fraction( 22, 525); 128 e7 = fraction( -1, 40); 129 } 130 131 /** {@inheritDoc} */ 132 @Override 133 public T[] getC() { 134 final T[] c = MathArrays.buildArray(getField(), 6); 135 c[0] = fraction(1, 5); 136 c[1] = fraction(3, 10); 137 c[2] = fraction(4, 5); 138 c[3] = fraction(8, 9); 139 c[4] = getField().getOne(); 140 c[5] = getField().getOne(); 141 return c; 142 } 143 144 /** {@inheritDoc} */ 145 @Override 146 public T[][] getA() { 147 final T[][] a = MathArrays.buildArray(getField(), 6, -1); 148 for (int i = 0; i < a.length; ++i) { 149 a[i] = MathArrays.buildArray(getField(), i + 1); 150 } 151 a[0][0] = fraction( 1, 5); 152 a[1][0] = fraction( 3, 40); 153 a[1][1] = fraction( 9, 40); 154 a[2][0] = fraction( 44, 45); 155 a[2][1] = fraction( -56, 15); 156 a[2][2] = fraction( 32, 9); 157 a[3][0] = fraction( 19372, 6561); 158 a[3][1] = fraction(-25360, 2187); 159 a[3][2] = fraction( 64448, 6561); 160 a[3][3] = fraction( -212, 729); 161 a[4][0] = fraction( 9017, 3168); 162 a[4][1] = fraction( -355, 33); 163 a[4][2] = fraction( 46732, 5247); 164 a[4][3] = fraction( 49, 176); 165 a[4][4] = fraction( -5103, 18656); 166 a[5][0] = fraction( 35, 384); 167 a[5][1] = getField().getZero(); 168 a[5][2] = fraction( 500, 1113); 169 a[5][3] = fraction( 125, 192); 170 a[5][4] = fraction( -2187, 6784); 171 a[5][5] = fraction( 11, 84); 172 return a; 173 } 174 175 /** {@inheritDoc} */ 176 @Override 177 public T[] getB() { 178 final T[] b = MathArrays.buildArray(getField(), 7); 179 b[0] = fraction( 35, 384); 180 b[1] = getField().getZero(); 181 b[2] = fraction( 500, 1113); 182 b[3] = fraction( 125, 192); 183 b[4] = fraction(-2187, 6784); 184 b[5] = fraction( 11, 84); 185 b[6] = getField().getZero(); 186 return b; 187 } 188 189 /** {@inheritDoc} */ 190 @Override 191 protected DormandPrince54FieldStepInterpolator<T> 192 createInterpolator(final boolean forward, T[][] yDotK, 193 final FieldODEStateAndDerivative<T> globalPreviousState, 194 final FieldODEStateAndDerivative<T> globalCurrentState, final FieldEquationsMapper<T> mapper) { 195 return new DormandPrince54FieldStepInterpolator<>(getField(), forward, yDotK, 196 globalPreviousState, globalCurrentState, 197 globalPreviousState, globalCurrentState, 198 mapper); 199 } 200 201 /** {@inheritDoc} */ 202 @Override 203 public int getOrder() { 204 return 5; 205 } 206 207 /** {@inheritDoc} */ 208 @Override 209 protected T estimateError(final T[][] yDotK, final T[] y0, final T[] y1, final T h) { 210 211 T error = getField().getZero(); 212 213 for (int j = 0; j < mainSetDimension; ++j) { 214 final T errSum = yDotK[0][j].multiply(e1). 215 add(yDotK[2][j].multiply(e3)). 216 add(yDotK[3][j].multiply(e4)). 217 add(yDotK[4][j].multiply(e5)). 218 add(yDotK[5][j].multiply(e6)). 219 add(yDotK[6][j].multiply(e7)); 220 221 final T yScale = RealFieldElement.max(y0[j].abs(), y1[j].abs()); 222 final T tol = (vecAbsoluteTolerance == null) ? 223 yScale.multiply(scalRelativeTolerance).add(scalAbsoluteTolerance) : 224 yScale.multiply(vecRelativeTolerance[j]).add(vecAbsoluteTolerance[j]); 225 final T ratio = h.multiply(errSum).divide(tol); 226 error = error.add(ratio.multiply(ratio)); 227 } 228 229 return error.divide(mainSetDimension).sqrt(); 230 } 231}