Modelling and Simulation in Python: Visualizing Natural Phenomena

5. Simulating a Swinging Pendulum

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\).

\(sin(\theta)\) is non-linear and is normally solved using a special case of the Taylor series, the Maclaurin series expansion (Knight & Adams, 1975) (Lam & Leary, 1979).

\[\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}\tag{3} \frac{d^2x_0}{dt^2} + \frac{g}{L}x_0 = 0 \end{align}\]

\[\begin{align*} \text{Period of oscillaiton, } T &= \frac{2\pi}{\omega} \\ \omega &= \sqrt{\frac{g}{L}} \\ T &= 2\pi \sqrt{\frac{L}{g}} \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

Initialize variables.

 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)

Prepare an empty plot.

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))

Animate your plot.

52anim = FuncAnimation(fig, 
53                    animate, 
54                    frames=len(angle), 
55                    interval=2,
56                    repeat=True,
57                    )
58plt.show()

Add Comment

Comments

No comments yet. Be the first!