Newtonian mechanics and equations of motion

← Previous

Homework assignment

This homework builds on the skeleton code to be found here:

• The file Cannonball.py contains the steering code and includes the parts specific to the physics problem;
• the file DEq_Solver.py is the file containing the solution algorithms for differential equations, which you worked on in the previous homework. On Thursday 28.1., after the submission deadline for this homework, we will publish this file as it should be with the Euler method successfully implemented;
• the file Interpolator.py, which we will use to implement methods for interpolation and root finding.

The assignment consists of three parts:

1. Implementing and enhancing the physics:
   This part relies on the Euler method already discussed and implemented in the homework and the lectures before. In terms of implementation it reduces to corresponding changes in methods in the physics class Cannonball.

   1. Transform the set of equations

      \[ \begin{eqnarray*} \dot{x}\,&=&\,\frac{\mathrm{d}x}{\mathrm{d}t}&=& v_x \\ \dot{y}\,&=&\,\frac{\mathrm{d}y}{\mathrm{d}y}&=& v_y \\ \dot{v_x}\,&=&\,\frac{\mathrm{d}v_x}{\mathrm{d}t}&=& -\frac{B_\mathrm{drag}vv_x}{m}\\ \dot{v_y}\,&=&\,\frac{\mathrm{d}v_y}{\mathrm{d}t}&=& -g-\frac{B_\mathrm{drag}vv_y}{m} \end{eqnarray*} \] describing the physical problem such that it can be written in the form \[\begin{eqnarray*} \underline x_{i+1} &=& \underline x_i+\Delta\underline x\\ t_{i+1} &=& t_i + \Delta t\,, \end{eqnarray*}\] where, as usual, $\Delta\underline x$ depends on $\underline f(\underline x_i,t_i)$, and $\underline x = (x,y,v_x,v_y)$. As before, implement the function $\underline f$ in the "dx_dt" method of your physics class, this time called "Cannonball".

   2. For realistic, long-range artillery, a significant effect consists of the inclusion of the effect of the air density on the drag effect, i.e. the resistance the air offers to the projectile. It is well-known that the density changes with altitude, making it a function of $y$ in our problem. The simplest approximation is to treat the atmosphere as an isothermal ideal gas, assuming that its temperature is constant. The pressure $p$ then depends on the altitude,

      \[ p(y)\,=\,p(0)\,\exp\left[- \frac{\mu gy}{k_BT}\right] \, , \] where $p(0)$ is the pressure at sea-level, where we placed the canon. $\mu$ is the average mass of the molecules in the atmosphere, $k_B$ is Boltzmann's constant, and $T$ is the atmosphere's absolute temperature. Since for an ideal gas pressure and density are proportional this implies that the density reads \[ \rho(y)\,=\, \rho(0)\,\exp\left[- \frac{y}{y_0}\right]\,\, , \,\,\,\,\, y_0 \,=\, \frac{k_BT}{\mu g} \approx 1.0\cdot10^4 \mathrm{m}\, . \] It is clear that the isothermal approximation is not very realistic, since we know that temperature may vary quite significantly with altitude. Assuming that the air is a poor heat conductor and that therefore convection is very slow leads to the adiabatic approximation, which ultimately leads to \[ \rho (y)\,=\, \rho(0) \left(1-\frac{ay}{T_0}\right)^\alpha \, ,\] where the parameters $a$ and $\alpha$ have to be taken from data, and where $T_0$ is the sea-level temperature.
      In the code you may fix $a=0.0065$ K/m, $\alpha=2.5$, and $T_0=300$ K.

      It is clear that with sufficiently large $y$ differences between these two approximations become apparent.
      However, in both cases the drag force is proportional to the air density, leading effectively to a drag force that depends not only on velocity but also on the altitude:

      \[{\mathrm{\underline F}}_\mathrm{drag}(y,{\mathrm{\underline v}}) \,=\, \frac{\rho (y)}{\rho_0}\,{\mathrm{\underline F}}_\mathrm{drag} ({\mathrm{\underline v}})\,=\, -\frac{\rho(y)}{\rho_0}\,B_\mathrm{drag}v^2{\mathrm{\underline e}}_v \, . \] So, obviously, the original force is to be modified by the extra factor $\rho(y)/\rho_0$.

      To quantify these differences, add both versions of the correction factor to the class "Cannonball". It is a good idea to add them in such a way that the original version of the drag force is multiplied with the result of a corresponding method, and to have two different methods. This has already been done, with obvious empty methods as place holders. To make this accessible for marking an input choice of altitude dependence after the input of the drag coefficient such that the answer "const" corresponds to no altitude dependence, "isothermal" refers to the isothermal approximation, and "adiabatic" refers to the adiabatic approximation. You'll find that this choice is already implemented as a place holder in the program from the lecture.

2. Implementation of the Runge-Kutta method:
   This part adds better methods to the DEq_Solver class.

   1. Start from the DEq_Solver class for the evaluation of first order differential equations and supplement methods for solutions employing the Runge-Kutta methods of order 2 and 4. To this end, use the equations in the section on Runge-Kutta methods, keeping in mind that the x now are vectors, $x\to\underline x$.

   2. Check that the results obtained with both methods converge to the same trajectory by decreasing the step size from 0.1 down to 0.001.
      To facilitate this a switch between the methods is provided selecting "RK2" for the second-order Runge-Kutta, "RK4" for the fourth-order Runge-Kutta, and "Euler" for the Euler method, respectively.

   3. You will need to modify the "solve" method in the DEq_Solver class to accomodate a condition to stop the stepping not in $t$ but in $\underline x$, in particular its $y$-component. The solve function accepts an optional argument, "positive_index" to indicate which component of $\underline x$ should be checked. For example, with "positive_index=3" the while-loop should interrupt as soon as $\underline{x}[3]$ turns negative. This first negative point should still be recorded.

3. Range of shot:
   

   1. Fill the interpolate method in Interpolator.py. It takes a value $x$, a left pair $x_i,f(x_i)$ and a right pair $x_{i+1},f(x_{i+1})$, where $x_i < x < x_{i+1}$. The function should return the linearly interpolated value $f(x)$ between these two points.

   2. Fill the $f$ function in Interpolator.py, The Interpolator class takes two same-length lists in the constructor, one with the $x_i$ values and one with the matching $f(x_i)$ values to be interpolated. The function should return $f(x)$ over the full range of these input lists.

   3. Use the interpolator to implement the range() function in Cannonball. It should return the point along x where the height becomes zero again.