Let the initial angle \(\theta_0 = \theta_{old}\) and initial velocity \(= v_{old}\). When the bob is released, \(\theta_{old}\) will reduce to \(\theta_{new}\) and \(v_{old}\) will change to \(v_{new}\), and so on.
\[\begin{align*}
W_{bob} &= mg \\
sin(\theta_{old}) &= \frac{\text{opp-side}}{mg}\\
\text{opp-side} &= mg \times sin(\theta_{old}) \\
\text{opp-side} &= -mg \times sin(\theta_{old}) \\
\end{align*}\]
The component \(mg \times sin(\theta)\) causes oscillation.
\[\begin{align*}
\text{adj-side} &= mg \times cos(\theta_{old}) \\
F &= ma \\
ma &= -mg \times sin(\theta_{old}) \\
m\frac{d^2x_0}{dt^2} &= -mg \times sin(\theta_{old}) \\
\frac{d^2x_0}{dt^2} &= -g \times sin(\theta_{old})
\end{align*}\]
The equation of motion here is independent of mass, \(m\).
\[\sum_{n=0}^{\infty}\frac{(-1)^n(\theta^{2n+1})}{(2n+1)!} \\
sin (\theta) = \theta - \frac{\theta^3}{3!} + \frac{\theta^5}{5!} - \cdots\]
For small \(\theta\), \(\approx \frac{\pi}{12} \text{ radians}\), \(sin(\theta)\) can then be approximated to \(sin(\theta) \approx \theta\), resulting to \(\approx 1\% \text{ error}\).
However, considering a very small \(\theta\), we can approximate \(sin(\theta)\).
\[\begin{align*}
sin(\theta_{old}) \approx \frac{x_0}{L}
\end{align*}\]
\[\begin{align*}
x_0 &= L \times sin(\theta_{old}) \\
y_0 &= -L \times cos(\theta_{old}) \\
acc &= -g \times sin(\theta_{old}) \\
v_{new} &= v_{old} - acc \times dt \\
\theta_{new} &= \theta_{old} - \theta_{change} \\
\theta &= \frac{\text{arc length}} {\text{radius}} \\
\theta_{change} &= \frac{\text{arc length}} {L} \\
\theta_{change} &= \frac{v_{old} \times dt} {L}
\end{align*}\]
The following code simulates the motion of a swinging pendulum using numerical integration, based on the principles of Newtonian mechanics. The code defines the initial parameters of the pendulum, including its mass, length, initial angle, and velocity. The loop iterates over a range of time steps, calculating the position and velocity of the pendulum at each step. The plot displays the pendulum’s motion using matplotlib, with a red circle representing the bob and a gray line indicating the rod. The simulation could be improved by using higher-order integration schemes, adding nonlinear damping terms, simulating multiple pendulums, and allowing user input. The code is useful for mechatronics engineers, who can use it to design and optimize control systems, simulate the motion of robotic arms, test and calibrate sensors, and harvest energy from mechanical systems. Check this code on GitHub and you are welcome to improve it.
Here is a swinging single pendulum simulation.
Import the required libraries.
1import numpy as np
2import matplotlib.pyplot as plt
3from matplotlib.animation import FuncAnimation
|
4m = 1 # kg (bob mass)
5L = 1 # m (pendulum lenth)
6g = 9.81 # m/s/s
7angle_0 = np.radians(60) # deg to rad (initial angular displacement)
8v_0 = 0 # m/s (tangential vel)
9step = 0 # starting simulation time step
10dt = 0.01 # time step for integration calculation
11angle = [angle_0] # prepare array for angle data
12v = [v_0] # prepare array for vel data
13angle_old = angle_0
|
Solve equations involved and store the data in an array.
14T = 2 * np.pi * np.sqrt(L / g)
15
16while 1==1:
17 step += 1
18 t = step * dt
19
20 angle_old, v_old = angle[-1], v[-1]
21 arc_lenth = v_old * dt
22 angle_change = arc_lenth/L
23 angle_new = angle_old - angle_change
24
25 acc = -g * np.sin(angle_old)
26 v_new = v_old - acc * dt
27
28 if t > T:
29 # end
30 break
31
32 angle.append(angle_new)
33 v.append(v_new)
|
34x0 = L * np.sin(angle_old)
35y0 = -L * np.cos(angle_old)
36
37fig, ax = plt.subplots()
38ax.set_xlim(-L*1.15, L*1.15)
39ax.set_ylim(-L*1.15, L*1.15)
40ax.set_xlabel('X')
41ax.set_ylabel('Y')
42ax.set_aspect('equal'); ax.text(0.95, 0.05, 'SiliconWit.com', transform=ax.transAxes, fontsize=10, color='gray', alpha=0.5, ha='right', va='bottom')
43graph = ax.plot([0, x0], [0, y0], lw=2, color='grey')[0]
44
45bob = ax.add_patch(plt.Circle((x0,y0), 0.08,
46 color='red', zorder=2))
|
Make a fuction to update the plot.
48def animate(frame_i):
49 x, y = L * np.sin(angle[frame_i]), -L * np.cos(angle[frame_i])
50 graph.set_data([0, x], [0, y])
51 bob.set_center((x, y))
|
52anim = FuncAnimation(fig,
53 animate,
54 frames=len(angle),
55 interval=2,
56 repeat=True,
57 )
58plt.show()
|
|
