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.analysis.function;
24
25 import org.hipparchus.analysis.UnivariateFunction;
26 import org.hipparchus.analysis.differentiation.DSFactory;
27 import org.hipparchus.analysis.differentiation.DerivativeStructure;
28 import org.hipparchus.exception.MathIllegalArgumentException;
29 import org.hipparchus.exception.NullArgumentException;
30 import org.hipparchus.util.FastMath;
31 import org.hipparchus.util.Precision;
32 import org.junit.jupiter.api.Test;
33
34 import static org.junit.jupiter.api.Assertions.assertEquals;
35 import static org.junit.jupiter.api.Assertions.assertThrows;
36
37
38
39
40 class HarmonicOscillatorTest {
41 private final double EPS = Math.ulp(1d);
42
43 @Test
44 void testSomeValues() {
45 final double a = -1.2;
46 final double w = 0.34;
47 final double p = 5.6;
48 final UnivariateFunction f = new HarmonicOscillator(a, w, p);
49
50 final double d = 0.12345;
51 for (int i = 0; i < 10; i++) {
52 final double v = i * d;
53 assertEquals(a * FastMath.cos(w * v + p), f.value(v), 0);
54 }
55 }
56
57 @Test
58 void testDerivative() {
59 final double a = -1.2;
60 final double w = 0.34;
61 final double p = 5.6;
62 final HarmonicOscillator f = new HarmonicOscillator(a, w, p);
63
64 for (int maxOrder = 0; maxOrder < 6; ++maxOrder) {
65 final DSFactory factory = new DSFactory(1, maxOrder);
66 final double d = 0.12345;
67 for (int i = 0; i < 10; i++) {
68 final double v = i * d;
69 final DerivativeStructure h = f.value(factory.variable(0, v));
70 for (int k = 0; k <= maxOrder; ++k) {
71 final double trigo;
72 switch (k % 4) {
73 case 0:
74 trigo = +FastMath.cos(w * v + p);
75 break;
76 case 1:
77 trigo = -FastMath.sin(w * v + p);
78 break;
79 case 2:
80 trigo = -FastMath.cos(w * v + p);
81 break;
82 default:
83 trigo = +FastMath.sin(w * v + p);
84 break;
85 }
86 assertEquals(a * FastMath.pow(w, k) * trigo,
87 h.getPartialDerivative(k),
88 Precision.EPSILON);
89 }
90 }
91 }
92 }
93
94 @Test
95 void testParametricUsage1() {
96 assertThrows(NullArgumentException.class, () -> {
97 final HarmonicOscillator.Parametric g = new HarmonicOscillator.Parametric();
98 g.value(0, null);
99 });
100 }
101
102 @Test
103 void testParametricUsage2() {
104 assertThrows(MathIllegalArgumentException.class, () -> {
105 final HarmonicOscillator.Parametric g = new HarmonicOscillator.Parametric();
106 g.value(0, new double[]{0});
107 });
108 }
109
110 @Test
111 void testParametricUsage3() {
112 assertThrows(NullArgumentException.class, () -> {
113 final HarmonicOscillator.Parametric g = new HarmonicOscillator.Parametric();
114 g.gradient(0, null);
115 });
116 }
117
118 @Test
119 void testParametricUsage4() {
120 assertThrows(MathIllegalArgumentException.class, () -> {
121 final HarmonicOscillator.Parametric g = new HarmonicOscillator.Parametric();
122 g.gradient(0, new double[]{0});
123 });
124 }
125
126 @Test
127 void testParametricValue() {
128 final double amplitude = 2;
129 final double omega = 3;
130 final double phase = 4;
131 final HarmonicOscillator f = new HarmonicOscillator(amplitude, omega, phase);
132
133 final HarmonicOscillator.Parametric g = new HarmonicOscillator.Parametric();
134 assertEquals(f.value(-1), g.value(-1, new double[] {amplitude, omega, phase}), 0);
135 assertEquals(f.value(0), g.value(0, new double[] {amplitude, omega, phase}), 0);
136 assertEquals(f.value(2), g.value(2, new double[] {amplitude, omega, phase}), 0);
137 }
138
139 @Test
140 void testParametricGradient() {
141 final double amplitude = 2;
142 final double omega = 3;
143 final double phase = 4;
144 final HarmonicOscillator.Parametric f = new HarmonicOscillator.Parametric();
145
146 final double x = 1;
147 final double[] grad = f.gradient(1, new double[] {amplitude, omega, phase});
148 final double xTimesOmegaPlusPhase = omega * x + phase;
149 final double a = FastMath.cos(xTimesOmegaPlusPhase);
150 assertEquals(a, grad[0], EPS);
151 final double w = -amplitude * x * FastMath.sin(xTimesOmegaPlusPhase);
152 assertEquals(w, grad[1], EPS);
153 final double p = -amplitude * FastMath.sin(xTimesOmegaPlusPhase);
154 assertEquals(p, grad[2], EPS);
155 }
156 }