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