Dynamics

• Author:

• Date:

• Time spent on this assignment:

Overview

In this assignment you are going to first

• write a first-order integrator to solve differential equations (like dynamics)

• improve this to be a second-order integrator

Then you will

• compare the path of balls thrown with and without air resistance

• measure the terminal velocity of a falling object.

Show code cell content

Hide code cell content

import numpy as np
import matplotlib.pyplot as plt
import scipy
def resetMe(keepList=[]):
    ll=%who_ls
    keepList=keepList+['resetMe','np','plt','scipy']
    for iiii in keepList:
        if iiii in ll:
            ll.remove(iiii)
    for iiii in ll:
        jjjj="^"+iiii+"$"
        %reset_selective -f {jjjj}
    ll=%who_ls
    plt.rcParams.update({"font.size": 14})
    return
import datetime;datetime.datetime.now()

---------------------------------------------------------------------------
ModuleNotFoundError                       Traceback (most recent call last)
Cell In[1], line 3
      1 import numpy as np
      2 import matplotlib.pyplot as plt
----> 3 import scipy
      4 def resetMe(keepList=[]):
      5     ll=get_ipython().run_line_magic('who_ls', '')

ModuleNotFoundError: No module named 'scipy'

Warmups

Lists to numpy arrays

In a number of cases you will hav a list like a=[1.2,3.2,5.4] and want to convert it to a numpy array. To do this, you can do a=np.array(a)

resetMe()
a=[1.2,3.2,5.4]
print(type(a))
a=np.array(a)
print(type(a))

<class 'list'>
<class 'numpy.ndarray'>

numpy arrays are useful because we do basic math operations on them (+,-,*,/) among other things we’ll see.

---

In this exercise, we are going to have a list of numpy arrays.

positions=[np.array([0.0,1.0]),
           np.array([0.1,2.0]),
           np.array([0.2,3.0])]

We are going to convert that into a two dimensional array using np.array(positions).

🦉Go ahead and make this conversion and then figure out how to seperately access the x positions and y positions (i.e. positions[:,1])

positions=[[0.0,1.0],[0.1,2.0],[0.2,3.0]] # this is an array of arrays
positions=np.array(positions) # convert to numpy array

print("shape of positions:", positions.shape, " first dimension is number of points, second dimension is x/y")


print(positions[:,1]) # y positions for all the points
print("The x positions should be [0.,0.1,0.2]")
print("The y positions should be [1.,2.,3.]")

shape of positions: (3, 2)  first dimension is number of points, second dimension is x/y
[1. 2. 3.]
Thex x positions should be [0.,0.1,0.2]
The y positions should be [1.,2.,3.]

Math with numpy arrays

Imagine that I have a numpy array representing a velocity: velocity=np.array([0.2,2.5]).

I can now do math with this array in an intuitive way. I could for example take a time dt=0.1 and multiply it by dt. This will multiply both the “0.2” and “2.5”.

• Multiply by a scalar value dt: velocity*dt

• Square: velocity**2

• logarithm: np.log(velocity)

• sin: np.sin(velocity)

• Add two vectors: x+velocity*dt etc:

All these operations are done element-wise; meaning that they are applied to each element in turn.

resetMe()
dt=0.1
velocity=np.array([0.2,2.5])
print("dt =",dt)
print("velocity is\t",velocity)
print("dt*velocity is\t" ,dt*velocity)
print("velocity**2 is \t",velocity**2)

dt = 0.1
myVel is	 [0.2 2.5]
dt*myVel is	 [0.02 0.25]
myVel**2 is 	 [0.04 6.25]

Fitting Lines

Suppose you have two arrays, x and y and you expect that there should be some linear relationship between \(x\) and \(y\).

slope=1.35
x = np.linspace(start=0,stop=2,num=40) # This will create a 1D array of shape (40,)
print("size of x:",x.shape)
# you can also call this as just np.linspace(0,2,40)
# slope*x is the underlying data, and then we add uniform noise from -0.2 to +0.2
y = slope*x + 0.4*(np.random.random(x.shape[0])-0.5) 

plt.plot(x,y,'.')
plt.xlabel("x")
plt.ylabel("y")
plt.title("Some Data")

# an alternative way to plot using objects
fig,ax=plt.subplots()
ax.plot(x,y,'.')
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_title("Data plotted using objects")

size of x: (40,)

Text(0.5, 1.0, 'Data plotted using objects')

Now let’s pretend that we didn’t know the original slope and we’d like to extract it.

myLine=np.polyfit(x,y,1). If you then print out myLine you should see that it gives you a list where the first value is the slope and the second value the y-intercept.

By the way, if you want to understand the arguments of the

# try it here
fit = np.polyfit(
    x,
    y,
    1 # degree of the polynomial -- this is a linear fit
) # returns the polynomial coefficients in descending order
print(f"The slope is {fit[0]} and the intercept is {fit[1]}")

help(np.polyfit) # help lets you understand how to use a function

plt.plot(x,y,'.')
plt.plot(x,fit[0]*x+fit[1])
# you can also do np.polyval(fit,x) # evaluates the polynomial at the points x

plt.xlabel("x")
plt.ylabel("y")
plt.title("Some Data")
plt.show()

The slope is 1.402181819755196 and the intercept is -0.06059016172536342
Help on _ArrayFunctionDispatcher in module numpy:

polyfit(x, y, deg, rcond=None, full=False, w=None, cov=False)
    Least squares polynomial fit.

    .. note::
       This forms part of the old polynomial API. Since version 1.4, the
       new polynomial API defined in `numpy.polynomial` is preferred.
       A summary of the differences can be found in the
       :doc:`transition guide </reference/routines.polynomials>`.

    Fit a polynomial ``p(x) = p[0] * x**deg + ... + p[deg]`` of degree `deg`
    to points `(x, y)`. Returns a vector of coefficients `p` that minimises
    the squared error in the order `deg`, `deg-1`, ... `0`.

    The `Polynomial.fit <numpy.polynomial.polynomial.Polynomial.fit>` class
    method is recommended for new code as it is more stable numerically. See
    the documentation of the method for more information.

    Parameters
    ----------
    x : array_like, shape (M,)
        x-coordinates of the M sample points ``(x[i], y[i])``.
    y : array_like, shape (M,) or (M, K)
        y-coordinates of the sample points. Several data sets of sample
        points sharing the same x-coordinates can be fitted at once by
        passing in a 2D-array that contains one dataset per column.
    deg : int
        Degree of the fitting polynomial
    rcond : float, optional
        Relative condition number of the fit. Singular values smaller than
        this relative to the largest singular value will be ignored. The
        default value is len(x)*eps, where eps is the relative precision of
        the float type, about 2e-16 in most cases.
    full : bool, optional
        Switch determining nature of return value. When it is False (the
        default) just the coefficients are returned, when True diagnostic
        information from the singular value decomposition is also returned.
    w : array_like, shape (M,), optional
        Weights. If not None, the weight ``w[i]`` applies to the unsquared
        residual ``y[i] - y_hat[i]`` at ``x[i]``. Ideally the weights are
        chosen so that the errors of the products ``w[i]*y[i]`` all have the
        same variance.  When using inverse-variance weighting, use
        ``w[i] = 1/sigma(y[i])``.  The default value is None.
    cov : bool or str, optional
        If given and not `False`, return not just the estimate but also its
        covariance matrix. By default, the covariance are scaled by
        chi2/dof, where dof = M - (deg + 1), i.e., the weights are presumed
        to be unreliable except in a relative sense and everything is scaled
        such that the reduced chi2 is unity. This scaling is omitted if
        ``cov='unscaled'``, as is relevant for the case that the weights are
        w = 1/sigma, with sigma known to be a reliable estimate of the
        uncertainty.

    Returns
    -------
    p : ndarray, shape (deg + 1,) or (deg + 1, K)
        Polynomial coefficients, highest power first.  If `y` was 2-D, the
        coefficients for `k`-th data set are in ``p[:,k]``.

    residuals, rank, singular_values, rcond
        These values are only returned if ``full == True``

        - residuals -- sum of squared residuals of the least squares fit
        - rank -- the effective rank of the scaled Vandermonde
           coefficient matrix
        - singular_values -- singular values of the scaled Vandermonde
           coefficient matrix
        - rcond -- value of `rcond`.

        For more details, see `numpy.linalg.lstsq`.

    V : ndarray, shape (deg + 1, deg + 1) or (deg + 1, deg + 1, K)
        Present only if ``full == False`` and ``cov == True``.  The covariance
        matrix of the polynomial coefficient estimates.  The diagonal of
        this matrix are the variance estimates for each coefficient.  If y
        is a 2-D array, then the covariance matrix for the `k`-th data set
        are in ``V[:,:,k]``


    Warns
    -----
    RankWarning
        The rank of the coefficient matrix in the least-squares fit is
        deficient. The warning is only raised if ``full == False``.

        The warnings can be turned off by

        >>> import warnings
        >>> warnings.simplefilter('ignore', np.exceptions.RankWarning)

    See Also
    --------
    polyval : Compute polynomial values.
    linalg.lstsq : Computes a least-squares fit.
    scipy.interpolate.UnivariateSpline : Computes spline fits.

    Notes
    -----
    The solution minimizes the squared error

    .. math::
        E = \sum_{j=0}^k |p(x_j) - y_j|^2

    in the equations::

        x[0]**n * p[0] + ... + x[0] * p[n-1] + p[n] = y[0]
        x[1]**n * p[0] + ... + x[1] * p[n-1] + p[n] = y[1]
        ...
        x[k]**n * p[0] + ... + x[k] * p[n-1] + p[n] = y[k]

    The coefficient matrix of the coefficients `p` is a Vandermonde matrix.

    `polyfit` issues a `~exceptions.RankWarning` when the least-squares fit is
    badly conditioned. This implies that the best fit is not well-defined due
    to numerical error. The results may be improved by lowering the polynomial
    degree or by replacing `x` by `x` - `x`.mean(). The `rcond` parameter
    can also be set to a value smaller than its default, but the resulting
    fit may be spurious: including contributions from the small singular
    values can add numerical noise to the result.

    Note that fitting polynomial coefficients is inherently badly conditioned
    when the degree of the polynomial is large or the interval of sample points
    is badly centered. The quality of the fit should always be checked in these
    cases. When polynomial fits are not satisfactory, splines may be a good
    alternative.

    References
    ----------
    .. [1] Wikipedia, "Curve fitting",
           https://en.wikipedia.org/wiki/Curve_fitting
    .. [2] Wikipedia, "Polynomial interpolation",
           https://en.wikipedia.org/wiki/Polynomial_interpolation

    Examples
    --------
    >>> import numpy as np
    >>> import warnings
    >>> x = np.array([0.0, 1.0, 2.0, 3.0,  4.0,  5.0])
    >>> y = np.array([0.0, 0.8, 0.9, 0.1, -0.8, -1.0])
    >>> z = np.polyfit(x, y, 3)
    >>> z
    array([ 0.08703704, -0.81349206,  1.69312169, -0.03968254]) # may vary

    It is convenient to use `poly1d` objects for dealing with polynomials:

    >>> p = np.poly1d(z)
    >>> p(0.5)
    0.6143849206349179 # may vary
    >>> p(3.5)
    -0.34732142857143039 # may vary
    >>> p(10)
    22.579365079365115 # may vary

    High-order polynomials may oscillate wildly:

    >>> with warnings.catch_warnings():
    ...     warnings.simplefilter('ignore', np.exceptions.RankWarning)
    ...     p30 = np.poly1d(np.polyfit(x, y, 30))
    ...
    >>> p30(4)
    -0.80000000000000204 # may vary
    >>> p30(5)
    -0.99999999999999445 # may vary
    >>> p30(4.5)
    -0.10547061179440398 # may vary

    Illustration:

    >>> import matplotlib.pyplot as plt
    >>> xp = np.linspace(-2, 6, 100)
    >>> _ = plt.plot(x, y, '.', xp, p(xp), '-', xp, p30(xp), '--')
    >>> plt.ylim(-2,2)
    (-2, 2)
    >>> plt.show()

Fitting Polynomials

Suppose that you think you have some data which goes as \(y = Ax^b\), but you don’t know \(b\). This is typically found using a log-log plot.

Note that \(\log(y) = \log(A) + b\log(x)\).
Plotting \(\log(x)\) vs \(\log(y)\) is called plotting the the data on a log-log scale. Let’s look at how that might work.

Below is some data which has a polynomial relationship. Notice how it looks quadratic.

slope=1.35
x = np.linspace(2,30,40)
y = x**2+2.0*(np.random.random(len(x))-0.5)
plt.plot(x,y,'.')

plt.xlabel("x")
plt.ylabel("y")
plt.title("Quadratic data on a linear plot")
plt.show()

Let’s go ahead and plot x and y on a log-log scale. To do this redo the plot above but add

plt.xscale('log')
plt.yscale('log')

before the plt.show().

###  plot things on a log-scale here.
plt.plot(x,y,'.')
plt.xscale('log')
plt.yscale('log')
plt.xlabel("x")
plt.ylabel("y")
plt.title("Quadratic data on a log-log plot")
plt.show()

You should notice your plot looks roughly linear. We can go get the slope by doing a linear fit to the log-log plot. This can be done by doing

myLine=np.polyfit(np.log(x),np.log(y),1)
print(myLine)

Go ahead and try this out.

myLine=np.polyfit(np.log(x[1:]),np.log(y[1:]),1) #!#
print(myLine) #!#

[1.98207642 0.05154248]

You should see the slope recovers the \(b\) of \(y=Ax^b\).

Let’s go ahead and plot to see how well this works

plt.plot(x,y,'.')
plt.plot(x,np.exp(myLine[1])* x**myLine[0])
plt.xscale('log')
plt.yscale('log')
plt.xlabel("x")
plt.ylabel("y")
plt.show()

---

---

---

Exercise 1: Euler Integration

• List of collaborators:

• References you used in developing your code:

In this exercise we are going to learn how to simulate a ball thrown into the air. To do this, we will need to learn how to simulate dynamics.

The secret to much of physics is differential equations. Differential equations answer the following question: given the current state of your physical system, what is its state in a moment of time (\(dt\) seconds) later?

We recall that in Newtonian dynamics the following things are true:

\[F_y(t) = ma_y(t) \rightarrow a_y(t) = \frac{F_y(t)}{m} \]

\[ \frac{dv_y(t)}{dt} = a_y(t) \rightarrow dv_y(t) = a_y(t) dt \]

\[ \frac{dy(t)}{dt} = v_y(t) \rightarrow dy(t) = v_y(t) dt \]

Taking the first order in a Taylor series, we get for finite change in time \(\Delta t\):

\[v_y(t+\Delta t) \simeq v_y(t) + a_y(t) \Delta t\]

\[y(t + \Delta t) \simeq y(t) + v_y(t) \Delta t\]

a. Euler Integration in one-dimension (first order integrator)

We can use this to write a function that takes the current time, position, velocity and dt and then gives back out the new time, position, and velocity - i.e

def step_gravity_1d(t,pos,vel,dt):
  """
  Starting with the position, velocity, and current time, return the next time, position, and velocity at time t+dt.
  The force applied is -9.8 kg m/s^2 (assuming a mass of 1 and gravity). 
  Args:
     t: scalar, current time (s)
     pos: scalar, position at time t (m)
     vel: scalar, velocity at time t (m/s)
     dt: scalar, timestep (s)
  Returns
    tuple containing new time, new position, new velocity
  """
  return (new_t,new_pos,new_vel)

As our force, we will use gravity - i.e. \(F_y=-9.8m\) - and choose a mass of \(m = 1\). 🦉Write your step function. We are then going to call it five times after throwing a ball into the air and make sure you get the correct result.

The velocity you should use to update the position is not the new velocity. It is the current velocity.

🦉 Put your code here including your step function.

To test your code, use these initial parameters:

t   = 0    # set the initial time. 
vel = 2.0  # throwing it into the air at 2 m/s
pos = 1.5  #  We threw it while standing from 1.5 meters tall.
dt  = 0.01 # We will take time steps of 0.01 seconds.

and then step it five times in a row to get the answer after \(T=0.05\) seconds

t,pos,vel=Step(t,pos,vel,dt) # After 0.01 seconds
t,pos,vel=Step(t,pos,vel,dt) # After 0.02 seconds
t,pos,vel=Step(t,pos,vel,dt) # After 0.03 seconds
t,pos,vel=Step(t,pos,vel,dt) # After 0.04 seconds
t,pos,vel=Step(t,pos,vel,dt) # After 0.05 seconds

At the end you should get matching results

print("t should be 0.05,\t t="      ,round(t,4))
print("pos should be 1.5902;\t pos=",round(pos,4))
print("vel should be 1.51,\t vel="  ,round(vel,4))

(You will also include your answer here for switching the 5 steps into a loop)

#### ANSWER HERE

🦉 Modify your code so that it runs in a loop, and verify that you still get the right answers. Your loop should look something like

T=0.05
nsteps = int(np.ceil(T/dt))
for step in range(0,nsteps):
  # do stuff

b. Euler Integration in two dimensions (first order integrator)

We now have a step function that allows us to work in one-dimension. Balls are sometimes thrown in directions other then straight up, though. Let’s modify our code to also work in two dimensions. In addition to the equation for \(a_y,v_y,y\) we can also have a similar set of equations for \(a_x,v_x,x\):

• \(v_x(t+\Delta t) = v_x(t) + a_x(t) \Delta t\)

• \(x(t+\Delta t) = x(t) + v_x(t) \Delta t\)

It makes sense to write the equations for x and y in a uniform “vector” or “array” format -i.e.

• [\(v_x(t+\Delta t),v_y(t+\Delta t)\)] = [\(v_x(t),v_y(t)\)] + [\(a_x(t),a_y(t)\)] \(\Delta t\)

• [\(x(t+\Delta t),y(t+\Delta t)\)] = [\(x(t),y(t)\)] + [\(v_x(t),v_y(t)\)] \(\Delta t\)

We are now going to modify our step function to work in in two dimensions, so pos and vel will be changed to be numpy arrays of, for example, [xPosition,yPosition].

You could implement this by doing something like this:

newpos = pos.copy()
newpos[0]= pos[0]+ vel[0]*dt
newpos[1]= pos[1]+ vel[1]*dt

but equivalently you could also do

newpos = pos + vel*dt

which would work in both 2d and 3d.

🦉Rewrite your step function to work in two dimensions and try to avoid ever doing vel[0] or pos[0] to explicitly take a velocity or position in the x-direction. It should look like

def step_gravity_2d(t, pos, vel, dt):
  """
Starting with the position, velocity, and current time, return the next time, position, and velocity at time t+dt.
  The force applied is -9.8 kg m/s^2 (assuming a mass of 1 and gravity). 
  Args:
     t: scalar, current time (s)
     pos: np.array of shape (2,), position at time t (m)
     vel: np.array of shape (2,), velocity at time t (m/s)
     dt: scalar, timestep (s)
  Returns
    tuple containing new time, new position, new velocity

  """
  assert isinstance(pos, np.ndarray), "pos must be a np.ndarray"

🦉 Test use the following initial conditions

t   = 0 # set t0
dt  = 0.01 # time step 
T = 0.05 # total time to run 
vel = np.array([0.1,10.])  # initial velocity in m/s
pos = np.array([0.0,0.0])  # m

You should get the following:

print("t should be 0.05,\t t=" ,round(t,4))
print("[x,y] should be [0.005,0.4902];\t [x,y]=" ,pos)
print("[vx,vy] should be [0.1,9.51],\t [vx,vy]=" ,vel)

###ANSWER HERE

c. A full simulation

In this part we are going to write a full simulation and compute the path of a ball until it hits the ground.

🦉Write a function as follows:

def throw_ball(init_pos, init_vel, dt):
    """
    Run a simulation of a thrown ball until the ball hits the ground (at y=0).
    Gravity of -9.8 m/s^2 is applied
    Arguments:
      init_pos: (2,) array given the initial positions.
      init_velocity: (2,) array giving the initial velocity.
      dt: scalar timestep
    Returns:
     (ts, vs, pos): (nsteps,), (nsteps, 2), (nsteps,2)
    """

• Initialize your time, height and velocity

• Intialize some lists to store your time, height and velocity at each step (i.e. velocities=[])

• Loop over time steps

  • Take a step

  • Put the velocity and position into lists (i.e. velocities.append(v))

  • If the height is lower then zero, then break out of the loop (if y<0: break)

• Return np.array of time, position, and velocity

🦉 Write a function def throw_ball_exact(initPos,initVel,ts) that returns the exact positions and velocities for all times in the array ts. Use the analytic solution from introductory mechanics.

🦉 Test your functions as follows:

• Run the same five steps and initial conditions from part 1b, and confirm that throw_ball gets the same numbers.

• Plot the output of throw_ball and throw_ball_exact, and confirm that they are very close. As dt becomes very small, they should agree.

Put your code here. Using the same initial conditions as above, print the position and velocities of five steps so we can check that you have the correct answer.

###ANSWER HERE

d. Making Plots

We want to plot

• \(y\) vs. \(t\) (i.e. plt.plot(ts,positions[:,1],label='first order integrator')). Note that x=positions[:,0] and y=positions[:,1].

• \(v_y\) vs \(t\)

• \(y\) vs. \(x\)

for both the exact and approximate curves. Use \(dt=0.05\) s.

You must label your axis (plt.xlabel("put actual label here...not what's written in this text"))

Put on each of your plots the exact answer computed using formulas you know from Physics 211 (i.e. plt.plot(ts,exact,label="Exact Answer"). Use plt.legend() to generate a legend so we know which line is which.

Your exact and integrator curves should look very similar but not exactly on top of each other (for the position).

Put the code which generates the plots (and the plots) here.

### ANSWER HERE

e. Animation

If you send animateMe the list of positions it will animate them. You shouldn’t have to change any code here as long as you’ve got a list of positions in the array positions. 🦉Go ahead and watch your thrown ball!

Use it as follows:

anim=animateMe(positions,True)
HTML(anim.to_jshtml())

Show code cell content

Hide code cell content

from matplotlib import animation
from IPython.display import HTML

def animateMe(positions,fullLine=True):
    # First set up the figure, the axis, and the plot element we want to animate
    fig = plt.figure();
    x_min=np.min(positions[:,0])
    x_max=np.max(positions[:,0])*1.1
    y_min=np.min(positions[:,1])
    y_max=np.max(positions[:,1])*1.1

    ax = plt.axes(xlim=(x_min, x_max), ylim=(y_min, y_max));
    #line, = ax.plot([], [], lw=2);
    line, = ax.plot([], [],'bo', ms=10);
    ax.set_xlabel("x")
    ax.set_ylabel("y")
    plt.close(fig)
    # initialization function: plot the background of each frame
    def init():
        line.set_data([], [])
        return line,

    # animation function.  This is called sequentially
    def animate(i):
        x = positions[:,0]
        y = positions[:,1]
        if fullLine:
            line.set_data(positions[0:i,0], positions[0:i,1])
        else:
            line.set_data(positions[i,0], positions[i,1])
#        line.set_markersize(10)
        return line,

    # call the animator.  blit=True means only re-draw the parts that have changed.
    anim = animation.FuncAnimation(fig, animate, init_func=init,
                                   frames=len(positions), interval=20, 
                                   blit=True,repeat=False);
    return anim

f. Errors

In this section, we’d like to understand the error of your integrator. For this problem, you can use the same initial conditions you’ve been using in the previous sections. We will assess the error in two ways.

First, how does the error accumulate as a function of time (when we fix \(dt\))?
We will only concern ourselves with the difference between the Euler-integrator \(y\) position and the exact \(y\) position (i.e. we will ignore the \(x\) position and velocity \(v\)). What you will find is that the error grows linearly with \(t\). This makes sense because each step makes an error and so cumulatively the total amount of error grows linearly. That means the longer you run your system, the larger error that you get. This linear growth will be true for almost any integrator that we use (although something special will go on for symplectic integrators).

🦉 Plot the difference \(y_\textrm{Euler}(t) - y_\textrm{exact}(t)\) as a function of time \(t\) for several values of dt (try at least 0.005 and 0.32). Check that the smaller value of dt is closer to the exact path.

Secondly, as \(dt \rightarrow 0\), we should obtain the same position at some time \(T\). That is, if we run a bunch of simulations with different \(dt\) to a total time \(T=1.28\) how does the final position change with \(dt\)?

🦉Write a function as follows:

def compute_position(dt, target_time):
  """
  For a fixed initial condition, run at timestep dt until the target time.
  Args: 
    dt: timestep (s)
    target_time:  time at which to return the position (s)
  Returns: 
    position at target_time
  """
  

Plot the \(x\) and \(y\) position for each value of \(dt \in \{0.32,0.16,0.08,0.04,0.02,0.01,0.005\}\)

While the error with time is always linear, the error with \(dt\) can change depending on the quality of your integrator. For the Euler-integrator you will find that the error with \(dt\) is also linear telling you that you are working with a first order integrator.

Put the two plots (and the code to generate them) here.

### ANSWER HERE

g. Energy

Another important way of assessing errors is to verify the laws of physics are followed in our simulation.

As we know from classical mechanics, the total energy is a constant of motion. Recall that the total energy is

\[ E = K + U = \frac{1}{2}mv^2 + mgy\]

where \(g=9.8\)

🦉Write a function that takes as arguments list of positions and velocities and returns a list of energy. You can assume that \(m=1.\) Then plot the energy vs time of the ball thrown in the air with your integrator. You should see that the energy is nearly constant over time. Because our approach is not perfect, there is a tiny bit of drift in the energy. Run it for a couple different timesteps and check that it gets better with smaller timestep.

Hint: the energy is around 50. Pay close attention to the scale on your plot; it may zoom in a lot.

### ANSWER HERE

---

---

---

Exercise 2: Air resistance

• List of collaborators:

• References you used in developing your code:

In this exercise we have two goals.

• One goal is to clean up our code so that it is more general (we will be using your 2d code). We have a lot of hard-coded things running around and so we want to get rid of that.

• Our second goal is to modify our code so that it can work with air resistance. We will be able to see how air resistance affects balls flying through the air. We can also see how a dropped ball reaches terminal velocity.

a. Code Modifications

In this section you are going to modify your code to compute dynamics with air resistance. The force should be of the form np.array([-b*v_x,-b*v_y-9.8*m]), where the constant b determines the strength of the air resistance. Because of the air resistance, the mass now matters, so that also needs to be an argument.

🦉Write a new step function to include air resistance.

def step_air_resistance_2d(t, pos, vel, dt, m, b):
  """
Starting with the position, velocity, and current time, return the next time, position, and velocity at time t+dt.
  The force applied is (-9.8 m/s^2)*m - b * v  
  Args:
     t: scalar, current time (s)
     pos: np.array of shape (2,), position at time t (m)
     vel: np.array of shape (2,), velocity at time t (m/s)
     dt: scalar, timestep (s)
     m: mass (kg)
     b: air resistance (kg /s)
  Returns
    tuple containing new time, new position, new velocity

  """
  assert isinstance(pos, np.ndarray), "pos must be a np.ndarray"

🦉Write a new throw_ball function.

def throw_ball_air_resistance(init_pos, init_vel, dt, m, b):
    """
    Run a simulation of a thrown ball until the ball hits the ground (at y=0).
    Arguments:
      init_pos: (2,) array given the initial positions.
      init_velocity: (2,) array giving the initial velocity.
      dt: scalar timestep
    Returns:
     (ts, vs, pos): (nsteps,), (nsteps, 2), (nsteps,2)
    """

Be careful to avoid global variables; your function must take in the parameters as arguments here.

Put your new code here. You don’t have to produce any output (you will produce graphs below)

resetMe()
### ANSWER HERE

b. Running and making Plots for Throwing a Ball

🦉For a ball in the air and make plots

• \(y\) vs \(t\)

• \(y\) vs \(x\) that includes (on the same plot) both the ball with air resistance and without air resistance. To label which is which you can do plt.plot(x,y,label="Air Resistance") and then call plt.legend(loc=0)

Use the following parameters:

bs = [0, 0.1] # you can loop over these
m = 1
initial_pos = np.array([0.0,0.0])
initial_velocity = np.array([0.1, 10.0])
dt = 0.01

Does it make sense? Under which conditions does the ball travel further?

### ANSWER HERE

c. Running and making plots for dropping a ball.

Now we would like you to set up a simulation which has a ball being dropped. The ball should have a mass of 1.0 kg and be dropped from a height of 1000 meters. The air resistance value of \(b\) should be 0.3. Use again a time-step of 0.01.

🦉We want you to plot

• \(y\) vs. \(t\)

• \(v_y\) vs \(t\).

We want you to see that the \(y\)-velocity reaches a terminal velocity (i.e. you stop picking up speed as the object is falling).

🦉Calculate by hand what the terminal velocity should be (consider when the force due to gravity and due to air reistance is the same) and draw a dotted line at this value on the plot (i.e. plt.axhline(terminalVelocity,linestyle='--').

Plots and graph here.

### ANSWER HERE

d. Hitting a target

Let’s suppose we want to throw a ball starting at position \((0,0)\) with some initial velocity v0 (which you will find) and have it hit a target on the ground at position \((10.0,0.0)\) m. Let \(b=3.0\) and \(m=1.0\). You may again use \(dt=0.01\). We are going to use scipy.optimize.minimize to find the initial velocity. The library scipy.optimize is very useful throughout this course.

To achieve this, first write a function def distance_from_target(init_vel,init_pos, dt, m, b): which takes an initial velocity and the parameters. The function should return how far (in absolute magnitude) a ball thrown with the parameters lands from the target.

Then you can call

ans=scipy.optimize.minimize(distance_from_target,
                            initial_velocity_guess, # guess for the starting velocity
                            args=(init_pos, dt, m, b),
                            method='Nelder-Mead',
                            options={'disp': True})

and the values in ans.x should be the velocity you found.

🦉Graph the curve you find to make sure it hits close to the target. Report your initial velocity.

Plot and graph here. You want to have print(ans.x) at the end so we can see what your results are.

### ANSWER HERE

e. Hitting a moving target (extra credit: 5 points)

Here we want to launch two balls. One has fixed parameters

b= 3.0
m = 1.0
initial_pos= np.array([0.0,0.0])
initial_velocity= np.array([10.0,10.0])
dt = 0.01,
target = 14.5

The other has the following parameters but you can optimize the initial velocity

b= 3.0
m = 1.0
initial_pos= np.array([10.0,2.0])
initial_velocity= np.array([-1.0,1.0])
dt = 0.01,
target = 14.5

🦉Figure out how to tune the initial velocity of the second ball so that the two balls collide in the air. Report your initial velocity and convince us that it actually worked.

Put your code and graph here. This is a little tricky since both balls are flying.

### ANSWER HERE

Exercise 3: Second Order Integrators

• List of collaborators:

• References you used in developing your code:

We saw that the first-order integrator we wrote to get to a time \(T\) had errors that scaled as \(dt\). In this exercise we are going to develop a integrator which scales as \(dt^2\).

In our first order technique we had used the position and the velocity at the start of an interval to estimate the force that would act during the interval. We then pretended that the force \(F\) and velocity \(v\) would remain constant during the interval so that we could calculate the increments to position and velocity as \(dx = v(t) dt\) and \(dv = [F(x,v,t)/m]dt.\) It was only at the end of each interval that we updated the position and the velocity by adding the increments to them.

We can improve the precision of our integrations by estimating the position, velocity, and force in the middle of our step (at time \(t+dt/2\)), then using that to decide how much the position and velocity will change during the entire time interval \(dt\). Let’s refer to this as a midpoint integration technique. You will likely want to pull your force calculation out into anotehr

If the position, velocity, and force at the start of the interval are \(x,v\), and \(F\), easonable estimates of the position and velocity at the midpoint of the time interval are

\[\begin{split} \begin{align} x_\textrm{midpoint} &= x + v \frac{\Delta t}{2}\\ v_\textrm{midpoint} &= v+ \frac{F(x,v,t)}{m} \frac{\Delta t}{2} \end{align}\end{split}\]

We can use these to estimate the changes \(\Delta x\) in position and \(\Delta v\) in velocity during the entire interval of duration \(\Delta t\), where \(t_\textrm{midpoint}\) is \(t + \Delta t/2\):

\[\begin{split} \begin{align} \Delta x &= v_\textrm{midpoint} \Delta t\\ \Delta v &= \frac{F(x_\textrm{midpoint},v_\textrm{midpoint},t_\textrm{midpoint})}{m} \Delta t \end{align} \end{split}\]

The approximate values of \(x\) and \(v\) at the end of the interval (and the start of the next) are

\[ x \leftarrow x+\Delta x\]

\[ v \leftarrow v+\Delta v\]

\[ t \leftarrow t+\Delta t\]

Note that \(x\), \(v\), and \(F\) are all vectors, even though we haven’t written the vector symbol over them

a. Midpoint Method

Implement this midpoint method replacing your previous step function. Test that it works by verifying you are essentially getting the same results as earlier (it might be slightly different because you have a smaller time step error).

Code for the midpoint method.

#ANSWER HERE

b. Measuring the error.

🦉Using your second order method, measure the time-step error of our integrator for a ball with friction in the same way that you did for the first order method.

Use initial conditions of \(x_0=[0,0],v_0=[0,10]\) m/s. Be sure to include air resistance, as without friction the second order method will be exact (so you won’t see a trend).

Unfortunately we don’t have a simple formula to get the exact result. What we will do instead is measure the answer for a time step that is very small (\(\Delta t=0.00001\)) and use that as our “exact” answer. Then we can compare it to larger time steps.

In the end, we want to determine \(\alpha\) in the equation \(\textrm{Error} \propto \Delta t^\alpha\) (the symbol \(\propto\) means “proportional to”). To get the \(\alpha\) we will use the process described in the warmup - i.e. plotting on a log-log plot and computing the slope.

Also note that you need to compare the positions at the same time, even though you are running at different timesteps. I’d suggest a time of 1.28, and one way to find the index is to do np.argmin(np.abs(ts-1.28)), where ts are the output of your integrator.

You should have your plots here. They should start with the dropping ball initial velocities. You should show that you are getting error which is second order.

### ANSWER HERE

Exercise 4: Higher Order Integrators (Extra Credit: 15 points total - 5 points/section)

• List of collaborators:

• References you used in developing your code:

a. Runge-Kutta (5 points)

We have implemented a first and second order integrator. We can even use higher order integrators. The typical one everyone uses is the fourth order Runge-Kutta (RK) integrator - see wikipedia.

🦉 Go ahead and implement a fourth-order RK integrator and show (using the same time-step techniques above) that it is fourth order.

b. scipy ODE (5 points)

🦉 Rewrite ThrowBall to use scipy.integrate.ode.

c. symplectic integrator (5 points)

Implement a symplectic integrator at any order. Show what order it works at and show that the energy stays close for a much longer period of time then Euler.

---

• Ex 1: George Gollin (original); Bryan Clark and Ryan Levy (modifications)

• Ex 2: Bryan Clark (original)

• Ex 3: George Gollin (original); Bryan Clark and Ryan Levy (modifications)

© Copyright 2021

from requests import get
from socket import gethostname, gethostbyname
ip = gethostbyname(gethostname()) 
filename = get(f"http://{ip}:9000/api/sessions").json()[0]['name']
filename = filename.replace('%20',' ')
print("CHECK THAT THE FILENAME IS correct", filename)
from google.colab import files
from google.colab import drive
drive.mount('/content/drive')
filepath = "/content/drive/MyDrive/Colab Notebooks/" + filename
!cp "$filepath" ./
!jupyter nbconvert --to HTML "$filename"
files.download(filename.replace("ipynb","html"))