Consider the problem of finding a curve that connects two points A and B, such that an object rolling down the curve from A to B takes the least amount of time to make the trip
This is called the Brachistochrone problem, a classic problem in the Calculus of Variations.
This notebook contains functions which will let you set up your own version of the problem, sketch out your own solution using your mouse, then generate a 3D-printable ramp to roll marbles down in order to test your design.
Source code for all func
The functions of this notebook require:
The python libraries can be installed by running pip install -r requirements.txt within the directory that contains this notebook file.
The computer vision library OpenCV can be installed by following their installation instructions. OpenCV can be tricky to install. I'd recommend using a system package manager like apt-get or homebrew if possible. On Windows, the easiest way to get started with OpenCV and Python is to use either Python(x,y) or the Anaconda Python Distribution
The CAD program OpenSCAD can be downloaded from here. This notebook involves calling OpenSCAD from the command line. To make this work, the path to the OpenSCAD binary will need to be given in the next cell.
First, let's make all necessary imports and set up the conditions (the location of the start and end points A and B) for the problem. Feel free to modify these points if you'd like.
The path to the OpenSCAD (including the name of the binary) needs to be given below. The default location given is the path on my Mac, you'll need to change this if necessary. If OpenSCAD is already in your system's path, the binary name should be all that is needed here.
In [ ]:
from lib.ramp import ramp, render_stl
from lib.interface import get_sketch, smooth_plot, show_stl, Animate_frame
from lib.fit import fit_multiple, rescale
from lib.cycloid import cycloid, total_time
from IPython.html.widgets import interactive, FloatSliderWidget
from IPython.display import display
from IPython.display import HTML, Image
import os
openscad_path = "/Applications/OpenSCAD.app/Contents/MacOS/OpenSCAD"
# tests the given path
if not os.path.exists(openscad_path):
raise Exception("OpenSCAD path is incorrect!")
A = [0, 95.0] # Start point
B = [210, 55.0] # End point
n_points = 250 # resolution
Run the next cell, and draw a curve with your mouse between the two dots shown in the window.
(Optional) Use the widget below to smooth out your drawn curve after it has been collected. If sigma = 0 (default), no smoothing will be performed.
In [ ]:
x,y = get_sketch(A, B, n_points)
y_smooth = y
def f(sigma):
global y_smooth
y_smooth = smooth_plot(x, -y, sigma, A, B)
sigma_slider = FloatSliderWidget(min=0, max=15, step=0.3, value=0)
w=interactive(f,sigma=sigma_slider)
display(w)
Now collect some more curves that connect A to B to compare against, plot them, and compute a prediction of descent time for each curve.
Move the slider to simulate a small object rolling down each curve over time.
In [ ]:
# Generate an inverted cycloid between A and B with a given number of sample points
x, y0 = cycloid(A, B, n=y_smooth.size)
# Generate a parabola
y1 = (x- B[0]*0.75)**2
# Generate a wavey curve
y2 = 0.5*x + 2.5*np.sin(0.1*x + 1*np.pi/2)
# Now fit everything together between A and B
y0, y_smooth, y1, y2, y3 = fit_multiple([y0, y_smooth,y1, y2, x], A, B)
t_drawn = total_time(x, y_smooth)
t_wavey = total_time(x, y2)
t_parabola = total_time(x, y1)
t_straight = total_time(x, y3)
t_cycloid = total_time(x, y0)
t_max = max([t_drawn, t_wavey, t_parabola, t_straight, t_cycloid])
# plots
names = ["Drawn","Wavey","Parabola","Straight","Inv. Cycloid"]
Ys = [y_smooth, y2, y1, y3, y0]
def animate(t):
Animate_frame(Ys, x, t, names)
text(A[0]+5, A[1]+10, "A", fontsize=20)
text(B[0]+5, B[1]+10, "B", fontsize=20)
t_slider = FloatSliderWidget(min=0, max=t_max+0.001, step=t_max/100, value=0)
w=interactive(animate, t=t_slider)
display(w)
In [ ]:
print
print "Estimated total descent times:"
print 30*"-"
print "Drawn : ", t_drawn
print "Wavey : ", t_wavey
print "Parabola :", t_parabola
print "Straight :", t_straight
print "Cycloid : ", t_cycloid
print
Finally, generate a 3D STL file of marble ramps based on the above curves
In [ ]:
# -- this will generate the raw OpenSCAD code that will create the marble ramps
# unless re-scaling in the 3D slicing software before printing, all units are in mm
name = "my_ramps" # name to use for the OpenSCAD modules and stl file
thickness = 8 # thickness of the base
length = B[0] - A[0] # length of the ramps
n = 250 # resolution of the ramp surfaces
back_height = 10 # height of the backing on the starting (A) side
back_depth = 15 # depth of the backing on the starting (A) side
channel_width = 20 # width of the marble channels
channel_depth = 10 # depth of the marble channels
wall_width = 1.5 # thickness of the channel walls
endstop_thickness = 5 # thickness of the end stop block
endstop_height = 30 # height of the end stop block
code = ramp([y_smooth,y2, y1, y3, y0],
thickness=thickness,
length=length,
start_backdepth = back_depth,
endstop_thickness=endstop_thickness,
endstop_height=endstop_height,
channel_width = channel_width,
wall_width = wall_width,
channel_depth = channel_depth,
n=n,
start_backheight=back_height,
name="ramp")
# Now, save the scad file, and export to .stl from the command line (may take a minute)
render_stl(openscad_path, code, name)
# Show the STL file below
show_stl(name)
The last steps above (STL generation and preview) may take a minute or so to complete.
The ramps seen are, starting on the left (facing the object from the front end/lower part):
If you print out these ramps and test them with marbles, you'll find that a marble rolling down the last curve (the inverted cycloid) will tend to reach the end point the fastest.
When prepping the STL file for printing, you can scale the object to fit the volume of your printer as long as you scale each axis uniformly (together). If the length vs height ratio of the ramps is modified, then the descent properties of the curve may be very different than what was computed in this notebook.
Neglecting friction and drag effects, the total time for an object to roll from A$=\left(x_0, y_0\right)$ to B$=\left(x_1, y_1\right)$ along a curve $y\left(x\right)$ in a uniform graviational field can be computed as
$$T = \int_{x_0}^{x_1} \left( \left(y^\prime\left(x\right) \right)^2 + 1 \right)^{\frac{1}{2}} \left(2 g y\left(x\right) \right)^{-\frac{1}{2}} dx$$where $y'\left(x\right) = \frac{dy}{dx}$ and $g$ is a gravitational acceleration constant. For simplicity, let $x_0 = y_0 = 0$ and let $x_1$ and $y_1$ shift appropriately to $x_n,y_n$ in order to represent the same total length and change of height. So the Brachistochrone problem then has the form:
$$\min_{y\left(x\right) \in C^2[0,x_1]} \int_{0}^{x_1} \left( \left(y' \right)^2 + 1 \right)^{\frac{1}{2}} \left(2 g y \right)^{-\frac{1}{2}} dx$$$$ =\sqrt{2g }^{-1} \min_{y\left(x\right) \in C^2[0,x_1]} \int_{0}^{x_1} F\left(y'\left(x\right), y\left(x\right), x\right) dx $$subject to the boundary conditions $y\left(0 \right ) = 0$ and $y\left(x_n \right ) = y_n$, where $$F\left(y'\left(x\right), y\left(x\right), x\right) = \sqrt{\frac{ y'\left(x \right)^2 + 1 }{y}}$$
The Calculus of Variations is a subject concerned with solving optimization problems of this type, where the values of an entire function between two points collectively serve as the parameter to be varied.
As it turns out, a first-order necessary condition for a local minimizer $y$ is that it satisifies the associated Euler-Lagrange equation:
$$\frac{\partial F}{\partial y} - \frac{d}{dx} \frac{\partial F}{\partial y'} = 0$$and if $F$ has no explicit dependance on $x$, then the E-L equation takes the form
$$F - y'\frac{\partial F}{\partial y'} = C $$for a constant $C$. For the Brachistochrone problem
$$ y'\frac{\partial F}{\partial y'} = y'^2 \left(y\left(1+y'^2\right) \right)^{-\frac{1}{2}} $$which leads to
$$ \sqrt{\frac{ y'^2 + 1 }{y}} - y'^2 \left(y\left(1+y'^2\right) \right)^{-\frac{1}{2}} = C $$$$\Rightarrow \frac{1}{\sqrt{y \left(1 + y'^2 \right)}} = C$$$$\Rightarrow y \left(1 + y'^2 \right) = \frac{1}{C^2}$$which is solved parametrically by the equations
$$x = \frac{1}{2C^2}\left(\theta - \sin\theta \right) = a\left(\theta - \sin\theta \right) $$$$y = \frac{1}{2C^2}\left(1 - \cos\theta \right) = a\left(1 - \cos\theta \right) $$for a constant $a$, which are recognized as the equations of a cycloid.
Shifting the problem back to A$=\left(x_0, y_0\right)$ to B$=\left(x_1, y_1\right)$, the equations are
$$x = x_0 + a\left(\theta - \sin\theta \right) $$$$y = y_0 + a\left(1 - \cos\theta \right) $$The constant $a$ can be computed based on known information. First, we know that at the end point B
$$\frac{y_0 - y_1}{x_1 - x_0} =\frac{\left(1 - \cos\theta \right)}{\left(\theta - \sin\theta \right)}$$which we can solve this numerically for $\theta \in [0, 2\pi]$. Then we can plugin $\theta$ into
$$x = x_0 + a\left(\theta - \sin\theta \right) $$to determine $a$.
The cycloid function in the file lib/cycloid.py computes a Cycloid between the two points given in the first cell of this notebook in exactly this manner.
In [ ]: