View Javadoc
1   /*
2    * Licensed to the Apache Software Foundation (ASF) 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 ASF 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  /*
19   * This is not the original file distributed by the Apache Software Foundation
20   * It has been modified by the Hipparchus project
21   */
22  
23  package org.hipparchus.ode.nonstiff;
24  
25  import org.hipparchus.CalculusFieldElement;
26  import org.hipparchus.Field;
27  import org.hipparchus.ode.FieldEquationsMapper;
28  import org.hipparchus.ode.FieldODEStateAndDerivative;
29  import org.hipparchus.util.FastMath;
30  import org.hipparchus.util.MathArrays;
31  
32  
33  /**
34   * This class implements the 5(4) Dormand-Prince integrator for Ordinary
35   * Differential Equations.
36  
37   * <p>This integrator is an embedded Runge-Kutta integrator
38   * of order 5(4) used in local extrapolation mode (i.e. the solution
39   * is computed using the high order formula) with stepsize control
40   * (and automatic step initialization) and continuous output. This
41   * method uses 7 functions evaluations per step. However, since this
42   * is an <i>fsal</i>, the last evaluation of one step is the same as
43   * the first evaluation of the next step and hence can be avoided. So
44   * the cost is really 6 functions evaluations per step.</p>
45   *
46   * <p>This method has been published (whithout the continuous output
47   * that was added by Shampine in 1986) in the following article :</p>
48   * <pre>
49   *  A family of embedded Runge-Kutta formulae
50   *  J. R. Dormand and P. J. Prince
51   *  Journal of Computational and Applied Mathematics
52   *  volume 6, no 1, 1980, pp. 19-26
53   * </pre>
54   *
55   * @param <T> the type of the field elements
56   */
57  
58  public class DormandPrince54FieldIntegrator<T extends CalculusFieldElement<T>>
59      extends EmbeddedRungeKuttaFieldIntegrator<T> {
60  
61      /** Simple constructor.
62       * Build a fifth order Dormand-Prince integrator with the given step bounds
63       * @param field field to which the time and state vector elements belong
64       * @param minStep minimal step (sign is irrelevant, regardless of
65       * integration direction, forward or backward), the last step can
66       * be smaller than this
67       * @param maxStep maximal step (sign is irrelevant, regardless of
68       * integration direction, forward or backward), the last step can
69       * be smaller than this
70       * @param scalAbsoluteTolerance allowed absolute error
71       * @param scalRelativeTolerance allowed relative error
72       */
73      public DormandPrince54FieldIntegrator(final Field<T> field,
74                                            final double minStep, final double maxStep,
75                                            final double scalAbsoluteTolerance,
76                                            final double scalRelativeTolerance) {
77          super(field, DormandPrince54Integrator.METHOD_NAME, 6,
78                minStep, maxStep, scalAbsoluteTolerance, scalRelativeTolerance);
79      }
80  
81      /** Simple constructor.
82       * Build a fifth order Dormand-Prince integrator with the given step bounds
83       * @param field field to which the time and state vector elements belong
84       * @param minStep minimal step (sign is irrelevant, regardless of
85       * integration direction, forward or backward), the last step can
86       * be smaller than this
87       * @param maxStep maximal step (sign is irrelevant, regardless of
88       * integration direction, forward or backward), the last step can
89       * be smaller than this
90       * @param vecAbsoluteTolerance allowed absolute error
91       * @param vecRelativeTolerance allowed relative error
92       */
93      public DormandPrince54FieldIntegrator(final Field<T> field,
94                                            final double minStep, final double maxStep,
95                                            final double[] vecAbsoluteTolerance,
96                                            final double[] vecRelativeTolerance) {
97          super(field, DormandPrince54Integrator.METHOD_NAME, 6,
98                minStep, maxStep, vecAbsoluteTolerance, vecRelativeTolerance);
99      }
100 
101     /** {@inheritDoc} */
102     @Override
103     public T[] getC() {
104         final T[] c = MathArrays.buildArray(getField(), 6);
105         c[0] = fraction(1,  5);
106         c[1] = fraction(3, 10);
107         c[2] = fraction(4,  5);
108         c[3] = fraction(8,  9);
109         c[4] = getField().getOne();
110         c[5] = getField().getOne();
111         return c;
112     }
113 
114     /** {@inheritDoc} */
115     @Override
116     public T[][] getA() {
117         final T[][] a = MathArrays.buildArray(getField(), 6, -1);
118         for (int i = 0; i < a.length; ++i) {
119             a[i] = MathArrays.buildArray(getField(), i + 1);
120         }
121         a[0][0] = fraction(     1,     5);
122         a[1][0] = fraction(     3,    40);
123         a[1][1] = fraction(     9,    40);
124         a[2][0] = fraction(    44,    45);
125         a[2][1] = fraction(   -56,    15);
126         a[2][2] = fraction(    32,     9);
127         a[3][0] = fraction( 19372,  6561);
128         a[3][1] = fraction(-25360,  2187);
129         a[3][2] = fraction( 64448,  6561);
130         a[3][3] = fraction(  -212,   729);
131         a[4][0] = fraction(  9017,  3168);
132         a[4][1] = fraction(  -355,    33);
133         a[4][2] = fraction( 46732,  5247);
134         a[4][3] = fraction(    49,   176);
135         a[4][4] = fraction( -5103, 18656);
136         a[5][0] = fraction(    35,   384);
137         a[5][1] = getField().getZero();
138         a[5][2] = fraction(   500,  1113);
139         a[5][3] = fraction(   125,   192);
140         a[5][4] = fraction( -2187,  6784);
141         a[5][5] = fraction(    11,    84);
142         return a;
143     }
144 
145     /** {@inheritDoc} */
146     @Override
147     public T[] getB() {
148         final T[] b = MathArrays.buildArray(getField(), 7);
149         b[0] = fraction(   35,   384);
150         b[1] = getField().getZero();
151         b[2] = fraction(  500, 1113);
152         b[3] = fraction(  125,  192);
153         b[4] = fraction(-2187, 6784);
154         b[5] = fraction(   11,   84);
155         b[6] = getField().getZero();
156         return b;
157     }
158 
159     /** {@inheritDoc} */
160     @Override
161     protected DormandPrince54FieldStateInterpolator<T>
162         createInterpolator(final boolean forward, T[][] yDotK,
163                            final FieldODEStateAndDerivative<T> globalPreviousState,
164                            final FieldODEStateAndDerivative<T> globalCurrentState, final FieldEquationsMapper<T> mapper) {
165         return new DormandPrince54FieldStateInterpolator<T>(getField(), forward, yDotK,
166                                                            globalPreviousState, globalCurrentState,
167                                                            globalPreviousState, globalCurrentState,
168                                                            mapper);
169     }
170 
171     /** {@inheritDoc} */
172     @Override
173     public int getOrder() {
174         return 5;
175     }
176 
177     /** {@inheritDoc} */
178     @Override
179     protected double estimateError(final T[][] yDotK, final T[] y0, final T[] y1, final T h) {
180 
181         final StepsizeHelper helper = getStepSizeHelper();
182         double error = 0;
183 
184         for (int j = 0; j < helper.getMainSetDimension(); ++j) {
185             final double errSum = DormandPrince54Integrator.E1 * yDotK[0][j].getReal() +  DormandPrince54Integrator.E3 * yDotK[2][j].getReal() +
186                                   DormandPrince54Integrator.E4 * yDotK[3][j].getReal() +  DormandPrince54Integrator.E5 * yDotK[4][j].getReal() +
187                                   DormandPrince54Integrator.E6 * yDotK[5][j].getReal() +  DormandPrince54Integrator.E7 * yDotK[6][j].getReal();
188             final double tol = helper.getTolerance(j, FastMath.max(FastMath.abs(y0[j].getReal()), FastMath.abs(y1[j].getReal())));
189             final double ratio  = h.getReal() * errSum / tol;
190             error += ratio * ratio;
191         }
192 
193         return FastMath.sqrt(error / helper.getMainSetDimension());
194 
195     }
196 
197 }