#include <iostream>
#include <cmath>
#include <vector>
// Constants
const double PI = 3.141592653589793;
const double G = 9.81; // Acceleration due to gravity (m/s^2)
const double DT = 0.01; // Time step for the simulation (s)
// Structure to hold the state of the pendulum
struct PendulumState {
double theta1; // Angle of the first pendulum (radians)
double theta2; // Angle of the second pendulum (radians)
double omega1; // Angular velocity of the first pendulum (rad/s)
double omega2; // Angular velocity of the second pendulum (rad/s)
};
// Function to compute the derivatives for the double pendulum
PendulumState computeDerivatives(const PendulumState& state, double m1, double m2, double l1, double l2) {
PendulumState derivatives;
double delta = state.theta2 - state.theta1;
double den1 = (m1 + m2) * l1 - m2 * l1 * std::cos(delta) * std::cos(delta);
double den2 = (l2 / l1) * den1;
derivatives.omega1 = ((m2 * l1 * state.omega1 * state.omega1 * std::sin(delta) * std::cos(delta) +
m2 * G * std::sin(state.theta2) * std::cos(delta) +
m2 * l2 * state.omega2 * state.omega2 * std::sin(delta) -
(m1 + m2) * G * std::sin(state.theta1)) / den1);
derivatives.omega2 = ((-m2 * l2 * state.omega2 * state.omega2 * std::sin(delta) * std::cos(delta) +
(m1 + m2) * G * std::sin(state.theta1) * std::cos(delta) -
(m1 + m2) * l1 * state.omega1 * state.omega1 * std::sin(delta) -
(m1 + m2) * G * std::sin(state.theta2)) / den2);
derivatives.theta1 = state.omega1;
derivatives.theta2 = state.omega2;
return derivatives;
}
// Function to update the state using Euler's method
void updateState(PendulumState& state, const PendulumState& derivatives) {
state.theta1 += derivatives.theta1 * DT;
state.theta2 += derivatives.theta2 * DT;
state.omega1 += derivatives.omega1 * DT;
state.omega2 += derivatives.omega2 * DT;
}
int main() {
// Initial conditions
PendulumState state;
state.theta1 = PI / 2; // 90 degrees
state.theta2 = PI / 4; // 45 degrees
state.omega1 = 0.0;
state.omega2 = 0.0;
double m1 = 1.0; // Mass of the first pendulum (kg)
double m2 = 1.0; // Mass of the second pendulum (kg)
double l1 = 1.0; // Length of the first pendulum (m)
double l2 = 1.0; // Length of the second pendulum (m)
for (int i = 0; i < 1000; ++i) {
PendulumState derivatives = computeDerivatives(state, m1, m2, l1, l2);
updateState(state, derivatives);
std::cout << "Time: " << i * DT << " s, "
<< "Theta1: " << state.theta1 << " rad, "
<< "Theta2: " << state.theta2 << " rad\n";
}
return 0;
}
