The Brachistochrone

Consider two points located near the surface of the earth--some point A and another point B--such that neither point is directly above the other.  With respect to Newtonian physics, the shortest path from A to B is the line segment that connects the two points.  Now, suppose that we connect A and B with a frictionless wire, and then allow a bead to slide along the wire from A to B, moved along only by the pull of gravity.  For the bead to be able to move from A to B, we would need the point A to be as far as or even further from the surface of the earth than B.  So, let us assume this condition.  With this restriction we should be able to bend the wire in some manner so as to enable the bead to slide from A to B.  We now wish to consider how the wire should be bent in order to create a path for which the time of travel from A to B is minimized.

Note that it is tempting to assume that this path of minimum time is the same as the path of minimum length--the line connecting the points A and B.  However, as we shall see, this is not correct.  The solution is a curve, known as the brachistochrone.  The actual means of deriving the solution of this problem is beyond our grasp, and so we shall approach it by constructing a structure that approximates the solution.

The solution to our problem is some curve that passes from the point A to the point B, and since A and B can both be placed in a plane that is "vertical" with respect to the surface of the earth, we can argue that the solution also lies in this plane since the pull of gravity is vertical and not horizontal.  We can designate a set of coordinate axes for this plane, taking its intersection with the (immediate) earth's surface as the x-axis, and the vertical line through A as the y-axis.  If (0, y_A) represents the coordinates of A and (x_B, y_B) represents the coordinates of B, then the statement of our problem implies that we can designate x_B > 0 and 0    y_B    y_A.  Note that we would expect that the solution curve would move horizontally away from A, towards B, without moving horizontally back towards A again, so we shall assume this.  If we further assume that the solution does not consist of any vertical "drops" in position, then we can describe the solution as a function in the variable x, y = f(x), having domain [0, x_B].
We shall now approximate this function with a collection of line segments, approximate the time of travel along each of these line segments, and then adjust the line segments so as to minimize the total approximate time of travel along them.  To insure a smooth transition of motion, we shall assume that these line segments are connected by means of a curved joint instead of the sharp corner where they intersect.

Let A = (x₀, y₀), (x₁, y₁), ..., (x_n, y_n) = B be n + 1 points that lie along the solution curve y = f(x), with 0 = x₀ < x₁ < ··· < x_n = x_B and 0    y_B = y_n    y₀ = y_A.  Suppose that we have a collection of n line segments that connect successive pairs of points, forming a path connecting the points A and B.  Our goal is to approximate the values of the coordinates of our points by minimizing the time of travel along such a path.

Let us consider how to find the time of travel, t_i, along the i^th line segment, from the point (x_i-1, y_i-1) to the point (x_i, y_i).  If the velocity along the i^th line segment is constant, then we can find the time it takes for the bead to traverse this part of the path by means of the relation

                distance = rate * time.

However, since our velocity is not constant along this line segment, we can at best use this relationship to approximate t_i:

                d_i    v_ave,i t_i,

where d_i is the length of the i^th line segment and v_ave,i is an average velocity of the bead along this line segment.  With this relation we shall derive a rough means of producing a simple approximation of t_i from the values of the coordinates of the endpoints of the i^th line segment.

The value d_i is just the distance from (x_i-1, y_i-1) to (x_i, y_i), and is given by the quantity

                d_i = [(x_i - x_i-1)² + (y_i - y_i-1)²]^1/2.

By placing smooth joints between successive line segments, the momentum of the bead will be conserved, and so will be transferred from one line segment to the next in such a way that the tangential direction of velocity will be directed along the length of the current line segment.  Now, the magnitude of the velocity v_i of the bead at the point (x_i, y_i) can be found from the magnitude of the kinetic energy of the bead at this point:

                K = .5 m v_i²,

where m is the mass of the bead.  Since the bead is released from rest, the gain in kinetic energy of the bead must be the same as the bead's loss in potential energy, the magnitude of which is given by

                P = m g h_i,

where g = 9.81 m/sec² is the acceleration due to gravity, and  h_i is the change in height of the bead from its initial position (0, y₀) to its current position (x_i, y_i):

                 h_i = y₀ - y_i.

Thus, setting these two equations equal to each other and simplifying, we find that the magnitude of the velocity of the bead at the point (x_i, y_i) is determined by

                v_i² = 2 g (y₀ - y_i).

By taking the square root of both sides we obtain

                v_i = [2 g (y₀ - y_i)]^1/2.

The velocity of the bead along the i^th line segment, from the point (x_i-1, y_i-1) to the point (x_i, y_i), does not change in a linear manner with respect to height y, and so finding the average velocity along this line segment, v_ave,i, is not a simple task.  However, if the change in height across this wire is small, then we can approximate the average velocity of the bead along the length of the wire with the average of the initial and final velocities of the bead on the wire:

                v_ave,i    (v_i-1 + v_i) / 2.

Now, we approximate the time of travel along the i^th line segment with the relation

                t_i    d_i / v_ave,i    2 d_i / (v_i-1 + v_i),

which can be expressed entirely in terms of the known value y₀ = y_A, and the additional values x_i-1, y_i-1, x_i, and y_i.  The total amount of time required to traverse the n line segments that make up our path is then

                t₁ + t₂ + ··· + t_n,

which we approximate with the sum

2 d₁ / (v₀ + v₁) + 2 d₂ / (v₁ + v₂) + ··· + 2 d_n / (v_n-1 + v_n).
By finding the values of x₁, x₂, ..., x_n-1 and y₁, y₂, ..., y_n-1 that minimize this sum, we find a path that approximates the solution of the original problem.

Given a function f(z₁, z₂, ..., z_n) of the variables z₁, z₂, ..., z_n, the variable z₁ lying in the interval [zmin₁, zmax₁], the variable z₂ lying in the interval [zmin₂, zmax₂], ..., and the variable z_n lying in the interval [zmin_n, zmax_n], Mathematica can be utilized to find the values of each of z₁, z₂, ..., z_n, in their respective intervals, such that the function f(z₁, z₂, ..., z_n) yields a minimum value.  This is done by means of the command FindMinimum as follows:

   FindMinimum[f(z₁,z₂,...,z_n),{z₁,zsd₁,zmin₁,zmax₁},
    {z₂,zsd₂,zmin₂,zmax₂},...,{z_n,zsd_n,zmin_n,zmax_n}]

Here we use the values zsd₁, zsd₂, ..., zsd_n as (optional) seed values--initial values for each of the z₁, z₂, ..., z_n from which Mathematica begins its search.  Note that this command will yield the output of a minimum value of f(z₁, z₂, ..., z_n) on its domain, and the values of z₁, z₂, ..., z_n at which this minimum occurs, but it will not necessarily yield the overall minimum value of the function on its domain.



Example 1:    Let us approximate the shape of the brachistochrone from the point (0 meters, 2 meters) to the point (2 meters, 1 meter) by representing the path of minimum time as four line segments that span equal horizontal intervals.

For this problem, we will find four line segments of equal (horizontal) width--meaning that our line segments will span from x = 0 to x = .5, from x = .5 to x = 1, from x = 1 to x = 1.5, and from x = 1.5 to x = 2--such that the approximate time of travel from (0, 2) to (2, 1) is minimized.  We know that the first line segment begins at the point (0, 2) and the last line segment ends at (2, 1).  However we do not know the y-coordinates of the endpoints that correspond to the values x = .5, x = 1, or x = 1.5.  Let us designate these respective values as y₁, y₂, and y₃.  Thus the endpoints of our line segments are (0, 2), (.5, y₁), (1, y₂), (1.5, y₃), and (2, 1).  The goal of the problem is then to find the values of y₁, y₂, and y₃ such that our approximate time of travel across these line segments, from (0, 2) to (2, 1), is minimized.

We shall be utilizing the Mathematica command FindMinimum to generate our approximate solution.  However, before we do this, we need to define some functions.  First we shall specify the constant g to represent the acceleration due to gravity:

   In[1]:= g=9.81;

We also input the known values of the endpoints of our line segments:

   In[2]:= x0=0.0;x1=.5;x2=1.0;x3=1.5;x4=2.0;
           y0=2.0;y4=1.0;

We now define the function that represents the approximate average velocity over the interval of height from y = s₁ to y = s₂:

   In[3]:= vave[s1_,s2_]=(Sqrt[2*g*(y0-s1)]
            +Sqrt[2*g*(y0-s2)])/2;

The following function is our distance function, finding the distance between the points (r₁, s₁) and (r₂, s₂):

   In[4]:= d[r1_,s1_,r2_,s2_]=Sqrt[(r2-r1)^2+(s2-s1)^2];

Finally, we input the function that approximates the time of travel between the points (r₁, s₁) and (r₂, s₂):

   In[5]:= t[r1_,s1_,r2_,s2_]=d[r1,s1,r2,s2]/vave[s1,s2];

Now that these functions have been input, we find the values of our unknowns that minimize the approximate time of travel from (0, 2) to (2, 1) as follows:

   In[6]:= FindMinimum[t[x0,y0,x1,y1]+t[x1,y1,x2,y2]
            +t[x2,y2,x3,y3]+t[x3,y3,x4,y4],{y1,1.75,0.0,2.0},
            {y2,1.5,0.0,2.0},{y3,1.25,0.0,2.0}]

   Out[6]= {0.828647,{y1->1.20939,y2->1.00609,y3->0.942089}}

We have specified the seed values of the unknowns in such a way that the endpoints of our line segments would all lie along the line from (0, 2) to (2, 1).  Note that this solution indicates that we need y₁ = 1.20939, y₂ = 1.00609, and y₃ = 0.942089 in order to minimize the approximate time of travel from (0, 2) to (2, 1).  This minimum approximate time is given as .828647 seconds, but we cannot be certain about its accuracy since we do not know about the accuracy of our method of approximation.

We obtain a graph of our solution by first specifying our points

   In[7]:= points={{0.0,2.0},{.5,1.20939},{1.0,1.00609},
            {1.5,0.942089},{2.0,1.0}};

and then by utilizing a ListPlot command with the PlotJoined option set at True to obtain the appropriate collection of line segments.

   In[8]:= ListPlot[points,PlotJoined->True,
            PlotRange->{0,2.1}]
   Out[8]= -Graphics-

Note that the actual brachistochrone that satisfies this problem is given (numerically) by the set of parametric equations

                x(s) = 0.517195 (s - sin(s)),

                y(s) = 1.4828 + 0.517195 cos(s),

for 0    s    3.50837.  When we plot this curve along with our solution, we obtain the following:
Thus it appears that our approximation is a fairly good representation of the actual solution.





Example 2:    Let us now approximate the shape of the brachistochrone from the point (0 meters, 2 meters) to the point (2 meters, 1 meter) by representing the path of minimum time as four line segments that span equal vertical intervals.

We now seek four line segments that span from y = 2 to y = 1.75, from y = 1.75 to y = 1.5, from y = 1.5 to y = 1.25, and from y = 1.25 to y = 1 such that the time of travel from (0, 2) to (2, 1) is minimized.  The first line segment begins at (0, 2) and the last line segment ends at (2, 1).  However we do not know the x-coordinates of the endpoints that correspond to the values y = 1.75, y = 1.5, or y = 1.25.  Let us designate these respective values as x₁, x₂, and x₃.  The endpoints of our line segments are then (0, 2), (x₁, 1.75), (x₂, 1.5), (x₃, 1.25), and (2, 1).  The goal of the problem is then to find the values of x₁, x₂, and x₃ such that our approximate time of travel across these line segments, from (0, 2) to (2, 1), is minimized.

Before we begin the computation we must clear the values of the x- and y-coordinates of the endpoints that we specified in Example 1.  We can clear all variables that begin with the letters x or y with the command

   In[9]:= Clear["x*","y*"]

Now we specify our known x and y values

  In[10]:= x0=0.0;x4=2.0;
           y0=2.0;y1=1.75;y2=1.5;y3=1.25;y4=1.0;

Assuming that our functions have already been defined from the previous example, we can now obtain our solution with the command

  In[11]:= FindMinimum[t[x0,y0,x1,y1]+t[x1,y1,x2,y2]
            +t[x2,y2,x3,y3]+t[x3,y3,x4,y4],{x1,.5,0.0,2.0},
            {x2,1.0,0.0,2.0},{x3,1.5,0.0,2.0}]

  Out[11]= {0.816944,{x1->0.0683342,x2->0.27427,
            x3->0.644697}}

This yields the values x₁ = 0.0683342, x₂ = 0.27427, and x₃ = 0.644697.  Thus the endpoints of our line segments are (0, 2), (0.0683342, 1.75), (0.27427, 1.5), (0.644697, 1.25), and (2, 1).  Note that the approximate time of travel is smaller than the value that we obtained in Example 1, but since both values are approximations we do not know which path actually yields a lower value.  With an input of our points

  In[12]:= points={{0.0,2.0},{0.0683342,1.75},
            {0.27427,1.5},{0.644697,1.25},{2.0,1.0}};

we can obtain a graph of the solution

  In[13]:= ListPlot[points,PlotJoined->True,
            PlotRange->{0,2.1}]
  Out[13]= -Graphics-

Note that, by specifying the y-coordinates of the endpoints of our line segments, we gave ourselves more control over the shape of the result.  In this case we kept the solution from dropping below the height y = 1.  Because we have no information about the vertical range of the solution curve before we do the computation, specifying the y-coordinates of the endpoints of our line segments, such as we did with this problem, is not advisable.





Example 3:    Let us approximate the shape of the brachistochrone from the point (0 meters, 2 meters) to the point (2 meters, 1 meter) by representing the path of minimum time as four line segments, but with no specifications on either the x- or y-coordinates of the endpoints of the line segments.

In this case we have four line segments whose endpoints we take to be (0, 2), (x₁, y₁), (x₂, y₂), (x₃, y₃), and (2, 1), with 0 < x₁ < x₂ < x₃ < 2.  Here we specify none of the values x₁, x₂, x₃, y₁, y₂, or y₃.  The goal of the problem is then to find values for these unknowns such that our approximate time of travel across these line segments, from (0, 2) to (2, 1), is minimized.

Once again we can assume that we have the appropriate functions already defined.  After we clear the values of the endpoints specified in the previous example

  In[14]:= Clear["x*","y*"]

we specify the coordinates that we are assuming:

  In[15]:= x0=0.0;x4=2.0;
           y0=2.0;y4=1.0;

We now utilize the command FindMinimum to find values for our unknowns that will yield a minimum approximate time.  Note the choice of seed values for our unknowns--we have chosen equally spaced values for the x-coordinates and a constant value for the y-coordinates.  It may be that one set of seed values will produce a solution and another set will not, so a given example may need to be reworked with different seed values if the computation does not successfully yield a solution.

  In[16]:= FindMinimum[t[x0,y0,x1,y1]+t[x1,y1,x2,y2]
            +t[x2,y2,x3,y3]+t[x3,y3,x4,y4],{x1,0.5,0.0,2.0},
            {x2,1.0,0.0,2.0},{x3,1.5,0.0,2.0},{y1,1.0,0.0,2.0},
            {y2,1.0,0.0,2.0},{y3,1.0,0.0,2.0}]

  Out[16]= {0.812228,{x1->0.0390706,x2->0.358182,x3->1.05619,
            y1->1.82105,y2->1.39645,y3->1.02852}}

Thus we see that the values x₁ = 0.0390706, x₂ = 0.358182, x₃ = 1.05619, y₁ = 1.82105, y₂ = 1.39645, and y₃ = 1.02852 yield a minimum approximate time of travel, with a value that is lower than the values obtained in Examples 1 and 2.  To obtain a plot of the graph of this solution, we input our points

  In[17]:= points={{0.0,2.0},{0.0390706,1.82105},
            {0.358182,1.39645},{1.05619,1.02852},{2.0,1.0}};

and then utilize the ListPlot command:

  In[18]:= ListPlot[points,PlotJoined->True,
            PlotRange->{0,2.1}]
  Out[18]= -Graphics-




Exercises

---

Last modified: Thu Feb 15 2001