mirror of
https://github.com/TheAlgorithms/Java.git
synced 2026-03-13 08:40:43 +08:00
Added SimplePendulumRK4 (#6800)
* Added SimplePendulumRK4 * Fixed build issue. * Fixed build issue. * Fixed build issue.
This commit is contained in:
Yash Rajput
126
src/main/java/com/thealgorithms/physics/SimplePendulumRK4.java
Normal file
126
src/main/java/com/thealgorithms/physics/SimplePendulumRK4.java
Normal file
@@ -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;
|
||||
}
|
||||
}
|
||||
@@ -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");
|
||||
}
|
||||
}
|
||||
Reference in New Issue
Block a user