In [ ]:
# Analyzing the cycloid graph from "The Brachistochrone" research paper
File description: This ipython notebook attempts to graph and analyze a
simple cycloid curve, from a set of parametric polar
equations derived via the "Euler-Lagrange method" in
the paper "The Brachistochrone", which attempted to
prove, guided by the principle of Fermat's Theorem of
Least Time and building upon Newton-Galileo's solution to
Bernoulli's problem in Acta Eruditorum: given two points
A and B in the vertical plane, the curve traced out by a
point acted on only by gravity, which starts at A and reaches
B in the shortest time is a cycloidial path called a
Brachistochrone. Specifically, the paper shows that the
Brachistochrone curve, which is essentially an inverted
cycloid graph drawn by any circle rolling along an even plane,
can be graphed with parametric functions taking in polar parameters
r and theta. This notebook uses python to graph the parametric
functions to further analyze their meaning, as an extension
and continuation of past research.
In [19]:
# import numpy for math symbols and matplotlib for plotting data
import numpy as np
import matplotlib.pyplot as plt
In [20]:
# define range for theta [0, 8pi]
theta = np.linspace(0, 8*np.pi, 1000)
In [21]:
# define a radius for the cycloid-drawing circle
r = 50
In [22]:
# define parametric functions x and y
x = r * (theta - np.sin(theta))
y = r * (1 - np.cos(theta))
In [38]:
# Get subplot, and select axis
fig, ax = plt.subplots()
# Set grids to true, and set origin of axes
ax.grid(True, which='both')
ax.axhline(y=0, color='k')
ax.axvline(x=0, color='k')
Out[38]:
In [42]:
# Give graph a title
plt.title('Parametric Functions of a Cycloid Curve')
Out[42]:
In [44]:
# define range for x-axis [0, 1300]
left, right = plt.xlim()
plt.xlim(right=1300)
# define range for y-axis [0, 200]
bottom, top = plt.ylim()
plt.ylim(top=200)
# plot the cycloid path in black
plt.plot(x, y, color='black')
Out[44]:
In [ ]:
As we can see, the graph describes the rolling motion of a circle of
radius 50 along an even plane, successfully drawing out multiple cycloid
curves within the range of [0, 4*pi].
Overall, this gives us a visual of a cycloid curve, but since a Brachistochrone
is an inverted cycloid, we must flip the graph by reflecting it upon
the y-axis. To do this, we must change the range of the y-axis, and
apply the reflective transformation to the y parameter of the
parametric equation upon calling the plot function.
In [46]:
# define range for y-axis [-200, 100]
bottom, top = plt.ylim()
plt.ylim(bottom=-100, top=100)
# plot the cycloid path in black
plt.plot(x, -y, color='black')
Out[46]:
In [ ]:
Now, the graph correctly demonstrates the brachistochrone curve, the path
of shortest time between two points A and B on a level plane.
In [ ]:
Let's say we want to determine if the brachistochrone curve is indeed
faster. To do this, we will test a range of paths drawn on a coordinate
axis starting from the origin. Just for simplicity's sake, we will use
the coordinates of the edges of a circle graph, with radius 1 around the
origin, and limit the range of coordinates to Quadrant IV (+, -). This is
because intuitively, we are drawing paths which are acted on by gravity
in a positive direction and neglecting paths which simply allow an
object to free fall once put on the cycloid path. The purpose of this
experiment is ultimately to determine if a cycloid curve can best a
straight line in terms of the time needed to travel down its length.
The equations we will use are given by the research paper, which has
already proven their validity.
In [ ]:
For example, take a point (1, -1) and draw a path to it from the origin.
This is the straight line path we will be comparing our findings to:
In [102]:
# define coordinates
x = np.arange(0, 1, 0.01)
y = -x
# Get subplot, and select axis
fig, ax = plt.subplots()
# Set grids to true, and set origin of axes
ax.grid(True, which='both')
ax.axhline(y=0, color='k')
ax.axvline(x=0, color='k')
ax.plot(x, y, color="black")
Out[102]:
In [ ]:
Then, using the parametric functions we derived from
our research, we can solve for the unknown parameter r
and determine the theta graph to draw the cycloid.
In [163]:
# Consider the parameteric functions
X: 1 = r(THETA - sin(THETA))
Y: -1 = r(1 - cos(THETA))
Y: r = -1/(1-cos(THETA))
# Hence, substituting r into X,
1 = -(THETA - sin(THETA)) / (1 - cos(THETA))
# Graphing this function:
# define function, set x to THETA to determine x-intercept
x = np.linspace(-5, 5, 100)
y = -((x - np.sin(x))/(1 - np.cos(x))) - 1
# Ignore divide by 0
np.seterr(divide='ignore', invalid='ignore')
# Get subplot, and select axis
fig, ax = plt.subplots()
# Set grids to true, and set origin of axes
ax.grid(True, which='both')
ax.axhline(y=0, color='k')
ax.axvline(x=0, color='k')
ax.plot(x, y, color="black")
Out[163]:
In [ ]:
From this graph, it is clear that the x-intercept occurs at approximately
x = -2.4. Translating this to radians, THETA = -2.4 Rad. Now, we can solve
for the radius r.
In [ ]:
# Recall from before,
Y: r = -1/(1-cos(THETA))
# Hence,
r = -1/(1 - cos(-2.4))
r = -0.58 units approximately
In [ ]:
Now, we can graph the cycloid
In [125]:
# define range for theta [0, 8pi]
theta = np.linspace(0, -2.4, 100)
r = -0.58
# using r = -0.58, define parametric functions x and y
x = r * (theta - np.sin(theta))
y = r * (1 - np.cos(theta))
# Ignore divide by 0
np.seterr(divide='ignore', invalid='ignore')
# Get subplot, and select axis
fig, ax = plt.subplots()
# Set grids to true, and set origin of axes
ax.grid(True, which='both')
ax.axhline(y=0, color='k')
ax.axvline(x=0, color='k')
ax.plot(x, y, color="black")
Out[125]:
In [ ]:
Let's compare the graphs side-by-side
In [126]:
# define coordinates
x = np.arange(0, 1, 0.01)
y = -x
# Get subplot, and select axis
fig, ax = plt.subplots()
# Set grids to true, and set origin of axes
ax.grid(True, which='both')
ax.axhline(y=0, color='k')
ax.axvline(x=0, color='k')
ax.plot(x, y, color="black")
# define range for theta [0, 8pi]
theta = np.linspace(0, -2.4, 100)
r = -0.58
# using r = -0.58, define parametric functions x and y
x = r * (theta - np.sin(theta))
y = r * (1 - np.cos(theta))
# Ignore divide by 0
np.seterr(divide='ignore', invalid='ignore')
ax.plot(x, y, color="black")
Out[126]:
In [ ]:
We now have an accurate visual depiction of the
fastest path between points (0, 0) and point(1, -1). Not
convinced? There is still more to analyze. Using our equation
for least time traveled along this Brachistochrone path,
we can populate a table comparing the times taken by a
straight path and a cycloid path.
In [ ]:
To calculate the time taken by the straight path, we must
use Galileo's equation, thanks to him, we can use the following
equation to determine the time taken for an object at rest
on top of an inclined plane to reach the bottom.
t = square-root[d * (2*l) / (g*h)],
Where t is the time in seconds, d is the length of the ramp,
g is gravitational acceleration, and l and h are the width
and height of the ramp, respectively.
From our own research, we have determined that the time taken
for an object to travel down a cycloid ramp is:
t = square-root[r/g] * THETA,
Where r and theta are the parametric function parameters.
In [ ]:
We will now use pandas, a dataframe library, to table all
of our data comparisons. Although we must first calculate
the necessary parameters to determine the time taken, for
this, we will use the following set of coordinates in
Quadrant IV of the coordinate axis:
{(x,y) | x^2 + y^2 = 1, x:[0, 1], y:[0, -1]}
This simplifies the data collection process. For straight
path, we have set d = 1, l = x, h = -y. Thus,
t = square-root[(2*x) / (g*(-y))]
For the cycloid path, we have to calculate r and theta for
each coordinate. To discretize the graph, we will use
coordinates rounded to 4 decimal places
In [161]:
# begin by importing pandas and math
import math as math
import pandas as pd
# Gravitational acceleration (9.81m/s2)
g = 9.81
# Get set of test-coordinates
pi = np.pi
# function to return n (x,y) coordinates lining circumference
def PointsInCircum(x_or_y, r, n):
if x_or_y == 1:
return [round(math.cos(-pi/(2*n)*x)*r, 4) for x in range(1,n)]
else:
return [round(math.sin(-pi/(2*n)*x)*r, 4) for x in range(1,n)]
x_vals = PointsInCircum(1, 1, 50)
y_vals = PointsInCircum(0, 1, 50)
print("X: ")
print(x_vals)
print("Y: ")
print(y_vals)
In [316]:
# Determine the time values for straight graph
t_strt = [round(math.sqrt((2*x_vals[index])/(g*(-y_vals[index]))), 4) for index in range(0, 49)]
In [301]:
# Determine time values for cycloid graph: to solve for r, we can use the following:
y = r * (1 - cos(THETA))
r = y/(1 - cos(THETA))
# Array for storing theta values
t_theta = []
# calculate the y-coordinate of the graph
def calculate_y(x_cord, y_cord):
return y_cord*((x - np.sin(x))/(1 - np.cos(x))) - x_cord
# find x-intercept of graph via approximation
myrange = np.arange(-100, 100, 0.01)
index = 0
# Iterate through values and find closest x-intercept
while index<49:
for x in myrange:
y = calculate_y(x_vals[index], y_vals[index])
if 0.00<= y <=0.01:
t_theta.append(round(x,3))
break
index = index + 1
print(t_theta)
In [310]:
# once all the theta values are calculated, we must calculate the r values
t_rvals = []
index = 0
while index<49:
rvals = -y_vals[index]/(1-np.cos(t_theta[index]))
t_rvals.append(round(rvals,3))
index = index + 1
print (t_rvals)
In [324]:
# Hence, the cycloid graph time values are:
t_cycl = []
g = 9.81
index = 0
while index<49:
time = -math.sqrt(t_rvals[index]/g)*t_theta[index]
t_cycl.append(round(time,4))
index = index + 1
print (t_cycl)
In [325]:
# Comparing the data side-by-side,
print("Straight Path (s): ")
print(t_strt)
print()
print("Cycloid Path (s): ")
print(t_cycl)
In [326]:
# Start by defining the index for tabling
indicies = []
index = 1
while index<50:
indicies.append(index)
index += 1
print(indicies)
In [329]:
# Create data for tabling
data = {'Point number': indicies,
'Straight Path (s)': t_strt,
'Cycloid Path (s)': t_cycl}
df = pd.DataFrame(data)
df
Out[329]:
In [ ]:
As we can see, for the first 20 data points, the brachistochrone path
is indeed faster, but after a while, the data points plateau, and the
straight path decreases in time marginally. This is likely because the
path has become near perfectly vertical meaning there is no frictino
to prevent the full force of gravitational acceleration, whereas for the
cycloid curve there may yet exist a path. Let's graph this and get a
visual feel for this trend!
In [346]:
# plot both paths in different colors
plt.figure(figsize=(8, 6), dpi=80)
plt.title("Time Difference of Straight and Cycloid Paths")
plt.xlabel("Point Number")
plt.ylabel("Time taken(s)")
# Set intervals for axes
plt.xticks(np.arange(0, 51, 3))
plt.yticks(np.arange(0.0, 3.0, 0.25))
x = df.iloc[:,0]
y = df.iloc[:,1]
plt.plot(x, y, color='blue')
x = df.iloc[:,0]
y = df.iloc[:,2]
plt.plot(x, y, color='red')
Out[346]:
In [ ]:
And that concludes our analysis! It has been a very enlightening research
and analysis process being able to analyze real-time data with python
frameworks. We have not only concluded that cycloid paths are indeed faster
most of the time, but also had the opportunity to explore many different
ways of visualizing data. Overall, this has been a very enjoyable research
project for me, something I was not able to finish back in high school,
and I hope to continue in the future!