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.ClassicalRungeKuttaIntegrator;
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55 public class ClassicalRungeKuttaStateInterpolator extends RungeKuttaStateInterpolator {
56
57
58 private static final long serialVersionUID = 20160328L;
59
60
61
62
63
64
65
66
67
68
69 public ClassicalRungeKuttaStateInterpolator(final boolean forward,
70 final double[][] yDotK,
71 final ODEStateAndDerivative globalPreviousState,
72 final ODEStateAndDerivative globalCurrentState,
73 final ODEStateAndDerivative softPreviousState,
74 final ODEStateAndDerivative softCurrentState,
75 final EquationsMapper mapper) {
76 super(forward, yDotK, globalPreviousState, globalCurrentState, softPreviousState, softCurrentState, mapper);
77 }
78
79
80 @Override
81 protected ClassicalRungeKuttaStateInterpolator create(final boolean newForward, final double[][] newYDotK,
82 final ODEStateAndDerivative newGlobalPreviousState,
83 final ODEStateAndDerivative newGlobalCurrentState,
84 final ODEStateAndDerivative newSoftPreviousState,
85 final ODEStateAndDerivative newSoftCurrentState,
86 final EquationsMapper newMapper) {
87 return new ClassicalRungeKuttaStateInterpolator(newForward, newYDotK,
88 newGlobalPreviousState, newGlobalCurrentState,
89 newSoftPreviousState, newSoftCurrentState,
90 newMapper);
91 }
92
93
94 @Override
95 protected ODEStateAndDerivative computeInterpolatedStateAndDerivatives(final EquationsMapper mapper,
96 final double time, final double theta,
97 final double thetaH, final double oneMinusThetaH) {
98
99 final double oneMinusTheta = 1.0 - theta;
100 final double oneMinus2Theta = 1.0 - theta * 2.0;
101 final double coeffDot1 = oneMinusTheta * oneMinus2Theta;
102 final double coeffDot23 = theta * oneMinusTheta * 2;
103 final double coeffDot4 = -theta * oneMinus2Theta;
104 final double[] interpolatedState;
105 final double[] interpolatedDerivatives;
106
107 if (getGlobalPreviousState() != null && theta <= 0.5) {
108 final double fourTheta2 = theta * theta * 4;
109 final double s = thetaH / 6.0;
110 final double coeff1 = s * (fourTheta2 - theta * 9 + 6);
111 final double coeff23 = s * (theta * 6 - fourTheta2);
112 final double coeff4 = s * (fourTheta2 - theta * 3);
113 interpolatedState = previousStateLinearCombination(coeff1, coeff23, coeff23, coeff4);
114 interpolatedDerivatives = derivativeLinearCombination(coeffDot1, coeffDot23, coeffDot23, coeffDot4);
115 } else {
116 final double fourTheta = theta * 4;
117 final double s = oneMinusThetaH / 6.0;
118 final double coeff1 = s * (theta * (-fourTheta + 5) - 1);
119 final double coeff23 = s * (theta * ( fourTheta - 2) - 2);
120 final double coeff4 = s * (theta * (-fourTheta - 1) - 1);
121 interpolatedState = currentStateLinearCombination(coeff1, coeff23, coeff23, coeff4);
122 interpolatedDerivatives = derivativeLinearCombination(coeffDot1, coeffDot23, coeffDot23, coeffDot4);
123 }
124
125 return mapper.mapStateAndDerivative(time, interpolatedState, interpolatedDerivatives);
126
127 }
128
129 }