Mathematica(l) Roller Coasters

Selwyn Hollis
Department of Mathematics
Armstrong Atlantic State University
Savannah, GA 31419
shollis@armstrong.edu

Introduction

Mathematica animations that simulate roller coasters provide an interesting and fun way of putting basic ideas from multivariable calculus and elementary physics to work. In this article we will describe how to create “monorail” roller coaster animations in which the track is a surface in (ℝ^3) in the form of a narrow strip. The space curve comprising the center “rail” of the track may be specified either as a closed “coordinate curve” on a given parametric surface or as a closed space curve. In the first case, the orientation of the vehicle makes use of partial derivatives. In the second case, the orientation of the vehicle makes use of the curve’s classical TNB frame, which consists of the tangent, normal, and binormal vectors.

The Physics

The motion of our simulated roller coaster will be made to appear realistic by our use of the principle of conservation of energy. We will assume a frictionless track and that upward is the positive vertical direction. So at all times we have

Δ ( KE   +  PE ) = 0,

where KE = m (v^2)/2 is kinetic energy, and PE = m g z is potential energy. Here m is mass, g is the gravitational acceleration constant, z is vertical displacement from any reference point, and v is speed. Since there is no force—and consequently no change in potential energy—in any horizontal plane, we actually have

KE = m (Subscript[v, z]^2)/2+c

where (Subscript[v, z]) is the vertical component of the velocity and c is a constant. So from the conservation principle above, it follows easily that

Δ ((Subscript[v, z]^2)) ∝ Δ z.

Suppose that at the highest point on the roller coaster’s track we have (Subscript[v, z]) = (Subscript[v, top]) and z = (Subscript[z, top]). Then at any point on the track,

(Subscript[v, z]^2)-(Subscript[v, top]^2)  = k( (Subscript[z, top]) - z),

where k > 0. The result of all this is that if the position of our roller coaster as a function of time t is given by R(t) = (x(t), y(t), z(t)), then

(((x^′)(t))^2)+(((y^′)(t))^2) is constant, and  ((z'(t))^2)=(Subscript[v, z]^2) = k((Subscript[z, top]) - z(t)) + (Subscript[v, top]^2).

Determining  Position

The center rail of the track will be a closed space curve with a given parametrization

r(u)=((Subscript[r, 1])(u), (Subscript[r, 2])(u), (Subscript[r, 3])(u)).

This parameterization will only describe the path of the vehicle, not its position as a function of time. Since an animation will consist of a sequence of frames that show the roller coaster at equally-spaced times

0,Δ t,2Δ t, 3 Δ t, …,

we need to determine the position of the vehicle at each of these times. This will be done by computing a corresponding sequence of values of the parameter u, which will entail recomputing the parameter step Δ u for each successive frame. The method for computing each successive Δ u is developed as follows.

Suppose that the actual position of the vehicle at time t is given by R(t)=(x(t), y(t), z(t)). Then, thinking of u as a function of t, we have

R(t)=r(u(t)),

from which follows

(R^′)(t)=d(r/d)ud(u/d)t

by the Chain Rule. We assume that d(u/d)t>0, so it follows that

d(u/d)t=‖(R^′)(t)‖/‖d(r/d)u‖ .

Now, since (((x^′)(t))^2)+(((y^′)(t))^2) is constant and (((z^′)(t))^2)=(Subscript[v, z]^2)=k((Subscript[z, top]) - z(t))+(Subscript[v, top]^2), this becomes

d(u/d)t=(Sqrt[k((Subscript[z, top]) - z(t))+c ])/‖d(r/d)u‖ ,

where c is a constant. Thus, the variable parameter step corresponding to the fixed time step Δ t can be estimated as

Δ u≈((Sqrt[k((Subscript[z, top]) - (Subscript[r, 3])(u))  +  c ])/(‖ r'(u) ‖))Δ t ,

which can be computed without knowledge of R(t). Since Δ t is fixed, we will introduce new constants \(TraditionalForm\`k\&~\) and ε and rewrite this simply as

Δ u≈ \(k\&~\)((Sqrt[(Subscript[z, top])- (Subscript[r, 3])(
                      u) + ε ])/(‖ r'(u) ‖)) .

The apparent speed of the vehicle in our animations will be controlled by adjusting the constants \(TraditionalForm\`k\&~\) and ε.

The Vehicle and Its Orientation

The vehicle will be represented by a flexible tube connected to the center rail of the track. Its cross-sections orthogonal to the track will be circles that intersect the track’s center rail. At each point on the center rail there is a plane that is orthogonal to the path at that point and contains a circular cross-section of the tube. In that plane, imagine two coordinate axes, a “horizontal” axis that is parallel to the track surface and a “vertical” axis that is orthogonal to the track. These axes are indeed horizontal and vertical with respect to the orientation of the vehicle.

So in order to keep the vehicle oriented properly at each location along the track, we need to compute a pair of orthogonal vectors that point in the directions of the axes just described. These vectors are:

(i)  N(u), a unit vector that is orthogonal to the center rail and parallel to the track surface;

(ii)  B(u), a unit vector that is orthogonal to both the center rail and N(u).

With these vectors defined, we can now describe the surface we will use to represent the vehicle at any point r(u) on the track. Its parametrization is

Φ(s, θ) = 
        r(s) + ρ B(
              s) (1 + cos θ) + ρ sin θ N(
              s),  u - ℓ/2 ≤ s ≤ 
        u + ℓ/2,  
      0 ≤ θ ≤ 2  π .

where ρ is the radius of the vehicle’s circular cross-sections, and ℓ is the vehicle’s length.

A Track Given  by a Closed Coordinate  Curve on a Parametric Surface

A nice example is provided by a Möbius strip, an example of which is parametrized and plotted as follows.

q[u_, v_] := {(2 - v Sin[u/2]) Cos[u], 
          v Cos[u/2], (2 - v Sin[u/2]) Sin[u]};

track = 
        ParametricPlot3D[
          Evaluate[q[u, v]], {u, 0, 2  π}, {v, - .2,  .2}, 
          PlotPoints -> {60, 3}, Boxed -> False, 
          Axes -> False];

The Vector Computations

The track’s center rail is the “coordinate curve” parametrized by

r[u_] = q[u, 0]

(By “coordinate curve”  we mean a curve with parametrization given by setting one of the coordinates of the surface parametrization equal to a constant.)

The corresponding unit tangent vector  is

unitTan[
        u_] = (r^′)[u]/(Sqrt[(r^′)[u] . (r^′)[u]]) // Simplify

A vector that is parallel to the track surface and not parallel to the curve is given by ∂(q/∂)v evaluated at v = 0.

vec2[u_] = (∂_v) q[u, v] /. v -> 0 // Simplify

One step of the Gram-Schmidt process produces a vector that is parallel to the track surface and orthogonal to the curve:

nrml[u_] = vec2[u] - vec2[u] . unitTan[u] unitTan[u]

The corresponding unit vector is

unitNrml[u_] = nrml[u]/(Sqrt[nrml[u] . nrml[u]]) // TrigExpand // 
        Simplify

This is the vector denoted above by N(u). (Note that in this example, the last two steps produced no change. These steps are needed in general, however.)

Finally, the cross product T(u)×N(u) produces the vector denoted above by B(u):

biNrml[u_] = unitTan[u]\[Cross]unitNrml[u] // Simplify

The Vehicle

The following produces the vehicle at any point along the curve. Below that is a sample plot using u = 1.5.

tube[u_] := 
        ParametricPlot3D[
          Evaluate[
            r[s] +  .15 ((1 + Cos[v]) biNrml[s] + 
                    Sin[v] unitNrml[s])], {s, u -  .2, u +  .2}, {v, 0, 
            2  π}, PlotPoints -> {5, 9}, 
          DisplayFunction -> Identity];

Show[tube[1.5],DisplayFunction -> $DisplayFunction];

The Animation

Now we are ready to compute the frames for the animation. The first line below defines the function that computes the parameter step Δ u. The numbers  .02 and  .67 there were chosen experimentally. Forty-eight frames will be produced, the first of which is shown below.

ztop = 2;
uStep[
          u_] =  .67 ((Sqrt[ ztop - r[u][[3]] +  .02])/(Sqrt[(r^′)[u] . (r^′)[u]]));
u = 0;While[u ≤ 4  π,u = u + uStep[u];
Show[tube[u],track,
Axes->False,PlotRange -> {{-2.25, 2.25}, {-1.1, 
                      1.1}, {-2.25, 2.25}},
Boxed -> False,DisplayFunction -> $DisplayFunction]];Clear[u]

A Track Given  by a Closed  Space Curve

An interesting track is obtained from the curve parametrized by

r[u_] := {-(2 + Cos[u]) Sin[2  u], (2 + Sin[u]) Cos[2  u], 
        2  Cos[u]}

ParametricPlot3D[Evaluate[r[u]], {u, 0, 2  π}, 
        ViewPoint -> {3, -2, 1}];

The Vector Computations

In this situation the vectors of interest comprise the classical TNB frame at each point along the curve. Because even fairly simple curves give rise to very large and complicated symbolic expressions for the necessary vectors, we will not calculate them explicitly.

The tangent vector of interest here is

tanV[u_] = (r^′)[u] // TrigExpand // Simplify

The corresponding unit tangent vector T(u) is

unitTan[u_] := Module[{v = tanV[u]}, (v/(Sqrt[v . v]))]

We will approximate the derivative of this unit  tangent vector by a central difference approximation:

dT[u_] := (unitTan[u +  .0005] - 
            unitTan[u -  .0005])/ .001

Now, the principal unit normal vector N(u) is

unitNrml[u_] := Module[{n = dT[u]}, (n/(Sqrt[n . n]))]

Finally the cross product T(u)×N(u) produces the binormal vector B(u):

biNrml[u_] := unitTan[u]\[Cross]unitNrml[u]

The Track and the Vehicle

To produce the track from the curve, we use the vector N(u) to create a surface parametrized by

Γ(u, v) = 
        r(u) + v N(u),   0 ≤ u ≤ 
        2  π,  -w/2 ≤ 
        v ≤ w/2 ,

where w is the width of the track.

track = 
        ParametricPlot3D[
          r[u] + v unitNrml[u], {u, 0, 2  π}, {v, - .2,  .2}, 
          ViewPoint -> {3, -2, 1}, 
          PlotPoints -> {100, 3}];

The following produces the vehicle at any point along the curve. Below that is a sample plot using u = 0. (Note the absence of the Evaluate function inside of ParametricPlot3D here. This is because of the way the unit normal and binormal vectors were defined above with delayed evaluation.)

tube[u_] := 
        ParametricPlot3D[
          r[s] +  .2  biNrml[s] (1 + Cos[θ]) +  .2  Sin[θ] 
              unitNrml[s], {s, u -  .15, u +  .15}, {θ, 0, 2  π}, 
          PlotPoints -> {5, 9}, 
          DisplayFunction -> Identity];

Show[tube[0], DisplayFunction -> $DisplayFunction];

The Animation

This animation is created in essentially the same way as the previous one. Forty frames will be produced with the particular parameter values used. The GraphicsArray that follows contains frames 1, 14, 20, and 22.

ztop = 2;
uStep[
          u_] = ((Sqrt[ztop - r[u][[3]] +  .02])/(Sqrt[(r^′)[u] . (r^′)[u]]));
u = 0; 
    While[u ≤ 2  π, 
      u = u + uStep[u]; 
Show[tube[u], track, 
        ViewPoint -> {3, -2, 1}, Axes -> False, 
        PlotRange -> {{-3.6, 3}, {-3.3, 3.3}, {-2.1, 
              2.4}}, 
DisplayFunction -> $DisplayFunction]];  Clear[u]