1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18 package org.apache.commons.math4.legacy.ode;
19
20 import org.apache.commons.math4.legacy.core.Field;
21 import org.apache.commons.math4.legacy.core.RealFieldElement;
22 import org.apache.commons.math4.legacy.core.MathArrays;
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42 public class TestFieldProblem3<T extends RealFieldElement<T>>
43 extends TestFieldProblemAbstract<T> {
44
45
46 private T e;
47
48
49
50
51
52
53 public TestFieldProblem3(Field<T> field, T e) {
54 super(field);
55 this.e = e;
56 T[] y0 = MathArrays.buildArray(field, 4);
57 y0[0] = e.subtract(1).negate();
58 y0[1] = field.getZero();
59 y0[2] = field.getZero();
60 y0[3] = e.add(1).divide(y0[0]).sqrt();
61 setInitialConditions(convert(0.0), y0);
62 setFinalConditions(convert(20.0));
63 setErrorScale(convert(1.0, 1.0, 1.0, 1.0));
64 }
65
66
67
68
69
70 public TestFieldProblem3(Field<T> field) {
71 this(field, field.getZero().add(0.1));
72 }
73
74 @Override
75 public T[] doComputeDerivatives(T t, T[] y) {
76
77 final T[] yDot = MathArrays.buildArray(getField(), getDimension());
78
79
80 T r2 = y[0].multiply(y[0]).add(y[1].multiply(y[1]));
81 T invR3 = r2.multiply(r2.sqrt()).reciprocal();
82
83
84 yDot[0] = y[2];
85 yDot[1] = y[3];
86 yDot[2] = invR3.negate().multiply(y[0]);
87 yDot[3] = invR3.negate().multiply(y[1]);
88
89 return yDot;
90 }
91
92 @Override
93 public T[] computeTheoreticalState(T t) {
94
95 final T[] y = MathArrays.buildArray(getField(), getDimension());
96
97
98 T E = t;
99 T d = convert(0);
100 T corr = convert(999.0);
101 for (int i = 0; i < 50 && corr.abs().getReal() > 1.0e-12; ++i) {
102 T f2 = e.multiply(E.sin());
103 T f0 = d.subtract(f2);
104 T f1 = e.multiply(E.cos()).subtract(1).negate();
105 T f12 = f1.add(f1);
106 corr = f0.multiply(f12).divide(f1.multiply(f12).subtract(f0.multiply(f2)));
107 d = d.subtract(corr);
108 E = t.add(d);
109 }
110
111 T cosE = E.cos();
112 T sinE = E.sin();
113
114 y[0] = cosE.subtract(e);
115 y[1] = e.multiply(e).subtract(1).negate().sqrt().multiply(sinE);
116 y[2] = sinE.divide(e.multiply(cosE).subtract(1));
117 y[3] = e.multiply(e).subtract(1).negate().sqrt().multiply(cosE).divide(e.multiply(cosE).subtract(1).negate());
118
119 return y;
120 }
121 }