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 /** Name of integration scheme. */
40 public 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 }