Added SimplePendulumRK4

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;
+    }
+}

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");
+    }
+}