Text provided under a Creative Commons Attribution license, CC-BY, code provided under an MIT license. (c) Aron J. Ahmadia, Lorena A. Barba, 2014. Thanks: Gilbert Forsyth and Olivier Mesnard, and NSF for support via CAREER award #1149784.

The Parabolic Particle Path

Adapted from Tutorial for HW #4 of GWU MAE 2117 by Chunlei Liang

Let us consider a particle sliding on a parabola which has a profile of:

\begin{align} y = x^2 \end{align}

The coordinates of the particle are shown in Figure 1.

Figure 1: Coordinates



In [3]:

    
SVG("./Figure_1.svg")









    Out[3]:

The reference frames of the particle can be defined by an inertial reference frame $I = (O, \vec{e}_x, \vec{e}_y)$ and a local reference frame $B = (P, \vec{e}_t, \vec{e}_n)$ as shown in Figure 2.

Figure 2: Reference Frame



In [4]:

    
SVG("./Figure_2.svg")









    Out[4]:

We can then draw a free-body diagram for this particle as depicted in Figure 3.

Figure 3: Free-body diagram



In [5]:

    
SVG("./Figure_3.svg")









    Out[5]:

The total force includes the gravitational force $-m_pg\vec{e}_y$ and the resultant force $N\vec{e}_n$.

The kinematics of the particle can be written as the following.

The position vector of the particle in the inertial reference frame is:

\begin{align} \vec{r} = x\vec{e}_x + y\vec{e}_y = x\vec{e}_x + x^2\vec{e}_y \end{align}

The velocity vector of the particle is therefore:

\begin{align} \vec{v} = \dot{x}\vec{e}_x + \dot{y}\vec{e}_y = \dot{x}\vec{e}_x + 2x\dot{x}\vec{e}_y \end{align}

The acceleration vector of the particle can be derived as:

\begin{align} \vec{a} = \ddot{x}\vec{e}_x + 2\left(\dot{x}^2 + x\ddot{x}\right)\vec{e}_y \end{align}

Thus, according to Newton's Second Law, we have:

Equation 1

\begin{align} -m_pg\vec{e}_y + N\vec{e}_n = m_p\left\{{\ddot{x}\vec{e}_x + \left(\dot{x}^2 + x\ddot{x}\right)\vec{e}_y}\right\} \end{align}

The tangent direction unit vector shown in Figure 2 can be related to the velocity vector through $\vec{e}_t = \frac{\vec{v}}{||\vec{v}||}$

One can compute the unit normal vector:

\begin{align} \vec{e_n} = \vec{e_z} \times \vec{e}_t = \frac{-2x\vec{e}_x + \vec{e}_y}{\sqrt{1+4x^2}} \end{align}

Equation (1) can now be split into two scalar equations:

Equation 2

\begin{align} \frac{-2xN}{\sqrt{1+4x^2}} = m_p\ddot{x} \end{align}

Equation 3

\begin{align} \frac{N}{\sqrt{1+4x^2}} -m_pg = m_p2(\dot{x}^2 + x\ddot{x}) \end{align}

Combine Equations (2) and (3) and eliminate N to obtain the equation of motion for the particle:

Equation 4

\begin{align} \ddot{x} = \frac{-2x(g+2\dot{x}^2)}{1+4x^2} \end{align}

Equation (4) can be split into two first-order ODEs

Equation 5

\begin{align} \dot{x} &= u \end{align}

Equation 6

\begin{align} \dot{u} &= \frac{-2x(g+2u^2)}{1+4x^2} \end{align}

If we introduce the vector $\vec{q} = \left[\begin{array}{c}x\\u\end{array}\right] = \left[\begin{array}{c}q_0\\q_1\end{array}\right]$ we can write equations (5) and (6) in terms of $\vec{q}$.

Code Fragment 1



In [1]:

    
# As usual, we begin by importing the NumPy package
import numpy as np

# Then, we define the gravitational constant
g = 9.8

# Now, we define the derivative of our vector Q
def Q_Dot(Q):
    """Calculate the derivative of Q"""
    x, u = Q[0], Q[1]
    x_dot = u
    u_dot = -2*x*(g+2*u**2)/(1+4*x**2)
    return np.asarray([x_dot, u_dot])

Go ahead and execute the above cell fragment. None of the three actions (importing a package, defining a variable, defining a function) produce any output, so IPython lets us know that the cell has executed by adding a number between the In [X] marker to the left of the cell.

Recall that to execute a code fragment, we need to activate the cell by clicking on it, then either pressing the "Play" button, or holding down the Shift key and pressing Enter.

Let's check our derivative function at a few initial conditions. We'll do that with a simple test function that simply prints input and output values. Once we've defined the test function, we can easily use it for a range of inputs without having to type the same print statements again and again.

Code Fragment 2



In [2]:

    
def check_Q_Dot(x_0, u_0):
    """Test function for Q_Dot"""
    print "initial velocity                 : %g" % x_0
    print "initial acceleration             : %g" % u_0
    x_dot, u_dot = Q_Dot([x_0, u_0])
    print "computed velocity derivative     : %g" % x_dot
    print "computed acceleration derivative : %g" % u_dot

First, if the parabola is at rest ($x=0, u=0$), we don't expect there to be any changes in velocity or acceleration.



In [3]:

    
check_Q_Dot(0.0, 0.0)









    



initial velocity                 : 0
initial acceleration             : 0
computed velocity derivative     : 0
computed acceleration derivative : -0

Good! Note that the acceleration derivative is zero, but has a negative sign caused by the slight rounding in our computations. Let's slightly improve our test. The computed velocity derivative should always be equal to the current acceleration, so let's make it an error when it doesn't!

Code Fragment 3



In [4]:

    
def check_Q_Dot(x_0, u_0):
    """Test function for Q_Dot"""
    print "initial velocity                 : %g" % x_0
    print "initial acceleration             : %g" % u_0
    x_dot, u_dot = Q_Dot([x_0, u_0])
    print "computed velocity derivative     : %g" % x_dot
    print "computed acceleration derivative : %g" % u_dot
    assert(x_dot == u_0)
    
check_Q_Dot(0.0, 2.0)









    



initial velocity                 : 0
initial acceleration             : 2
computed velocity derivative     : 2
computed acceleration derivative : -0

Numerical Checkpoint 1

Use equations (5) and (6) to calculate $\dot{x}$ and $\dot{u}$, given $x_0=1$ and $u_0=1$. Check your work using the code cell below.



In [5]:

    
#check_Q_Dot( , )

Euler's Method

Now that we have our transformed our second-order scalar ODE into a system of first-order ODEs, we can solve this system using a numerical integrator.

The simplest integrator is Euler's method. Recall that if we know that:

\begin{align} \dot{y}(t) &= f(t,y(t)) \\ y(t_0) &= y_0 \end{align}

then we can solve for $y$ at a set of discrete time steps of size $h$, such that $t_n = t_0 + n*h$. Euler's method relies on the recurrence relation:

\begin{align} y_{n+1} = y_n + hf(t_n,y_n) \end{align}

Code Fragment 4



In [6]:

    
def euler_step(f, y_n, n, h):
    """Take a single step using Euler's method.

    y_n - The current value of the vector y at t=n
    f   - function with with calling form f(y,t), computes the derivative of y at t 
    """
    
    y_dot = f(y_n, n)
    y_n_plus_1 = y_n + h*y_dot
    return y_n_plus_1

As always, we verify with a test example. Let's assume we are observing simple particle motion, with initial position $y(0)=0$ and constant velocity $\dot{y}(t) = 10$. If we take a time step $h=0.5$, we expect $y(0.5) = 5$.

Code Fragment 5



In [7]:

    
def constant_velocity_10(y_n, n): 
    """Derivative of particle's position that is at a constant velocity of 10"""
    return 10     
euler_step(constant_velocity_10, 0, 0, 0.5)









    Out[7]:





5.0

We could have also used a lambda statement here to define the constant_velocity_10 function, but function definitions in Python are so easy that we just created another one!

Numerically Integrating

Now that we have a way to take a single time step, we'd like to be able to take several of them in order to calculate $y$ at some final time, $t_{final}$. To do that, we'll need to use a for loop, once for each time step we'd like to take. More sophisticated numerical integration techniques take time steps of varying length, but for this exercise we'll take $n=\frac{t_{final}-t_0}{h}$ Euler steps.

Code Fragment 6



In [8]:

    
def simple_integrator(f, y_0, t, h):
    """Compute trajectory of values of y using Euler's method
    
    f       - Computes the derivative of y at t, with calling form f(y,t)
    y_0     - Initial value of vector y
    t       - t_0, t_final (distance should be a multiple of h)
    h       - Requested distance between time steps
    """
    
    t_0, t_final = t
    # n - number of steps to take
    n = int((t_final-t_0)/h)

    # d - dimension of the problem, equal to number of values in vector y
    d = len(y_0)
    
    
    # y_traj - trajectory of y at each step i, represented as [i,:]
    # The trajectory stores y at every time, including the initial time t_0.
    # y_traj will store the output of $n$ steps + 1 for the initial value.

    y_traj = np.zeros((n+1,d))

    # t  - complete set of values of t at each step i
    t = np.linspace(t_0, t_final, n+1)
    
    y_traj[0,:] = y_0
    
    for i in range(n):
        y_traj[i+1,:] = euler_step(f, y_traj[i,:], t[i], h)
    return t, y_traj

Whew, that might seem like a lot of code, but the majority of the code fragment is comments, and as long as you're comfortable with for loops in Python and using np.zeros to create empty arrays and np.linspace to create sets of equally spaced points, this function should make sense.

Let's test it out, using a visual inspection this time. Again, the following two lines of code tell IPython that we want to view our output figures in the Notebook and that we'll be using plt as shorthand for matplotlib.pyplot.

Code Fragment 7



In [9]:

    
%matplotlib inline
import matplotlib.pyplot as plt

One final trick. Our general integrator calls functions that might do something with the current time (that is, they are functions of $y(t)$ and $t$). Since our Q_Dot only relies on $y$, we need to create an adapter function around Q_Dot that acts like it will use $t$, but really doesn't. Here's how to do that:

Code Fragment 8



In [10]:

    
def Q_Dot_t(y, t):
    """Adapter function to allow us to use general integrator with Q_Dot"""
    return Q_Dot(y)

The sign of well-designed code is that all of the complexity is hidden from the rest of the program. Here's our call to simple_integrator for a test problem with initial velocity of 0.5 and initial acceleration of 0.5 from time 0 to 3, with time steps of size 0.001.

Code Fragment 9



In [11]:

    
t, Q_Traj = simple_integrator(Q_Dot_t, [0.5, 0.5], [0, 3], 0.001)

plt.plot(t, Q_Traj[:,0]);

Instructor Notes / Code



In [12]:

    
from IPython.core.display import HTML
def css_styling():
    styles = open("../styles/custom.css", "r").read()
    return HTML(styles)
css_styling()









    Out[12]:



In [1]:

    
from IPython.core.display import SVG

If needed, you can modify any of the Figures used in this IPython Notebook by modifying the Asymptote code below. This code is not distributed to students.



In [13]:

    
%install_ext http://raw.github.com/jrjohansson/ipython-asymptote/master/asymptote.py









    



Installed asymptote.py. To use it, type:
  %load_ext asymptote



In [14]:

    
%load_ext asymptote



In [15]:

    
%%asymptote Figure_1

defaultpen(1.5);
import graph;
size(200,200);

real f(real x) {return x^2;}
pair F(real x) {return (x,f(x));}

xaxis();
yaxis();

draw(graph(f,-2.5,2.5,operator ..));

pair zx=(2,0);
pair zy=(0,4);
pair P=(2,4);

dot("$P$",P);
draw(zx--P,red+dashed);
draw(zy--P,red+dashed);
label("$x$",(0.5*zx.x,4),N);
label("$y$",(2,0.5*zy.y),E);









    Out[15]:



In [16]:

    
%%asymptote Figure_2

defaultpen(1.5);
import graph;
size(200,200);

real f(real x) {return x^2;}
pair F(real x) {return (x,f(x));}

xaxis("$\vec{e}_x$");
yaxis("$\vec{e}_y$");

draw(graph(f,-2.5,2.5,operator ..));

pair P=(2,4);
pair et=P+0.5*(1,4);
pair en=P+0.5*(-4,1);

dot("$P$",P);
draw(P--et,blue,EndArrow(size=10));
label("$\vec{e}_t$", et, SE, blue);
draw(P--en,blue,EndArrow(size=10));
label("$\vec{e}_n$", en, 3*E, blue);









    Out[16]:



In [17]:

    
%%asymptote Figure_3

defaultpen(1.5);
import graph;
size(200,200);

real f(real x) {return x^2;}
pair F(real x) {return (x,f(x));}

draw(graph(f,-2.5,2.5,operator ..));

pair P=(2,4);
pair eg=P+0.5*sqrt(17)*(0,-1);
pair en=P+0.5*(-4,1);

dot("$P$",P);
draw(P--eg,blue,EndArrow(size=10));
label("$-m_pg\vec{e}_y$", eg, SE, blue);
draw(P--en,blue,EndArrow(size=10));
label("$N\vec{e}_n$", en, 3*E, blue);









    Out[17]: