1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23 package org.hipparchus.ode.nonstiff;
24
25 import org.hipparchus.CalculusFieldElement;
26 import org.hipparchus.Field;
27 import org.hipparchus.ode.FieldEquationsMapper;
28 import org.hipparchus.ode.FieldODEStateAndDerivative;
29 import org.hipparchus.util.FastMath;
30 import org.hipparchus.util.MathArrays;
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
56
57
58 public class DormandPrince54FieldIntegrator<T extends CalculusFieldElement<T>>
59 extends EmbeddedRungeKuttaFieldIntegrator<T> {
60
61
62
63
64
65
66
67
68
69
70
71
72
73 public DormandPrince54FieldIntegrator(final Field<T> field,
74 final double minStep, final double maxStep,
75 final double scalAbsoluteTolerance,
76 final double scalRelativeTolerance) {
77 super(field, DormandPrince54Integrator.METHOD_NAME, 6,
78 minStep, maxStep, scalAbsoluteTolerance, scalRelativeTolerance);
79 }
80
81
82
83
84
85
86
87
88
89
90
91
92
93 public DormandPrince54FieldIntegrator(final Field<T> field,
94 final double minStep, final double maxStep,
95 final double[] vecAbsoluteTolerance,
96 final double[] vecRelativeTolerance) {
97 super(field, DormandPrince54Integrator.METHOD_NAME, 6,
98 minStep, maxStep, vecAbsoluteTolerance, vecRelativeTolerance);
99 }
100
101
102 @Override
103 public T[] getC() {
104 final T[] c = MathArrays.buildArray(getField(), 6);
105 c[0] = fraction(1, 5);
106 c[1] = fraction(3, 10);
107 c[2] = fraction(4, 5);
108 c[3] = fraction(8, 9);
109 c[4] = getField().getOne();
110 c[5] = getField().getOne();
111 return c;
112 }
113
114
115 @Override
116 public T[][] getA() {
117 final T[][] a = MathArrays.buildArray(getField(), 6, -1);
118 for (int i = 0; i < a.length; ++i) {
119 a[i] = MathArrays.buildArray(getField(), i + 1);
120 }
121 a[0][0] = fraction( 1, 5);
122 a[1][0] = fraction( 3, 40);
123 a[1][1] = fraction( 9, 40);
124 a[2][0] = fraction( 44, 45);
125 a[2][1] = fraction( -56, 15);
126 a[2][2] = fraction( 32, 9);
127 a[3][0] = fraction( 19372, 6561);
128 a[3][1] = fraction(-25360, 2187);
129 a[3][2] = fraction( 64448, 6561);
130 a[3][3] = fraction( -212, 729);
131 a[4][0] = fraction( 9017, 3168);
132 a[4][1] = fraction( -355, 33);
133 a[4][2] = fraction( 46732, 5247);
134 a[4][3] = fraction( 49, 176);
135 a[4][4] = fraction( -5103, 18656);
136 a[5][0] = fraction( 35, 384);
137 a[5][1] = getField().getZero();
138 a[5][2] = fraction( 500, 1113);
139 a[5][3] = fraction( 125, 192);
140 a[5][4] = fraction( -2187, 6784);
141 a[5][5] = fraction( 11, 84);
142 return a;
143 }
144
145
146 @Override
147 public T[] getB() {
148 final T[] b = MathArrays.buildArray(getField(), 7);
149 b[0] = fraction( 35, 384);
150 b[1] = getField().getZero();
151 b[2] = fraction( 500, 1113);
152 b[3] = fraction( 125, 192);
153 b[4] = fraction(-2187, 6784);
154 b[5] = fraction( 11, 84);
155 b[6] = getField().getZero();
156 return b;
157 }
158
159
160 @Override
161 protected DormandPrince54FieldStateInterpolator<T>
162 createInterpolator(final boolean forward, T[][] yDotK,
163 final FieldODEStateAndDerivative<T> globalPreviousState,
164 final FieldODEStateAndDerivative<T> globalCurrentState, final FieldEquationsMapper<T> mapper) {
165 return new DormandPrince54FieldStateInterpolator<T>(getField(), forward, yDotK,
166 globalPreviousState, globalCurrentState,
167 globalPreviousState, globalCurrentState,
168 mapper);
169 }
170
171
172 @Override
173 public int getOrder() {
174 return 5;
175 }
176
177
178 @Override
179 protected double estimateError(final T[][] yDotK, final T[] y0, final T[] y1, final T h) {
180
181 final StepsizeHelper helper = getStepSizeHelper();
182 double error = 0;
183
184 for (int j = 0; j < helper.getMainSetDimension(); ++j) {
185 final double errSum = DormandPrince54Integrator.E1 * yDotK[0][j].getReal() + DormandPrince54Integrator.E3 * yDotK[2][j].getReal() +
186 DormandPrince54Integrator.E4 * yDotK[3][j].getReal() + DormandPrince54Integrator.E5 * yDotK[4][j].getReal() +
187 DormandPrince54Integrator.E6 * yDotK[5][j].getReal() + DormandPrince54Integrator.E7 * yDotK[6][j].getReal();
188 final double tol = helper.getTolerance(j, FastMath.max(FastMath.abs(y0[j].getReal()), FastMath.abs(y1[j].getReal())));
189 final double ratio = h.getReal() * errSum / tol;
190 error += ratio * ratio;
191 }
192
193 return FastMath.sqrt(error / helper.getMainSetDimension());
194
195 }
196
197 }