Catalog
  1. 1. Purpose
  2. 2. Theory
  3. 3. The Results
  4. 4. The code
Python Plot (1) A small ball on a smooth semicircle (Animated)

All the python study code can be found in my own repository.

Purpose

In high school we learned how to use the energy conservation law to get the velocity for a ball slides down a 1/4 circle. However, the time needed in this process is still unknown. The period for a pendulum also uses a approximated expression. In this note, I will try to solve the time evolution for a ball slide down from a smooth semi-circle numerically via python. I will compare the oscillator approximation and accurate result in the same animated figure .

Theory

According to the energy conservation law

We need a linear differential form and the velocity $v$ should be used. The relation between the angle $\theta$ and the time is

We only consider the simple case that the small ball is initially rest, that is

The x and y coordinate can be expressed as

If $\theta$ is near 0,then $\sin(\theta)\approx \theta$,

This is the equation for a simple harmonic oscillator, when $\theta$ is small enough, this equation works well.

The Results

I used the function “odeint” to solve the differential equations numerically. I give two examples with the angle $\theta$ is very big or relatively small.

  • With the ball initially in the $\theta=\pi/6$, the evolution of oscillator approximations and accurate cases are different but not so big.

  • With the ball initially in the $\theta=\pi/2$, the oscillator approximation is not a good approximation now.

The code

1
# To plot a animated oscillator
2
# Use to spot, one spot is fixed at the origin and the
3
# other is moving
4
5
# import necessary library
6
import numpy as np
7
from matplotlib import pyplot as plt
8
from matplotlib import animation
9
from scipy.integrate import odeint
10
11
# Define the parameters
12
theta_0 = -np.pi/6
13
theta_0_degree = theta_0/np.pi*180
14
R = 1.0
15
g = 9.8
16
coe = 1
17
18
# Define the differential equations
19
20
21
def evo_accu(z, t, g, R, coe):
22
    v = z[0]
23
    theta = z[1]
24
    dtheta = v/R
25
    dv = -coe*g*np.sin(theta)
26
    dzdt = [dv, dtheta]
27
    return dzdt
28
29
30
def evo_osci(z, t, g, R, coe):
31
    v = z[0]
32
    theta = z[1]
33
    dtheta = v/R
34
    dv = -coe*g*theta
35
    dzdt = [dv, dtheta]
36
    return dzdt
37
38
39
t_max = 1.0/R*g*3
40
numt = 500
41
tmat = np.linspace(0, t_max, numt)
42
z0 = [0, theta_0]
43
theta_accu = odeint(evo_accu, z0, tmat, args=(g, R, coe))
44
theta_osci = odeint(evo_osci, z0, tmat, args=(g, R, coe))
45
# The coordinate expression
46
x_accu = R*np.sin(theta_accu[:, 1])
47
y_accu = -R*np.cos(theta_accu[:, 1])
48
49
x_osci = R*np.sin(theta_osci[:, 1])
50
y_osci = -R*np.cos(theta_osci[:, 1])
51
52
theta_accu_degree = theta_accu[:, 1]/np.pi*180
53
theta_osci_degree = theta_osci[:, 1]/np.pi*180
54
# Set the curve of the half circle
55
x_cir = np.linspace(-R, R, 100)
56
y_cir = -np.sqrt(R**2-x_cir**2)
57
# Count the loop
58
count_accu = np.zeros(numt)
59
count_osci = np.zeros(numt)
60
61
counter_accu = 0
62
counter_osci = 0
63
for i in range(numt-1):
64
    if(theta_accu[i, 1]*theta_accu[i+1, 1] < 0):
65
        counter_accu = counter_accu+1
66
        count_accu[i] = counter_accu
67
    else:
68
        count_accu[i] = counter_accu
69
70
    if(theta_osci[i, 1]*theta_osci[i+1, 1] < 0):
71
        counter_osci = counter_osci+1
72
        count_osci[i] = counter_osci
73
    else:
74
        count_osci[i] = counter_osci
75
        
76
# plot the two different evolution in the same figure
77
# First set up the figure, the axis, and the plot element we want to animate
78
79
fig = plt.figure()
80
ax = plt.axes(xlim=(-2, 2), ylim=(-2, 2), aspect='equal')
81
line_accu, = ax.plot([], [], 'r-o', label='Accurate')
82
line_osci, = ax.plot([], [], 'g-o', label='Oscillator Approx')
83
theta_accu_text = ax.text(0.02, 0.8, '', transform=ax.transAxes)
84
theta_osci_text = ax.text(0.02, 0.7, '', transform=ax.transAxes)
85
loop_accu_text = ax.text(0.02, 0.6, '', transform=ax.transAxes)
86
loop_osci_text = ax.text(0.02, 0.5, '', transform=ax.transAxes)
87
line_cir, = ax.plot(x_cir, y_cir, '--', label='Circle Half')
88
# initialization function: plot the background of each frame
89
90
91
def init():
92
    line_accu.set_data([], [])
93
    line_osci.set_data([], [])
94
    theta_accu_text.set_text('')
95
    theta_osci_text.set_text('')
96
    loop_accu_text.set_text('')
97
    loop_osci_text.set_text('')
98
    return line_accu, line_osci, theta_accu_text, theta_osci_text, loop_accu_text, loop_osci_text,
99
100
# animation function.  This is called sequentially
101
102
def animate(i):
103
    x_accu_mat = [0, x_accu[i]]
104
    y_accu_mat = [0, y_accu[i]]
105
    x_osci_mat = [0, x_osci[i]]
106
    y_osci_mat = [0, y_osci[i]]
107
    line_accu.set_data(x_accu_mat, y_accu_mat)
108
    line_osci.set_data(x_osci_mat, y_osci_mat)
109
    theta_accu_text.set_text(
110
        '$ \Theta_{accu} $ = %.2f $^o $' % theta_accu_degree[i])
111
    theta_osci_text.set_text(
112
        '$ \Theta_{osci} $ = %.2f $^o $' % theta_osci_degree[i])
113
    loop_accu_text.set_text('loop accu=%.1d' % count_accu[i])
114
    loop_osci_text.set_text('loop osci=%.1d' % count_osci[i])
115
    return line_accu, line_osci, theta_accu_text, theta_osci_text, loop_accu_text, loop_osci_text,
116
117
118
# call the animator.  blit=True means only re-draw the parts that have changed.
119
anim = animation.FuncAnimation(fig, animate, init_func=init,
120
                               frames=numt, interval=10, blit=False)
121
plt.legend(loc='best')
122
plt.xlabel('X Direction (m)')
123
plt.ylabel('Y Direction (m)')
124
ax.text(0.02, 0.02, 'g = %.2f m/s' % g, transform=ax.transAxes)
125
ax.text(0.32, 0.02, 'R = %.1f m' % R, transform=ax.transAxes)
126
ax.text(0.62, 0.02, '$\Theta_{0}$ = %.1d$^o $' %
127
        theta_0_degree, transform=ax.transAxes)
128
# plt.show()
129
anim.save("Ball_1.gif", writer='pillow')
Author: Knifelee
Link: https://knifelees3.github.io/2019/07/23/A_En_Python_PlotStudy1_SemiCircle_OscilatorApprox/
Copyright Notice: All articles in this blog are licensed under CC BY-NC-SA 4.0 unless stating additionally.
Donate
  • 微信
  • 支付寶

Comment