View Javadoc
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) Higham and Hall integrator for
27   * Ordinary 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.</p>
34   *
35   */
36  
37  public class HighamHall54Integrator extends EmbeddedRungeKuttaIntegrator {
38  
39      /** Integrator method name. */
40      static final String METHOD_NAME = "Higham-Hall 5(4)";
41  
42      /** Error weights Butcher array. */
43      static final double[] STATIC_E = {
44          -1.0/20.0, 0.0, 81.0/160.0, -6.0/5.0, 25.0/32.0, 1.0/16.0, -1.0/10.0
45      };
46  
47      /** Simple constructor.
48       * Build a fifth order Higham and Hall integrator with the given step bounds
49       * @param minStep minimal step (sign is irrelevant, regardless of
50       * integration direction, forward or backward), the last step can
51       * be smaller than this
52       * @param maxStep maximal step (sign is irrelevant, regardless of
53       * integration direction, forward or backward), the last step can
54       * be smaller than this
55       * @param scalAbsoluteTolerance allowed absolute error
56       * @param scalRelativeTolerance allowed relative error
57       */
58      public HighamHall54Integrator(final double minStep, final double maxStep,
59                                    final double scalAbsoluteTolerance,
60                                    final double scalRelativeTolerance) {
61          super(METHOD_NAME, -1,
62                minStep, maxStep, scalAbsoluteTolerance, scalRelativeTolerance);
63      }
64  
65      /** Simple constructor.
66       * Build a fifth order Higham and Hall integrator with the given step bounds
67       * @param minStep minimal step (sign is irrelevant, regardless of
68       * integration direction, forward or backward), the last step can
69       * be smaller than this
70       * @param maxStep maximal step (sign is irrelevant, regardless of
71       * integration direction, forward or backward), the last step can
72       * be smaller than this
73       * @param vecAbsoluteTolerance allowed absolute error
74       * @param vecRelativeTolerance allowed relative error
75       */
76      public HighamHall54Integrator(final double minStep, final double maxStep,
77                                    final double[] vecAbsoluteTolerance,
78                                    final double[] vecRelativeTolerance) {
79          super(METHOD_NAME, -1,
80                minStep, maxStep, vecAbsoluteTolerance, vecRelativeTolerance);
81      }
82  
83      /** {@inheritDoc} */
84      @Override
85      public double[] getC() {
86          return new double[] {
87              2.0/9.0, 1.0/3.0, 1.0/2.0, 3.0/5.0, 1.0, 1.0
88          };
89      }
90  
91      /** {@inheritDoc} */
92      @Override
93      public double[][] getA() {
94          return new double[][] {
95              {2.0/9.0},
96              {1.0/12.0, 1.0/4.0},
97              {1.0/8.0, 0.0, 3.0/8.0},
98              {91.0/500.0, -27.0/100.0, 78.0/125.0, 8.0/125.0},
99              {-11.0/20.0, 27.0/20.0, 12.0/5.0, -36.0/5.0, 5.0},
100             {1.0/12.0, 0.0, 27.0/32.0, -4.0/3.0, 125.0/96.0, 5.0/48.0}
101         };
102     }
103 
104     /** {@inheritDoc} */
105     @Override
106     public double[] getB() {
107         return new double[] {
108             1.0/12.0, 0.0, 27.0/32.0, -4.0/3.0, 125.0/96.0, 5.0/48.0, 0.0
109         };
110     }
111 
112     /** {@inheritDoc} */
113     @Override
114     protected HighamHall54StateInterpolator
115     createInterpolator(final boolean forward, double[][] yDotK,
116                        final ODEStateAndDerivative globalPreviousState,
117                        final ODEStateAndDerivative globalCurrentState,
118                        final EquationsMapper mapper) {
119         return new HighamHall54StateInterpolator(forward, yDotK,
120                                                 globalPreviousState, globalCurrentState,
121                                                 globalPreviousState, globalCurrentState,
122                                                 mapper);
123     }
124 
125     /** {@inheritDoc} */
126     @Override
127     public int getOrder() {
128         return 5;
129     }
130 
131     /** {@inheritDoc} */
132     @Override
133     protected double estimateError(final double[][] yDotK,
134                                    final double[] y0, final double[] y1,
135                                    final double h) {
136 
137         final StepsizeHelper helper = getStepSizeHelper();
138         double error = 0;
139 
140         for (int j = 0; j < helper.getMainSetDimension(); ++j) {
141             double errSum = STATIC_E[0] * yDotK[0][j];
142             for (int l = 1; l < STATIC_E.length; ++l) {
143                 errSum += STATIC_E[l] * yDotK[l][j];
144             }
145 
146             final double tol = helper.getTolerance(j, FastMath.max(FastMath.abs(y0[j]), FastMath.abs(y1[j])));
147             final double ratio  = h * errSum / tol;
148             error += ratio * ratio;
149 
150         }
151 
152         return FastMath.sqrt(error / helper.getMainSetDimension());
153 
154     }
155 
156 }