Skip to content

Commit 29c18a3

Browse files
Added SimplePendulumRK4
1 parent 3519a91 commit 29c18a3

File tree

2 files changed

+478
-0
lines changed

2 files changed

+478
-0
lines changed
Lines changed: 124 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,124 @@
1+
package com.thealgorithms.physics;
2+
3+
/**
4+
* Simulates a simple pendulum using the Runge-Kutta 4th order method.
5+
* The pendulum is modeled with the nonlinear differential equation.
6+
*/
7+
public final class SimplePendulumRK4 {
8+
9+
private SimplePendulumRK4() {
10+
throw new AssertionError("No instances.");
11+
}
12+
13+
private final double length; // meters
14+
private final double g; // acceleration due to gravity (m/s^2)
15+
16+
/**
17+
* Constructs a simple pendulum simulator.
18+
*
19+
* @param length the length of the pendulum in meters
20+
* @param g the acceleration due to gravity in m/s^2
21+
*/
22+
public SimplePendulumRK4(double length, double g) {
23+
if (length <= 0) {
24+
throw new IllegalArgumentException("Length must be positive");
25+
}
26+
if (g <= 0) {
27+
throw new IllegalArgumentException("Gravity must be positive");
28+
}
29+
this.length = length;
30+
this.g = g;
31+
}
32+
33+
/**
34+
* Computes the derivatives of the state vector.
35+
* State: [theta, omega] where theta is angle and omega is angular velocity.
36+
*
37+
* @param state the current state [theta, omega]
38+
* @return the derivatives [dtheta/dt, domega/dt]
39+
*/
40+
private double[] derivatives(double[] state) {
41+
double theta = state[0];
42+
double omega = state[1];
43+
double dtheta = omega;
44+
double domega = -(g / length) * Math.sin(theta);
45+
return new double[] {dtheta, domega};
46+
}
47+
48+
/**
49+
* Performs one time step using the RK4 method.
50+
*
51+
* @param state the current state [theta, omega]
52+
* @param dt the time step size
53+
* @return the new state after time dt
54+
*/
55+
public double[] stepRK4(double[] state, double dt) {
56+
if (state == null || state.length != 2) {
57+
throw new IllegalArgumentException("State must be array of length 2");
58+
}
59+
if (dt <= 0) {
60+
throw new IllegalArgumentException("Time step must be positive");
61+
}
62+
63+
double[] k1 = derivatives(state);
64+
double[] s2 = new double[] {state[0] + 0.5 * dt * k1[0], state[1] + 0.5 * dt * k1[1]};
65+
66+
double[] k2 = derivatives(s2);
67+
double[] s3 = new double[] {state[0] + 0.5 * dt * k2[0], state[1] + 0.5 * dt * k2[1]};
68+
69+
double[] k3 = derivatives(s3);
70+
double[] s4 = new double[] {state[0] + dt * k3[0], state[1] + dt * k3[1]};
71+
72+
double[] k4 = derivatives(s4);
73+
74+
double thetaNext = state[0] + (dt / 6.0) * (k1[0] + 2 * k2[0] + 2 * k3[0] + k4[0]);
75+
double omegaNext = state[1] + (dt / 6.0) * (k1[1] + 2 * k2[1] + 2 * k3[1] + k4[1]);
76+
77+
return new double[] {thetaNext, omegaNext};
78+
}
79+
80+
/**
81+
* Simulates the pendulum for a given duration.
82+
*
83+
* @param initialState the initial state [theta, omega]
84+
* @param dt the time step size
85+
* @param steps the number of steps to simulate
86+
* @return array of states at each step
87+
*/
88+
public double[][] simulate(double[] initialState, double dt, int steps) {
89+
double[][] trajectory = new double[steps + 1][2];
90+
trajectory[0] = initialState.clone();
91+
92+
double[] currentState = initialState.clone();
93+
for (int i = 1; i <= steps; i++) {
94+
currentState = stepRK4(currentState, dt);
95+
trajectory[i] = currentState.clone();
96+
}
97+
98+
return trajectory;
99+
}
100+
101+
/**
102+
* Calculates the total energy of the pendulum.
103+
* E = (1/2) * m * L^2 * omega^2 + m * g * L * (1 - cos(theta))
104+
* We use m = 1 for simplicity.
105+
*
106+
* @param state the current state [theta, omega]
107+
* @return the total energy
108+
*/
109+
public double calculateEnergy(double[] state) {
110+
double theta = state[0];
111+
double omega = state[1];
112+
double kineticEnergy = 0.5 * length * length * omega * omega;
113+
double potentialEnergy = g * length * (1 - Math.cos(theta));
114+
return kineticEnergy + potentialEnergy;
115+
}
116+
117+
public double getLength() {
118+
return length;
119+
}
120+
121+
public double getGravity() {
122+
return g;
123+
}
124+
}

0 commit comments

Comments
 (0)