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