Skip to content

Commit e1773e9

Browse files
Adding DampedOscillator code (TheAlgorithms#6801)
* Adding DampedOscillator code * Adding one more test case * Adding one more test case * Adding one more test case * Fixing build issues. * Fixing build issues.
1 parent fb5a765 commit e1773e9

File tree

2 files changed

+252
-0
lines changed

2 files changed

+252
-0
lines changed
Lines changed: 109 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,109 @@
1+
package com.thealgorithms.physics;
2+
3+
/**
4+
* Models a damped harmonic oscillator, capturing the behavior of a mass-spring-damper system.
5+
*
6+
* <p>The system is defined by the second-order differential equation:
7+
* x'' + 2 * gamma * x' + omega₀² * x = 0
8+
* where:
9+
* <ul>
10+
* <li><b>omega₀</b> is the natural (undamped) angular frequency in radians per second.</li>
11+
* <li><b>gamma</b> is the damping coefficient in inverse seconds.</li>
12+
* </ul>
13+
*
14+
* <p>This implementation provides:
15+
* <ul>
16+
* <li>An analytical solution for the underdamped case (γ < ω₀).</li>
17+
* <li>A numerical integrator based on the explicit Euler method for simulation purposes.</li>
18+
* </ul>
19+
*
20+
* <p><strong>Usage Example:</strong>
21+
* <pre>{@code
22+
* DampedOscillator oscillator = new DampedOscillator(10.0, 0.5);
23+
* double displacement = oscillator.displacementAnalytical(1.0, 0.0, 0.1);
24+
* double[] nextState = oscillator.stepEuler(new double[]{1.0, 0.0}, 0.001);
25+
* }</pre>
26+
*
27+
* @author [Yash Rajput](https://github.com/the-yash-rajput)
28+
*/
29+
public final class DampedOscillator {
30+
31+
/** Natural (undamped) angular frequency (rad/s). */
32+
private final double omega0;
33+
34+
/** Damping coefficient (s⁻¹). */
35+
private final double gamma;
36+
37+
private DampedOscillator() {
38+
throw new AssertionError("No instances.");
39+
}
40+
41+
/**
42+
* Constructs a damped oscillator model.
43+
*
44+
* @param omega0 the natural frequency (rad/s), must be positive
45+
* @param gamma the damping coefficient (s⁻¹), must be non-negative
46+
* @throws IllegalArgumentException if parameters are invalid
47+
*/
48+
public DampedOscillator(double omega0, double gamma) {
49+
if (omega0 <= 0) {
50+
throw new IllegalArgumentException("Natural frequency must be positive.");
51+
}
52+
if (gamma < 0) {
53+
throw new IllegalArgumentException("Damping coefficient must be non-negative.");
54+
}
55+
this.omega0 = omega0;
56+
this.gamma = gamma;
57+
}
58+
59+
/**
60+
* Computes the analytical displacement of an underdamped oscillator.
61+
* Formula: x(t) = A * exp(-γt) * cos(ω_d t + φ)
62+
*
63+
* @param amplitude the initial amplitude A
64+
* @param phase the initial phase φ (radians)
65+
* @param time the time t (seconds)
66+
* @return the displacement x(t)
67+
*/
68+
public double displacementAnalytical(double amplitude, double phase, double time) {
69+
double omegaD = Math.sqrt(Math.max(0.0, omega0 * omega0 - gamma * gamma));
70+
return amplitude * Math.exp(-gamma * time) * Math.cos(omegaD * time + phase);
71+
}
72+
73+
/**
74+
* Performs a single integration step using the explicit Euler method.
75+
* State vector format: [x, v], where v = dx/dt.
76+
*
77+
* @param state the current state [x, v]
78+
* @param dt the time step (seconds)
79+
* @return the next state [x_next, v_next]
80+
* @throws IllegalArgumentException if the state array is invalid or dt is non-positive
81+
*/
82+
public double[] stepEuler(double[] state, double dt) {
83+
if (state == null || state.length != 2) {
84+
throw new IllegalArgumentException("State must be a non-null array of length 2.");
85+
}
86+
if (dt <= 0) {
87+
throw new IllegalArgumentException("Time step must be positive.");
88+
}
89+
90+
double x = state[0];
91+
double v = state[1];
92+
double acceleration = -2.0 * gamma * v - omega0 * omega0 * x;
93+
94+
double xNext = x + dt * v;
95+
double vNext = v + dt * acceleration;
96+
97+
return new double[] {xNext, vNext};
98+
}
99+
100+
/** @return the natural (undamped) angular frequency (rad/s). */
101+
public double getOmega0() {
102+
return omega0;
103+
}
104+
105+
/** @return the damping coefficient (s⁻¹). */
106+
public double getGamma() {
107+
return gamma;
108+
}
109+
}
Lines changed: 143 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,143 @@
1+
package com.thealgorithms.physics;
2+
3+
import static org.junit.jupiter.api.Assertions.assertAll;
4+
import static org.junit.jupiter.api.Assertions.assertEquals;
5+
import static org.junit.jupiter.api.Assertions.assertThrows;
6+
7+
import org.junit.jupiter.api.DisplayName;
8+
import org.junit.jupiter.api.Test;
9+
10+
/**
11+
* Unit tests for {@link DampedOscillator}.
12+
*
13+
* <p>Tests focus on:
14+
* <ul>
15+
* <li>Constructor validation</li>
16+
* <li>Analytical displacement for underdamped and overdamped parameterizations</li>
17+
* <li>Basic numeric integration sanity using explicit Euler for small step sizes</li>
18+
* <li>Method argument validation (null/invalid inputs)</li>
19+
* </ul>
20+
*/
21+
@DisplayName("DampedOscillator — unit tests")
22+
public class DampedOscillatorTest {
23+
24+
private static final double TOLERANCE = 1e-3;
25+
26+
@Test
27+
@DisplayName("Constructor rejects invalid parameters")
28+
void constructorValidation() {
29+
assertAll("invalid-constructor-params",
30+
()
31+
-> assertThrows(IllegalArgumentException.class, () -> new DampedOscillator(0.0, 0.1), "omega0 == 0 should throw"),
32+
() -> assertThrows(IllegalArgumentException.class, () -> new DampedOscillator(-1.0, 0.1), "negative omega0 should throw"), () -> assertThrows(IllegalArgumentException.class, () -> new DampedOscillator(1.0, -0.1), "negative gamma should throw"));
33+
}
34+
35+
@Test
36+
@DisplayName("Analytical displacement matches expected formula for underdamped case")
37+
void analyticalUnderdamped() {
38+
double omega0 = 10.0;
39+
double gamma = 0.5;
40+
DampedOscillator d = new DampedOscillator(omega0, gamma);
41+
42+
double a = 1.0;
43+
double phi = 0.2;
44+
double t = 0.123;
45+
46+
// expected: a * exp(-gamma * t) * cos(omega_d * t + phi)
47+
double omegaD = Math.sqrt(Math.max(0.0, omega0 * omega0 - gamma * gamma));
48+
double expected = a * Math.exp(-gamma * t) * Math.cos(omegaD * t + phi);
49+
50+
double actual = d.displacementAnalytical(a, phi, t);
51+
assertEquals(expected, actual, 1e-12, "Analytical underdamped displacement should match closed-form value");
52+
}
53+
54+
@Test
55+
@DisplayName("Analytical displacement gracefully handles overdamped parameters (omegaD -> 0)")
56+
void analyticalOverdamped() {
57+
double omega0 = 1.0;
58+
double gamma = 2.0; // gamma > omega0 => omega_d = 0 in our implementation (Math.max)
59+
DampedOscillator d = new DampedOscillator(omega0, gamma);
60+
61+
double a = 2.0;
62+
double phi = Math.PI / 4.0;
63+
double t = 0.5;
64+
65+
// With omegaD forced to 0 by implementation, expected simplifies to:
66+
double expected = a * Math.exp(-gamma * t) * Math.cos(phi);
67+
double actual = d.displacementAnalytical(a, phi, t);
68+
69+
assertEquals(expected, actual, 1e-12, "Overdamped handling should reduce to exponential * cos(phase)");
70+
}
71+
72+
@Test
73+
@DisplayName("Explicit Euler step approximates analytical solution for small dt over short time")
74+
void eulerApproximatesAnalyticalSmallDt() {
75+
double omega0 = 10.0;
76+
double gamma = 0.5;
77+
DampedOscillator d = new DampedOscillator(omega0, gamma);
78+
79+
double a = 1.0;
80+
double phi = 0.0;
81+
82+
// initial conditions consistent with amplitude a and zero phase:
83+
// x(0) = a, v(0) = -a * gamma * cos(phi) + a * omegaD * sin(phi)
84+
double omegaD = Math.sqrt(Math.max(0.0, omega0 * omega0 - gamma * gamma));
85+
double x0 = a * Math.cos(phi);
86+
double v0 = -a * gamma * Math.cos(phi) - a * omegaD * Math.sin(phi); // small general form
87+
88+
double dt = 1e-4;
89+
int steps = 1000; // simulate to t = 0.1s
90+
double tFinal = steps * dt;
91+
92+
double[] state = new double[] {x0, v0};
93+
for (int i = 0; i < steps; i++) {
94+
state = d.stepEuler(state, dt);
95+
}
96+
97+
double analyticAtT = d.displacementAnalytical(a, phi, tFinal);
98+
double numericAtT = state[0];
99+
100+
// Euler is low-order — allow a small tolerance but assert it remains close for small dt + short time.
101+
assertEquals(analyticAtT, numericAtT, TOLERANCE, String.format("Numeric Euler should approximate analytical solution at t=%.6f (tolerance=%g)", tFinal, TOLERANCE));
102+
}
103+
104+
@Test
105+
@DisplayName("stepEuler validates inputs and throws on null/invalid dt/state")
106+
void eulerInputValidation() {
107+
DampedOscillator d = new DampedOscillator(5.0, 0.1);
108+
109+
assertAll("invalid-stepEuler-args",
110+
()
111+
-> assertThrows(IllegalArgumentException.class, () -> d.stepEuler(null, 0.01), "null state should throw"),
112+
()
113+
-> assertThrows(IllegalArgumentException.class, () -> d.stepEuler(new double[] {1.0}, 0.01), "state array with invalid length should throw"),
114+
() -> assertThrows(IllegalArgumentException.class, () -> d.stepEuler(new double[] {1.0, 0.0}, 0.0), "non-positive dt should throw"), () -> assertThrows(IllegalArgumentException.class, () -> d.stepEuler(new double[] {1.0, 0.0}, -1e-3), "negative dt should throw"));
115+
}
116+
117+
@Test
118+
@DisplayName("Getter methods return configured parameters")
119+
void gettersReturnConfiguration() {
120+
double omega0 = Math.PI;
121+
double gamma = 0.01;
122+
DampedOscillator d = new DampedOscillator(omega0, gamma);
123+
124+
assertAll("getters", () -> assertEquals(omega0, d.getOmega0(), 0.0, "getOmega0 should return configured omega0"), () -> assertEquals(gamma, d.getGamma(), 0.0, "getGamma should return configured gamma"));
125+
}
126+
127+
@Test
128+
@DisplayName("Analytical displacement at t=0 returns initial amplitude * cos(phase)")
129+
void analyticalAtZeroTime() {
130+
double omega0 = 5.0;
131+
double gamma = 0.2;
132+
DampedOscillator d = new DampedOscillator(omega0, gamma);
133+
134+
double a = 2.0;
135+
double phi = Math.PI / 3.0;
136+
double t = 0.0;
137+
138+
double expected = a * Math.cos(phi);
139+
double actual = d.displacementAnalytical(a, phi, t);
140+
141+
assertEquals(expected, actual, 1e-12, "Displacement at t=0 should be a * cos(phase)");
142+
}
143+
}

0 commit comments

Comments
 (0)