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