fork download
copy 
  1. # -*- coding: utf-8 -*-
  2. """
  3. Created on Thu Jun  5 10:48:06 2025
  4.  
  5. @author: eddyb
  6. """
  7. import numpy as np
  8. import matplotlib.pyplot as plt
  9. from math import sin, cos
  10.  
  11. # Physical parameters
  12. g = 9.81       # gravity (m/s²)
  13. L = 1.0        # pendulum length (m)
  14. theta0 = np.pi/4  # initial angle (radians)
  15. omega0 = 0.0      # initial angular velocity (rad/s)
  16.  
  17. # Simulation parameters
  18. t_max = 10.0    # simulation time (s)
  19. dt = 0.01       # time step (s)
  20.  
  21. # Time array
  22. t = np.arange(0, t_max, dt)
  23.  
  24. # Initialize arrays
  25. theta = np.zeros_like(t)
  26. omega = np.zeros_like(t)
  27. theta[0] = theta0
  28. omega[0] = omega0
  29.  
  30. # Energy arrays for verification
  31. E_pot = np.zeros_like(t)  # potential energy
  32. E_kin = np.zeros_like(t)  # kinetic energy
  33. E_total = np.zeros_like(t)  # total energy
  34.  
  35. # Verlet integration (better for oscillatory systems)
  36. for i in range(1, len(t)):
  37.     # Calculate acceleration
  38.     alpha = -(g/L) * sin(theta[i-1])
  39.  
  40.     # Verlet velocity update
  41.     omega[i] = omega[i-1] + alpha * dt
  42.     theta[i] = theta[i-1] + omega[i] * dt  # Note we use omega[i] here
  43.  
  44.     # Calculate energies for verification
  45.     E_pot[i] = g * L * (1 - cos(theta[i]))  # Potential energy
  46.     E_kin[i] = 0.5 * (L*omega[i])**2        # Kinetic energy
  47.     E_total[i] = E_pot[i] + E_kin[i]        # Total mechanical energy
  48.  
  49. # Convert to Cartesian coordinates
  50. x = L * np.sin(theta)
  51. y = -L * np.cos(theta)
  52.  
  53. # Create plots
  54. plt.figure(figsize=(14, 10))
  55.  
  56. # Angle vs Time
  57. plt.subplot(3, 2, 1)
  58. plt.plot(t, theta)
  59. plt.xlabel('Time (s)')
  60. plt.ylabel('Angle (rad)')
  61. plt.title('Angle vs Time')
  62. plt.grid()
  63.  
  64. # Angular Velocity vs Time
  65. plt.subplot(3, 2, 2)
  66. plt.plot(t, omega)
  67. plt.xlabel('Time (s)')
  68. plt.ylabel('Angular Velocity (rad/s)')
  69. plt.title('Angular Velocity vs Time')
  70. plt.grid()
  71.  
  72. # Phase Space
  73. plt.subplot(3, 2, 3)
  74. plt.plot(theta, omega)
  75. plt.xlabel('Angle (rad)')
  76. plt.ylabel('Angular Velocity (rad/s)')
  77. plt.title('Phase Space')
  78. plt.grid()
  79.  
  80. # Trajectory
  81. plt.subplot(3, 2, 4)
  82. plt.plot(x, y)
  83. plt.xlabel('x position (m)')
  84. plt.ylabel('y position (m)')
  85. plt.title('Pendulum Trajectory')
  86. plt.grid()
  87. plt.axis('equal')
  88.  
  89. # Energy plots
  90. plt.subplot(3, 2, 5)
  91. plt.plot(t, E_pot, label='Potential')
  92. plt.plot(t, E_kin, label='Kinetic')
  93. plt.xlabel('Time (s)')
  94. plt.ylabel('Energy (J)')
  95. plt.title('Potential and Kinetic Energy')
  96. plt.legend()
  97. plt.grid()
  98.  
  99. plt.subplot(3, 2, 6)
  100. plt.plot(t, E_total)
  101. plt.xlabel('Time (s)')
  102. plt.ylabel('Energy (J)')
  103. plt.title('Total Mechanical Energy (should be constant)')
  104. plt.grid()
  105.  
  106. plt.tight_layout()
  107. plt.show()
Success #stdin #stdout 1s 66752KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
Standard output is empty
https://ideone.com/noTS8K
language:
Python (cpython 2.7.16)
created:
3 days ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!