1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18 package org.hipparchus.ode.nonstiff.interpolators;
19
20 import org.hipparchus.ode.EquationsMapper;
21 import org.hipparchus.ode.ODEStateAndDerivative;
22 import org.hipparchus.ode.nonstiff.HighamHall54Integrator;
23
24
25
26
27
28
29
30
31
32 public class HighamHall54StateInterpolator extends RungeKuttaStateInterpolator {
33
34
35 private static final long serialVersionUID = 20111120L;
36
37
38
39
40
41
42
43
44
45
46 public HighamHall54StateInterpolator(final boolean forward,
47 final double[][] yDotK,
48 final ODEStateAndDerivative globalPreviousState,
49 final ODEStateAndDerivative globalCurrentState,
50 final ODEStateAndDerivative softPreviousState,
51 final ODEStateAndDerivative softCurrentState,
52 final EquationsMapper mapper) {
53 super(forward, yDotK, globalPreviousState, globalCurrentState, softPreviousState, softCurrentState, mapper);
54 }
55
56
57 @Override
58 protected HighamHall54StateInterpolator create(final boolean newForward, final double[][] newYDotK,
59 final ODEStateAndDerivative newGlobalPreviousState,
60 final ODEStateAndDerivative newGlobalCurrentState,
61 final ODEStateAndDerivative newSoftPreviousState,
62 final ODEStateAndDerivative newSoftCurrentState,
63 final EquationsMapper newMapper) {
64 return new HighamHall54StateInterpolator(newForward, newYDotK,
65 newGlobalPreviousState, newGlobalCurrentState,
66 newSoftPreviousState, newSoftCurrentState,
67 newMapper);
68 }
69
70
71 @Override
72 protected ODEStateAndDerivative computeInterpolatedStateAndDerivatives(final EquationsMapper mapper,
73 final double time, final double theta,
74 final double thetaH, final double oneMinusThetaH) {
75
76 final double bDot0 = 1 + theta * (-15.0/2.0 + theta * (16.0 - 10.0 * theta));
77 final double bDot1 = 0;
78 final double bDot2 = theta * (459.0/16.0 + theta * (-729.0/8.0 + 135.0/2.0 * theta));
79 final double bDot3 = theta * (-44.0 + theta * (152.0 - 120.0 * theta));
80 final double bDot4 = theta * (375.0/16.0 + theta * (-625.0/8.0 + 125.0/2.0 * theta));
81 final double bDot5 = theta * 5.0/8.0 * (2 * theta - 1);
82
83 final double[] interpolatedState;
84 final double[] interpolatedDerivatives;
85 if (getGlobalPreviousState() != null && theta <= 0.5) {
86 final double b0 = thetaH * (1.0 + theta * (-15.0/4.0 + theta * (16.0/3.0 - 5.0/2.0 * theta)));
87 final double b1 = 0;
88 final double b2 = thetaH * ( theta * (459.0/32.0 + theta * (-243.0/8.0 + theta * 135.0/8.0)));
89 final double b3 = thetaH * ( theta * (-22.0 + theta * (152.0/3.0 + theta * -30.0)));
90 final double b4 = thetaH * ( theta * (375.0/32.0 + theta * (-625.0/24.0 + theta * 125.0/8.0)));
91 final double b5 = thetaH * ( theta * (-5.0/16.0 + theta * 5.0/12.0));
92 interpolatedState = previousStateLinearCombination(b0 , b1, b2, b3, b4, b5);
93 interpolatedDerivatives = derivativeLinearCombination(bDot0 , bDot1, bDot2, bDot3, bDot4, bDot5);
94 } else {
95 final double theta2 = theta * theta;
96 final double h = thetaH / theta;
97 final double b0 = h * (-1.0/12.0 + theta * (1.0 + theta * (-15.0/4.0 + theta * (16.0/3.0 + theta * -5.0/2.0))));
98 final double b1 = 0;
99 final double b2 = h * (-27.0/32.0 + theta2 * (459.0/32.0 + theta * (-243.0/8.0 + theta * 135.0/8.0)));
100 final double b3 = h * (4.0/3.0 + theta2 * (-22.0 + theta * (152.0/3.0 + theta * -30.0)));
101 final double b4 = h * (-125.0/96.0 + theta2 * (375.0/32.0 + theta * (-625.0/24.0 + theta * 125.0/8.0)));
102 final double b5 = h * (-5.0/48.0 + theta2 * (-5.0/16.0 + theta * 5.0/12.0));
103 interpolatedState = currentStateLinearCombination(b0 , b1, b2, b3, b4, b5);
104 interpolatedDerivatives = derivativeLinearCombination(bDot0 , bDot1, bDot2, bDot3, bDot4, bDot5);
105 }
106
107 return mapper.mapStateAndDerivative(time, interpolatedState, interpolatedDerivatives);
108
109 }
110
111 }