Simulation of Ballistic Motion, Interpolation, Runge-Kutta, and Solve_ivp

B. Stratospheric skydive

Felix Baumgartner (Links to an external site.), and later Alan Eustace (Links to an external site.) have jumped from helium balloons at altitudes around 40km, free-falling to Earth at record-breaking speeds. At this altitude, gravity is still approximately the same as at sea level, but atmospheric drag is much reduced. Aerodynamic drag is modelled as where is the drag coefficient which we (crudely, but reasonably accurately) assume is constant. is the density (a function of altitude) and is the speed. 

 

Write a function to simulate the 1D dynamics of ballistic motion through the atmosphere under constant gravity . The simulation should terminate when the object hits the ground. Use this function to simulate a fall from 50km of an object with mass 100kg, crosssectional area 0.5m², and drag coefficient 0.5.

Plot a graph of speed vs time. Label your graph using plt.xlabel and plt.ylabel. Using plt.title, set the title of the graph to include the maximum speed attained during descent. (This value should be extracted using code, and not entered by hand.)

The Quadrupole ion trap (Links to an external site.) is a crucial technology from mass spectroscopy to quantum computing. Using oscillating electric fields, it confines charged particles in vacuum. A rotating saddle with a ball placed on top has long been used as a mechanical analog
of this trap. Video illustration (Links to an external site.) The saddle is described by where and are length scales. The acceleration for a ball resting on this surface is . Hence, acceleration as a function of position for a saddle rotating at angular frequency is where are coordinates i.e. the vector acceleration depends on position. For this question, we ignore the rotational inertia of the ball'

Simulate 10s of the dynamics for a ball released at (5cm,0cm), initally at rest, for rotational frequency 1Hz (i.e. ) and A=10cm. Use L=1m throughout. Plot the trajectory in the plane with plt.plot(xs,ys) (or similar)and use plt.axis('equal') to control the aspect ratio. Interpolating

Say we have some points on a curve and we want a model-independent way to approximate values in between. For example, we might have some altitude measurements on a hillside, and we want to estimate the altitude at intermediate points. This is called interpolation. Start with 1D. Imagine we have some points on a curve.
x = [-2. , -1.6, -1.2, -0.8, -0.4, 0. , 0.4, 0.8, 1.2, 1.6, 2. ]
y = [0.004, 0.100, 0.226, 0.566, 0.853, 0.965, 0.855, 0.539, 0.235, 0.095, 0.024]   import matplotlib.pyplot as plt  plt.plot(x,y,'.') # '.' means plot points. see plt.plot? for other formatting option.

C: Rotating saddle trap

We can interpolate to guess, in some reasonable way, the values between these known data points. from scipy.interpolate import interp1d
f = interp1d(x,y) # basic form. doesn't allow points outside. defaults to linear in terpolation

f_lin = interp1d(x,y,kind='linear', bounds_error=False, fill_value=0)

f_cub = interp1d(x,y,kind='cubic', bounds_error=False, fill_value=0)

x_int = np.linspace(-3,3,1001)

y_lin = f_lin(x_int)

y_cub = f_cub(x_int)

plt.figure()

plt.plot(x_int,y_lin,label='linear')

plt.plot(x_int,y_cub,label='cubic')

plt.plot(x,y,'.',label='data')

plt.legend()

ODEs and methods

def rk4(f, x0, dt):  """Solve differential equation x'(t) = f(t,x(t))  using 4th order Runge--Kutta with timestep dt. Returns a generator which yields each tn,xn.
 tn = 0
 xn = x0

 while True:
 # gives x0 first, then calculates xn for n>0
 yield tn, xn

 k1 = dt*f(tn,xn)
 k2 = dt*f(tn+dt/2,xn+k1/2)
 k3 = dt*f(tn+dt/2,xn+k2/2)
 k4 = dt*f(tn+dt,xn+k3)

 # use xn = xn + ... (not x += ...) so a numpy array is copied, not changed in
-place
 xn = xn + (k1+2*k2+2*k3+k4)/6
 tn = tn + dt

(use the code of Runge–Kutta above and then use the code below) def f(t,X): # function to compute rate-of-change
 x,v = X # for the damped harmonic oscillator xdot = v
 vdot = -a*v-b*x
 Xdot = np.array([xdot,vdot]) return Xdot a, b = 0.1, 1 # parameters of the damped oscillator
X0 = np.array([1,0]) # initial conditions: displaced, at rest solver = euler(f,X0,0.01) # glue it all togther (again, never use Euler for a real problem)
 ts = []
 Xs = [] for t,X in solver:  ts.append(t)
  Xs.append(X)
  if t > 100:  break Xs = np.array(Xs) # convert the list of state vectors into a 2D array import matplotlib.pyplot as plot plot.plot(ts, Xs[:,0]) # plot position vs time solve_ivp

Example

 

Here’s how to solve the undamped harmonic oscillator using solve_ivp. import numpy as np import matplotlib.pyplot as plt from scipy.integrate import solve_ivp # 1
def f(t, X):
x,v = X

xdot = v
vdot = -x  return np.array([xdot, vdot])
 X0 = np.array([1,0])
solution = solve_ivp(f, [0, 100], X0) # 2
ts = solution.t # 3
Xs = solution.y.T # 4
plt.plot(ts, np.cos(ts), label='Analytical')
plt.plot(ts, Xs[:,0], label='solve_ivp')
plt.legend(); plt.xlabel('time')Much of this is very familiar.

The new bits are numbered above:

1. Import the function
2. Compute the solution for time from [0, 100] and store in solution
3. The times at which solutions are computed are stored in solution.t
4. The values of the state-space vector at these times is stored here; .T takes the transpose Specifing times at which to compute solution When we plot this, we notice that it’s correct (i.e. agrees with analytical, as far as we can see) but it’s rather sparse (i.e. we’d like to know the points a bit more frequently). We can either 1) specify the times at which to find solutions, or 2) ask it for “dense” output.

1. We make an array using linspace and ask for solutions at these times: ts = np.linspace(0,100,10001) solution = solve_ivp(f, [min(ts), max(ts)], X0, t_eval=ts)
2. We pass the optional argument dense_output solution = solve_ivp(f, [0, 100], X0, dense_output=True) Specifying when to stop the simulation
solve_ivp has an optional argument events. One use of this is to stop the simulation when some condition is met. The output from the function should cross zero when the condition is met. An example of this is below, where a projectile will be solved until it hits the ground.
def trajectory(speed, angle, drag): from numpy import array, cos, sin from scipy.integrate import solve_ivp
g = 9.81

def f(t,X):
x,y,vx,vy = X
xdot,ydot = vx,vy
vxdot,vydot = -drag * vx, -drag * vy - g
return array([xdot, ydot, vxdot, vydot])
def hit_ground(t,X): # 1
x,y,vx,vy = X # 2
return y # 3
hit_ground.terminal = True # 4
hit_ground.direction = -1 # 5
X0 = [0,0,speed*cos(angle),speed*sin(angle)]
solution = solve_ivp(f, [0, 100], X0, events=hit_ground, max_step=0.01) # 6

return solution.t, solution.

 

1. We define a function which will test whether a condition is met; it has the same signature (t,X) as our function f.
2. We unpack values from the state space in the familiar way.
3. In this case we simply return the vertical position.
4. This is new. We are directly setting a property of the function. solve_ivp will look for this.
5. Similarly, we set a property to tell solve_ivp that this condition means “STOP!” when the particle is coming downwards.
6. We pass this function to solve_ivp as the events argument