diff --git a/src/main/java/com/thealgorithms/physics/SimplePendulumRK4.java b/src/main/java/com/thealgorithms/physics/SimplePendulumRK4.java new file mode 100644 index 000000000000..6de69c103b5a --- /dev/null +++ b/src/main/java/com/thealgorithms/physics/SimplePendulumRK4.java @@ -0,0 +1,126 @@ +package com.thealgorithms.physics; + +/** + * Simulates a simple pendulum using the Runge-Kutta 4th order method. + * The pendulum is modeled with the nonlinear differential equation. + * + * @author [Yash Rajput](https://github.com/the-yash-rajput) + */ +public final class SimplePendulumRK4 { + + private SimplePendulumRK4() { + throw new AssertionError("No instances."); + } + + private final double length; // meters + private final double g; // acceleration due to gravity (m/s^2) + + /** + * Constructs a simple pendulum simulator. + * + * @param length the length of the pendulum in meters + * @param g the acceleration due to gravity in m/s^2 + */ + public SimplePendulumRK4(double length, double g) { + if (length <= 0) { + throw new IllegalArgumentException("Length must be positive"); + } + if (g <= 0) { + throw new IllegalArgumentException("Gravity must be positive"); + } + this.length = length; + this.g = g; + } + + /** + * Computes the derivatives of the state vector. + * State: [theta, omega] where theta is angle and omega is angular velocity. + * + * @param state the current state [theta, omega] + * @return the derivatives [dtheta/dt, domega/dt] + */ + private double[] derivatives(double[] state) { + double theta = state[0]; + double omega = state[1]; + double dtheta = omega; + double domega = -(g / length) * Math.sin(theta); + return new double[] {dtheta, domega}; + } + + /** + * Performs one time step using the RK4 method. + * + * @param state the current state [theta, omega] + * @param dt the time step size + * @return the new state after time dt + */ + public double[] stepRK4(double[] state, double dt) { + if (state == null || state.length != 2) { + throw new IllegalArgumentException("State must be array of length 2"); + } + if (dt <= 0) { + throw new IllegalArgumentException("Time step must be positive"); + } + + double[] k1 = derivatives(state); + double[] s2 = new double[] {state[0] + 0.5 * dt * k1[0], state[1] + 0.5 * dt * k1[1]}; + + double[] k2 = derivatives(s2); + double[] s3 = new double[] {state[0] + 0.5 * dt * k2[0], state[1] + 0.5 * dt * k2[1]}; + + double[] k3 = derivatives(s3); + double[] s4 = new double[] {state[0] + dt * k3[0], state[1] + dt * k3[1]}; + + double[] k4 = derivatives(s4); + + double thetaNext = state[0] + dt / 6.0 * (k1[0] + 2 * k2[0] + 2 * k3[0] + k4[0]); + double omegaNext = state[1] + dt / 6.0 * (k1[1] + 2 * k2[1] + 2 * k3[1] + k4[1]); + + return new double[] {thetaNext, omegaNext}; + } + + /** + * Simulates the pendulum for a given duration. + * + * @param initialState the initial state [theta, omega] + * @param dt the time step size + * @param steps the number of steps to simulate + * @return array of states at each step + */ + public double[][] simulate(double[] initialState, double dt, int steps) { + double[][] trajectory = new double[steps + 1][2]; + trajectory[0] = initialState.clone(); + + double[] currentState = initialState.clone(); + for (int i = 1; i <= steps; i++) { + currentState = stepRK4(currentState, dt); + trajectory[i] = currentState.clone(); + } + + return trajectory; + } + + /** + * Calculates the total energy of the pendulum. + * E = (1/2) * m * L^2 * omega^2 + m * g * L * (1 - cos(theta)) + * We use m = 1 for simplicity. + * + * @param state the current state [theta, omega] + * @return the total energy + */ + public double calculateEnergy(double[] state) { + double theta = state[0]; + double omega = state[1]; + double kineticEnergy = 0.5 * length * length * omega * omega; + double potentialEnergy = g * length * (1 - Math.cos(theta)); + return kineticEnergy + potentialEnergy; + } + + public double getLength() { + return length; + } + + public double getGravity() { + return g; + } +} diff --git a/src/test/java/com/thealgorithms/physics/SimplePendulumRK4Test.java b/src/test/java/com/thealgorithms/physics/SimplePendulumRK4Test.java new file mode 100644 index 000000000000..bec4ad8cbbd7 --- /dev/null +++ b/src/test/java/com/thealgorithms/physics/SimplePendulumRK4Test.java @@ -0,0 +1,261 @@ +package com.thealgorithms.physics; + +import org.junit.jupiter.api.Assertions; +import org.junit.jupiter.api.DisplayName; +import org.junit.jupiter.api.Test; + +/** + * Test class for SimplePendulumRK4. + * Tests numerical accuracy, physical correctness, and edge cases. + */ +class SimplePendulumRK4Test { + + private static final double EPSILON = 1e-6; + private static final double ENERGY_DRIFT_TOLERANCE = 1e-3; + @Test + @DisplayName("Test constructor creates valid pendulum") + void testConstructor() { + SimplePendulumRK4 pendulum = new SimplePendulumRK4(1.5, 9.81); + Assertions.assertNotNull(pendulum); + Assertions.assertEquals(1.5, pendulum.getLength(), EPSILON); + Assertions.assertEquals(9.81, pendulum.getGravity(), EPSILON); + } + + @Test + @DisplayName("Test constructor rejects negative length") + void testConstructorNegativeLength() { + Assertions.assertThrows(IllegalArgumentException.class, () -> { new SimplePendulumRK4(-1.0, 9.81); }); + } + + @Test + @DisplayName("Test constructor rejects negative gravity") + void testConstructorNegativeGravity() { + Assertions.assertThrows(IllegalArgumentException.class, () -> { new SimplePendulumRK4(1.0, -9.81); }); + } + + @Test + @DisplayName("Test constructor rejects zero length") + void testConstructorZeroLength() { + Assertions.assertThrows(IllegalArgumentException.class, () -> { new SimplePendulumRK4(0.0, 9.81); }); + } + + @Test + @DisplayName("Test getters return correct values") + void testGetters() { + SimplePendulumRK4 pendulum = new SimplePendulumRK4(2.5, 10.0); + Assertions.assertEquals(2.5, pendulum.getLength(), EPSILON); + Assertions.assertEquals(10.0, pendulum.getGravity(), EPSILON); + } + + @Test + @DisplayName("Test single RK4 step returns valid state") + void testSingleStep() { + SimplePendulumRK4 pendulum = new SimplePendulumRK4(1.0, 9.81); + double[] state = {0.1, 0.0}; + double[] newState = pendulum.stepRK4(state, 0.01); + + Assertions.assertNotNull(newState); + Assertions.assertEquals(2, newState.length); + } + + @Test + @DisplayName("Test equilibrium stability (pendulum at rest stays at rest)") + void testEquilibrium() { + SimplePendulumRK4 pendulum = new SimplePendulumRK4(1.0, 9.81); + double[] state = {0.0, 0.0}; + + for (int i = 0; i < 100; i++) { + state = pendulum.stepRK4(state, 0.01); + } + + Assertions.assertEquals(0.0, state[0], EPSILON, "Theta should remain at equilibrium"); + Assertions.assertEquals(0.0, state[1], EPSILON, "Omega should remain zero"); + } + + @Test + @DisplayName("Test small angle oscillation returns to initial position") + void testSmallAngleOscillation() { + SimplePendulumRK4 pendulum = new SimplePendulumRK4(1.0, 9.81); + double initialAngle = Math.toRadians(5.0); + double[] state = {initialAngle, 0.0}; + double dt = 0.01; + + // Theoretical period for small angles + double expectedPeriod = 2 * Math.PI * Math.sqrt(1.0 / 9.81); + int stepsPerPeriod = (int) (expectedPeriod / dt); + + double[][] trajectory = pendulum.simulate(state, dt, stepsPerPeriod); + double finalTheta = trajectory[stepsPerPeriod][0]; + + // After one period, should return close to initial position + double error = Math.abs(finalTheta - initialAngle) / Math.abs(initialAngle); + Assertions.assertTrue(error < 0.05, "Small angle approximation error should be < 5%"); + } + + @Test + @DisplayName("Test large angle oscillation is symmetric") + void testLargeAngleOscillation() { + SimplePendulumRK4 pendulum = new SimplePendulumRK4(1.0, 9.81); + double[] state = {Math.toRadians(120.0), 0.0}; + + double[][] trajectory = pendulum.simulate(state, 0.01, 500); + + double maxTheta = Double.NEGATIVE_INFINITY; + double minTheta = Double.POSITIVE_INFINITY; + + for (double[] s : trajectory) { + maxTheta = Math.max(maxTheta, s[0]); + minTheta = Math.min(minTheta, s[0]); + } + + Assertions.assertTrue(maxTheta > 0, "Should have positive excursions"); + Assertions.assertTrue(minTheta < 0, "Should have negative excursions"); + + // Check symmetry + double asymmetry = Math.abs((maxTheta + minTheta) / maxTheta); + Assertions.assertTrue(asymmetry < 0.1, "Oscillation should be symmetric"); + } + + @Test + @DisplayName("Test energy conservation for small angle") + void testEnergyConservationSmallAngle() { + SimplePendulumRK4 pendulum = new SimplePendulumRK4(1.0, 9.81); + double[] state = {Math.toRadians(15.0), 0.0}; + + double initialEnergy = pendulum.calculateEnergy(state); + + for (int i = 0; i < 1000; i++) { + state = pendulum.stepRK4(state, 0.01); + } + + double finalEnergy = pendulum.calculateEnergy(state); + double drift = Math.abs(finalEnergy - initialEnergy) / initialEnergy; + + Assertions.assertTrue(drift < ENERGY_DRIFT_TOLERANCE, "Energy drift should be < 0.1%, got: " + (drift * 100) + "%"); + } + + @Test + @DisplayName("Test energy conservation for large angle") + void testEnergyConservationLargeAngle() { + SimplePendulumRK4 pendulum = new SimplePendulumRK4(1.0, 9.81); + double[] state = {Math.toRadians(90.0), 0.0}; + + double initialEnergy = pendulum.calculateEnergy(state); + + for (int i = 0; i < 1000; i++) { + state = pendulum.stepRK4(state, 0.01); + } + + double finalEnergy = pendulum.calculateEnergy(state); + double drift = Math.abs(finalEnergy - initialEnergy) / initialEnergy; + + Assertions.assertTrue(drift < ENERGY_DRIFT_TOLERANCE, "Energy drift should be < 0.1%, got: " + (drift * 100) + "%"); + } + + @Test + @DisplayName("Test simulate method returns correct trajectory") + void testSimulate() { + SimplePendulumRK4 pendulum = new SimplePendulumRK4(1.0, 9.81); + double[] initialState = {Math.toRadians(20.0), 0.0}; + int steps = 100; + + double[][] trajectory = pendulum.simulate(initialState, 0.01, steps); + + Assertions.assertEquals(steps + 1, trajectory.length, "Trajectory should have steps + 1 entries"); + Assertions.assertArrayEquals(initialState, trajectory[0], EPSILON, "First entry should match initial state"); + + // Verify state changes over time + boolean changed = false; + for (int i = 1; i <= steps; i++) { + if (Math.abs(trajectory[i][0] - initialState[0]) > EPSILON) { + changed = true; + break; + } + } + Assertions.assertTrue(changed, "Simulation should progress from initial state"); + } + + @Test + @DisplayName("Test energy calculation at equilibrium") + void testEnergyAtEquilibrium() { + SimplePendulumRK4 pendulum = new SimplePendulumRK4(1.0, 9.81); + double[] state = {0.0, 0.0}; + double energy = pendulum.calculateEnergy(state); + Assertions.assertEquals(0.0, energy, EPSILON, "Energy at equilibrium should be zero"); + } + + @Test + @DisplayName("Test energy calculation at maximum angle") + void testEnergyAtMaxAngle() { + SimplePendulumRK4 pendulum = new SimplePendulumRK4(1.0, 9.81); + double[] state = {Math.PI / 2, 0.0}; + double energy = pendulum.calculateEnergy(state); + Assertions.assertTrue(energy > 0, "Energy should be positive at max angle"); + } + + @Test + @DisplayName("Test energy calculation with angular velocity") + void testEnergyWithVelocity() { + SimplePendulumRK4 pendulum = new SimplePendulumRK4(1.0, 9.81); + double[] state = {0.0, 1.0}; + double energy = pendulum.calculateEnergy(state); + Assertions.assertTrue(energy > 0, "Energy should be positive with velocity"); + } + + @Test + @DisplayName("Test stepRK4 rejects null state") + void testStepRejectsNullState() { + SimplePendulumRK4 pendulum = new SimplePendulumRK4(1.0, 9.81); + Assertions.assertThrows(IllegalArgumentException.class, () -> { pendulum.stepRK4(null, 0.01); }); + } + + @Test + @DisplayName("Test stepRK4 rejects invalid state length") + void testStepRejectsInvalidStateLength() { + SimplePendulumRK4 pendulum = new SimplePendulumRK4(1.0, 9.81); + Assertions.assertThrows(IllegalArgumentException.class, () -> { pendulum.stepRK4(new double[] {0.1}, 0.01); }); + } + + @Test + @DisplayName("Test stepRK4 rejects negative time step") + void testStepRejectsNegativeTimeStep() { + SimplePendulumRK4 pendulum = new SimplePendulumRK4(1.0, 9.81); + Assertions.assertThrows(IllegalArgumentException.class, () -> { pendulum.stepRK4(new double[] {0.1, 0.2}, -0.01); }); + } + + @Test + @DisplayName("Test extreme condition: very large angle") + void testExtremeLargeAngle() { + SimplePendulumRK4 pendulum = new SimplePendulumRK4(1.0, 9.81); + double[] state = {Math.toRadians(179.0), 0.0}; + double[] result = pendulum.stepRK4(state, 0.01); + + Assertions.assertNotNull(result); + Assertions.assertTrue(Double.isFinite(result[0]), "Should handle large angles without NaN"); + Assertions.assertTrue(Double.isFinite(result[1]), "Should handle large angles without NaN"); + } + + @Test + @DisplayName("Test extreme condition: high angular velocity") + void testExtremeHighVelocity() { + SimplePendulumRK4 pendulum = new SimplePendulumRK4(1.0, 9.81); + double[] state = {0.0, 10.0}; + double[] result = pendulum.stepRK4(state, 0.01); + + Assertions.assertNotNull(result); + Assertions.assertTrue(Double.isFinite(result[0]), "Should handle high velocity without NaN"); + Assertions.assertTrue(Double.isFinite(result[1]), "Should handle high velocity without NaN"); + } + + @Test + @DisplayName("Test extreme condition: very small time step") + void testExtremeSmallTimeStep() { + SimplePendulumRK4 pendulum = new SimplePendulumRK4(1.0, 9.81); + double[] state = {Math.toRadians(10.0), 0.0}; + double[] result = pendulum.stepRK4(state, 1e-6); + + Assertions.assertNotNull(result); + Assertions.assertTrue(Double.isFinite(result[0]), "Should handle small time steps without NaN"); + Assertions.assertTrue(Double.isFinite(result[1]), "Should handle small time steps without NaN"); + } +}