|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 |
#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; } |
Explanation
- Constants:
PI,G, andDTare constants representing the value of π, the gravitational constant, and the time step for the simulation, respectively.
- PendulumState Structure:
- This structure holds the state of the double pendulum, including the angles (
theta1,theta2) and angular velocities (omega1,omega2) of the two pendulums.
- This structure holds the state of the double pendulum, including the angles (
- computeDerivatives Function:
- This function calculates the derivatives of the angles and angular velocities based on the current state using the equations of motion for the double pendulum. The calculations take into account the masses (
m1,m2) and lengths (l1,l2) of the pendulum arms.
- This function calculates the derivatives of the angles and angular velocities based on the current state using the equations of motion for the double pendulum. The calculations take into account the masses (
- updateState Function:
- This function updates the state of the pendulum by integrating the derivatives using Euler’s method. This method is a simple numerical technique for solving differential equations.
- Main Function:
- The main function initializes the state with specific angles and angular velocities.
- The simulation runs for 1000 iterations, updating and printing the state of the pendulum at each time step.
Usage
- Initial Angles and Velocities: The angles
theta1andtheta2are initialized to π/2 (90°) and π/4 (45°), respectively, representing the initial displacement of the pendulum arms. - Masses and Lengths: Both pendulums have equal masses and lengths, but these can be varied to observe different behaviors.
- Output: The program outputs the angles of the two pendulums at each time step, providing insight into the complex and chaotic motion of the double pendulum.