Homework 1: Computing and Data Analysis Review#
Learning Objectives#
In this homework assignment, you will review and practice the following skills from CBE 20258:
Plotting in Python including best practices
Nonlinear regression with weights based on measurement uncertainty
Error propagation coupled with nonlinear regression
Solving a boundary value problem by calculating the root of a nonlinear function
Each topic links to pages on a website with the CBE 20258 content. The search feature on the CBE 20258 is a great way to find more related examples.
Introduction#
This homework problem is inspired, in part, by a Wired blog post analyzing an incident with a cannonball on the popular show Myth Busters.
Consider \(n\) iron cannonballs of mass 12 lb (5.44 kg) fired at unknown angles \(\theta_i\) (where \(i \in \{1,...,n\}\)) with unknown initial velocities \(v_{0,i}\). The drag on the cannonball is \(F_D = \frac{1}{2} \rho C_D A v |v|\) where \(\rho\) is the density of air, \(C_D\) is the coefficient of drag, and \(A\) is the cross sectional area of the cannonball. We can lump these parameters together into a single unknown coefficient \(\bar{C}_D\):
It is important to use \(v |v|\) to account for the change in sign when the cannonball transitions from moving upward to moving downward.
Applying Newton’s laws of motion gives the following coupled differential equations:
We assume the cannonball is fired from the origin at time \(t=0\) where \(x\) and \(z\) represent horiztonal and vertical displacement, respectively. This gives the follow initial values:
The \(n\) cannonballs are fired on completely flat terrian. We observe the impact location (horizonal distance traveled), impact time, and crater size. The latter can be correlated with impact velocity:
Time (s) |
Horizontal Displacement (m) |
Impact Velocity (m/s) |
|---|---|---|
12.0 |
726 |
73 |
6.8 |
529 |
77 |
15.3 |
616 |
76 |
Given these data, our ultimate goals are:
Build a mathematical model
Estimate \(\theta_i\), \(v_{0,i}\) and \(\bar{C}_D\)
Predict the distance traveled and impact velocity of a cannonball shot at a new angle and initial velocity with uncertainty
We will assume the time, location, and velocity measurements are corrupted with normally distributed noise with a standard deviations of 0.2 s, 5 m, and 20 m/s, respectively. When we do regression, we need to account for the fact that the velocity estimates (based on crater size) have greater uncertainty.
Exercise 1#
On paper, derive the above ordinary differential equation model from first principles, similar to Physics I. Start with a free body diagram. Use trigonomety to justify the initial conditions.
Exercise 2#
Write a Python function to numerical simulate a cannon ball experiment. The inputs to your function should be the unknown parameters \(\theta_i\), \(v_{0,i}\), and \(\bar{C}_D\). Your function should return:
Final time \(t_f > 0\) when the cannonball returns to zero elevation, i.e., \(z(t_f) = 0\)
Horizontal distance traveled, \(x(t_f)\)
Impact velocity, \(v(t_f) = \sqrt{v_x(t_f)^2 + v_z(t_f)^2}\)
import numpy as np
import matplotlib.pyplot as plt
import scipy.integrate as integrate
%matplotlib inline
def ode_rhs(t, y, Cd):
"""ODE system for a cannonball with drag
Arguments:
t: time (s)
y: state vector [x, z, vx, vz] (m, m, m/s, m/s)
Cd: drag coefficient (kg/m)
Returns:
dy: time derivative of the state vector [vx, vz, ax, az] (m/s, m/s, m/s^2, m/s^2)
"""
x, z, vx, vz = y
vabs = np.sqrt(vx**2 + vz**2)
g = 9.81 # m/s^2
m = 5.44 # kg
# Add your solution here
def simulate_experiment(Cd, angle, v0=100, tmax=200, plot=True, tdata=None, xdata=None, vdata=None):
"""Simulate an experiment.
Arguments:
Cd: drag coefficient (kg/m)
angle: launch angle (degrees)
v0: initial speed (m/s)
tmax: maximum time (s)
plot: bool, if true, plot simulation results and data
tdata: measured impact time (s) [only used for plotting]
xdata: measured impact location (m) [only used for plotting]
vdata: inferred impect velocity (m/s) [only used for plotting]
Returns:
tfinal: final time (s)
xfinal: final x position (m)
vfinal: final speed (m/s)
"""
# Add your solution here
Exercise 3#
Test your function by simulating the provided data for the cannonball experiments. For each experiment, make a plot showing the model predictions (lines) and observed data (symbols). You can use the data given below as intial guesses for \(\theta_{i}\) and \(v_{0,i}\). You can use \(\bar{C_D} = 0.003\) as an intial guess for the coefficient of drag.
import pandas as pd
data = {}
data['time'] = [12.0, 6.9, 15.3]
data['distance'] = [726., 529., 616.]
data['velocity'] = [73., 77., 76.]
data['theta_guess'] = [40., 20., 60.]
data['v0_guess'] = [110., 100., 90.]
df = pd.DataFrame(data)
df.head()
# Confirm your plots show both the simulation results (lines)
# and experimental data (symbols)
for i, row in df.iterrows():
simulate_experiment(0.003,
row.theta_guess,
xdata=row.distance,
vdata=row.velocity,
tdata=row.time,
v0=row.v0_guess,
plot=True)
Exercise 4#
Perform nonlinear regression to estimate \(\bar{C}_D\), \(\theta_1\), …, \(\theta_n\), \(v_{0,1}\), …, \(v_{0,n}\). Report the estimate parameters with a reasonable number of significant digits. Include plots to compare the model predictions and experiment data.
from scipy.optimize import least_squares
def unpack(param):
""" Unpacks the vector "param"
Arguments:
param: vector [Cd, theta1, theta2, ..., thetaN, v01, v02, ..., v0N]
Returns,
Cd: lumped drag coefficient (kg/m)
theta: launch angles (degrees)
v0: initial speeds (m/s)
"""
n = (len(param) - 1 )//2
Cd = param[0]
theta = param[1:(n+1)]
v0 = param[(n+1):(2*n+1)]
return Cd, theta, v0
def pack(Cd, theta, v0):
""" Packs the parameters into a vector
Arguments:
Cd: lumped drag coefficient (kg/m)
theta: launch angles (degrees)
v0: initial speeds (m/s)
Returns:
param: vector [Cd, theta1, theta2, ..., thetaN, v01, v02, ..., v0N]
"""
return np.hstack([Cd, theta, v0])
# Assumed standard deviation of measurement error
sigma_t = 0.2 # s
sigma_x = 5 # m
sigma_v = 2 # m/s
def residuals_all(parameter_vec, data, plot=False):
"""
Computes the residuals for all experiments
Arguments:
parameter_vec: vector [Cd, theta1, theta2, ..., thetaN, v01, v02, ..., v0N]
data: dataframe with experimental data
plot: bool, if true, plot simulation results and data
Returns:
residuals: vector of residuals scaled by the measurement uncertainty
"""
n = len(data)
Cd, theta, v0 = unpack(parameter_vec)
# Allocate matrix for residuals, n experiments (rows) with 3 measurements per experiment (columns)
residuals = np.zeros((n,3))
# Loop over experiments
for i, row in data.iterrows():
# Simulate experiment
tsim, xsim, vsim = simulate_experiment(Cd,
theta[i],
v0 = v0[i],
xdata=row.distance,
vdata=row.velocity,
tdata=row.time,
plot=plot)
# Compute and record scaled residuals
residuals[i,0] = (tsim - row.time)/sigma_t
residuals[i,1] = (xsim - row.distance)/sigma_x
residuals[i,2] = (vsim - row.velocity)/sigma_v
# Flatten the residual matrix into a vector
return residuals.reshape(3*n)
# Define initial guess
initial_vector = pack([0.003], df.theta_guess.to_numpy(), df.v0_guess.to_numpy())
# Evaluate residuals at initial guess
print(residuals_all(initial_vector,df,plot=False))
# Perform nonlinear regression
# Add your solution here
# Determine the number of experiments
n = len(df)
# Extract the estimated parameters
Cd_est = sln.x[0]
theta_est = sln.x[1:(n+1)]
v0_est = sln.x[(n+1):(2*n+1)]
# Print the results (with rounding and units)
print("Cd =",round(Cd_est,4)," kg/m")
for i in range(0,n):
print("\nExperiment",i+1)
print("theta =",round(theta_est[i],1),"degrees")
print("v0 =",round(v0_est[i],1),"m/s")
# Plot the experimental data and the simulation results
residuals_all(sln.x, df, True)
Exercise 5#
Quantify the uncertainty in your parameter estimates. Hint: Calculate the covariance matrix and correlation matrix of the estimated parameters.
# Add your solution here
# Add your solution here
# Add your solution here
# Add your solution here
# Add your solution here
Interpret your uncertianty estimate in a few sentences.
Answer:
In a few sentences, describe and justify your choice of uncertainty quantification method. In other words, what are the pros and cons of the selected method and why it is a reasonable choice for this problem.
Answer:
Exercise 6#
Determine the firing angle \(\theta_j\) and intial velocity \(v_{0,j}\) to hit a target at 300 m with an impact velocity of 80 m/s.
from scipy.optimize import fsolve
def firing_residuals(guess):
""" Residuals for the firing experiment
Arguments:
guess: vector [angle in deg, v0 in m/s]
Returns:
residuals: vector [position, impact speed]
"""
# Add your solution here
Exercise 7#
Quantify and interpret the uncertainty in your prediction (distance, time, velocity) from the previous exercise. Assume you can specify the firing angle and initial velocity with uncertainties of 1 degree and 2 m/s, respectively. You also have some uncertainty in your estimate of \(C_D\) (quantified above). Hint: Calculate the covariance matrix and correlation matrix for the predicted displacement and impact velocity.
# Add your solution here
# Here is code copied from CBE 20258 notebooks that you may find helpful
# You can also answer this question without calculating a Jacobian
def Jacobian(f,x,delta = 1.0e-7):
'''Approximate Jacobian using forward finite difference
Args:
f: vector-valued function
x: point to build approximation J(x) around
delta: finite difference step size
Returns:
J: m x n Jacobian matrix (approximation), where
n = length of x
m = length of f(x)
Reference:
https://ndcbe.github.io/data-and-computing/notebooks/06/Newton-Raphson-Methods-for-Systems-of-Equations.html
'''
#Evaluate function f at x
fx = f(x) #only need to evaluate this once
# Determine size
N = x.size
M = fx.size
# Allocate empty matrix
J = np.zeros((M,N))
idelta = 1.0/delta #division is slower than multiplication
x_perturbed = x.copy() #copy x to add delta
# Loop over elements of x and columns of J
for i in range(N):
# Perturb (apply step) to element i of x
x_perturbed[i] += delta
# Approximate column in Jacobian
col = (f(x_perturbed) - fx) * idelta
# Reset element of x
x_perturbed[i] = x[i]
# Save results
J[:,i] = col
# end for loop
return J
# Add your solution here
Please interpret the uncertainty in your predictions. How do these results compare to the measurement uncertainty?
Answer:
Submission Instructions#
Please submit your answer to Exercise 1, neatly written on paper, to Gradescope as a PDF. Please submit your answers to the remaining exercises in this notebook (.ipynb file) to Gradescope. There are two separate assignments on Gradescope.