1 /*
2 * Licensed to the Hipparchus project under one or more
3 * contributor license agreements. See the NOTICE file distributed with
4 * this work for additional information regarding copyright ownership.
5 * The Hipparchus project licenses this file to You under the Apache License, Version 2.0
6 * (the "License"); you may not use this file except in compliance with
7 * the License. You may obtain a copy of the License at
8 *
9 * https://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 */
17
18 package org.hipparchus.ode.nonstiff;
19
20 import org.hipparchus.ode.EquationsMapper;
21 import org.hipparchus.ode.ODEStateAndDerivative;
22 import org.hipparchus.ode.nonstiff.interpolators.DormandPrince54StateInterpolator;
23 import org.hipparchus.util.FastMath;
24
25
26 /**
27 * This class implements the 5(4) Dormand-Prince integrator for Ordinary
28 * Differential Equations.
29
30 * <p>This integrator is an embedded Runge-Kutta integrator
31 * of order 5(4) used in local extrapolation mode (i.e. the solution
32 * is computed using the high order formula) with stepsize control
33 * (and automatic step initialization) and continuous output. This
34 * method uses 7 functions evaluations per step. However, since this
35 * is an <i>fsal</i>, the last evaluation of one step is the same as
36 * the first evaluation of the next step and hence can be avoided. So
37 * the cost is really 6 functions evaluations per step.</p>
38 *
39 * <p>This method has been published (whithout the continuous output
40 * that was added by Shampine in 1986) in the following article :</p>
41 * <pre>
42 * A family of embedded Runge-Kutta formulae
43 * J. R. Dormand and P. J. Prince
44 * Journal of Computational and Applied Mathematics
45 * volume 6, no 1, 1980, pp. 19-26
46 * </pre>
47 *
48 */
49
50 public class DormandPrince54Integrator extends EmbeddedRungeKuttaIntegrator {
51
52 /** Name of integration scheme. */
53 public static final String METHOD_NAME = "Dormand-Prince 5 (4)";
54
55 /** Error array, element 1. */
56 static final double E1 = 71.0 / 57600.0;
57
58 // element 2 is zero, so it is neither stored nor used
59
60 /** Error array, element 3. */
61 static final double E3 = -71.0 / 16695.0;
62
63 /** Error array, element 4. */
64 static final double E4 = 71.0 / 1920.0;
65
66 /** Error array, element 5. */
67 static final double E5 = -17253.0 / 339200.0;
68
69 /** Error array, element 6. */
70 static final double E6 = 22.0 / 525.0;
71
72 /** Error array, element 7. */
73 static final double E7 = -1.0 / 40.0;
74
75 /** Simple constructor.
76 * Build a fifth order Dormand-Prince integrator with the given step bounds
77 * @param minStep minimal step (sign is irrelevant, regardless of
78 * integration direction, forward or backward), the last step can
79 * be smaller than this
80 * @param maxStep maximal step (sign is irrelevant, regardless of
81 * integration direction, forward or backward), the last step can
82 * be smaller than this
83 * @param scalAbsoluteTolerance allowed absolute error
84 * @param scalRelativeTolerance allowed relative error
85 */
86 public DormandPrince54Integrator(final double minStep, final double maxStep,
87 final double scalAbsoluteTolerance,
88 final double scalRelativeTolerance) {
89 super(METHOD_NAME, 6,
90 minStep, maxStep, scalAbsoluteTolerance, scalRelativeTolerance);
91 }
92
93 /** Simple constructor.
94 * Build a fifth order Dormand-Prince integrator with the given step bounds
95 * @param minStep minimal step (sign is irrelevant, regardless of
96 * integration direction, forward or backward), the last step can
97 * be smaller than this
98 * @param maxStep maximal step (sign is irrelevant, regardless of
99 * integration direction, forward or backward), the last step can
100 * be smaller than this
101 * @param vecAbsoluteTolerance allowed absolute error
102 * @param vecRelativeTolerance allowed relative error
103 */
104 public DormandPrince54Integrator(final double minStep, final double maxStep,
105 final double[] vecAbsoluteTolerance,
106 final double[] vecRelativeTolerance) {
107 super(METHOD_NAME, 6,
108 minStep, maxStep, vecAbsoluteTolerance, vecRelativeTolerance);
109 }
110
111 /** {@inheritDoc} */
112 @Override
113 public double[] getC() {
114 return new double[] {
115 1.0 / 5.0, 3.0 / 10.0, 4.0 / 5.0, 8.0 / 9.0, 1.0, 1.0
116 };
117 }
118
119 /** {@inheritDoc} */
120 @Override
121 public double[][] getA() {
122 return new double[][] {
123 { 1.0 / 5.0 },
124 { 3.0 / 40.0, 9.0 / 40.0 },
125 { 44.0 / 45.0, -56.0 / 15.0, 32.0 / 9.0 },
126 { 19372.0 / 6561.0, -25360.0 / 2187.0, 64448.0 / 6561.0, -212.0 / 729.0 },
127 { 9017.0 / 3168.0, -355.0 / 33.0, 46732.0 / 5247.0, 49.0 / 176.0, -5103.0 / 18656.0 },
128 { 35.0 / 384.0, 0.0, 500.0 / 1113.0, 125.0 / 192.0, -2187.0 / 6784.0, 11.0 / 84.0 }
129 };
130 }
131
132 /** {@inheritDoc} */
133 @Override
134 public double[] getB() {
135 return new double[] {
136 35.0 / 384.0, 0.0, 500.0 / 1113.0, 125.0 / 192.0, -2187.0 / 6784.0, 11.0 / 84.0, 0.0
137 };
138 }
139
140 /** {@inheritDoc} */
141 @Override
142 protected DormandPrince54StateInterpolator
143 createInterpolator(final boolean forward, double[][] yDotK,
144 final ODEStateAndDerivative globalPreviousState,
145 final ODEStateAndDerivative globalCurrentState,
146 final EquationsMapper mapper) {
147 return new DormandPrince54StateInterpolator(forward, yDotK,
148 globalPreviousState, globalCurrentState,
149 globalPreviousState, globalCurrentState,
150 mapper);
151 }
152
153 /** {@inheritDoc} */
154 @Override
155 public int getOrder() {
156 return 5;
157 }
158
159 /** {@inheritDoc} */
160 @Override
161 protected double estimateError(final double[][] yDotK,
162 final double[] y0, final double[] y1,
163 final double h) {
164
165 final StepsizeHelper helper = getStepSizeHelper();
166 double error = 0;
167
168 for (int j = 0; j < helper.getMainSetDimension(); ++j) {
169 final double errSum = E1 * yDotK[0][j] + E3 * yDotK[2][j] +
170 E4 * yDotK[3][j] + E5 * yDotK[4][j] +
171 E6 * yDotK[5][j] + E7 * yDotK[6][j];
172
173 final double tol = helper.getTolerance(j, FastMath.max(FastMath.abs(y0[j]), FastMath.abs(y1[j])));
174 final double ratio = h * errSum / tol;
175 error += ratio * ratio;
176
177 }
178
179 return FastMath.sqrt(error / helper.getMainSetDimension());
180
181 }
182
183 }