ECMM171P Programming for Engineers
Assignment 3: Projectiles
Hand-in date: 18:00 Friday 13th December 2019
❼ This assignment will count for 50% of your total grade in the course.
❼ Remember that all work you submit should be your own work. This assignment is not
meant to be a team exercise. The University treats plagiarism very seriously.
❼ Ensure that your programs run without any error on the Harrison computers - Matlab
version R2019a. Programs with errors may get zero marks. Additionally, ensure that
you pay careful attention to the instructions and marking scheme below, particularly in
regards to defining specific functions, commenting your code appropriately and ensuring
its readability. Programs with no comments will get 20 % less than what would normally
be awarded.
❼ You should submit your source code by 18:00 Friday 13th December 2019
❼ If you have more general or administrative problems please e-mail me. Always include
the course number (ECMM171P) in the subject of your e-mail.
The path of a projectile is an interesting, important, and sometimes challenging problem. It’s
very likely that you will have studied the path of a projectile in the absence of atmospheric
drag in basic mechanics, since this can be derived analytically. In the real world, however,
drag can significantly alter the trajectory. Unfortunately, in introducing this effect, we lose the
ability to analytically obtain solutions for the projectile’s trajectory: simulating the trajectory is
therefore required. In this assignment, you will therefore write a Matlab program that computes
trajectories under the influence of drag, and use this to solve a simple targeting problem.
1 Theory
We will neglect any crosswind effects, so that the projectile travels in a two-dimensional plane
with coordinates (x, y). The projectile’s position, velocity and acceleration can therefore be represented
as vector functions of time, which we denote by r(t) = [rx(t), ry(t)], v(t) = [vx(t), vy(t)]
and a(t) respectively. As the projectile travels through space, we will assume that it experiences
two forces: acceleration due to gravity, and a drag force Fd acting in the direction opposite to
its trajectory. A free body diagram of this scenario can be found in Figure 1. The drag force
can be roughly approximated by the drag coefficient cd through the relationship
Fd =12ρcdA|v|v,
which also involves the velocity v, air density ρ and the frontal area A (i.e. the two-dimensional
cross-section exposed to the drag force). cd is dimensionless and is used to quantify the drag
of an object in a fluid: the lower the value, the lower the drag force.
To simulate the trajectory of the projectile, we can use Newton’s second law:
F = ma ⇒ a(t) = 1
Figure 1: Free body diagram of the projectile.
where m is the mass of the projectile and g = (0, g) is acceleration due to gravity with g =9.81 ms−2
. Since we are interested in the projectile’s trajectory r, we can then utilise the fact
that a(t) = dvdt= v˙(t) and v(t) = drdt= r˙(t)
to give us a ordinary differential equationr¨(t) = 1
(1)Notice that this is a non-linear equation that cannot be solved analytically – we therefore
need a computer to simulate it. It is also second-order: that is, it depends on the second
derivative of r.
All that is left to do is equip the above with an appropriate set of initial conditions at time
t = 0: we will imagine that the projectile is fired from an initial starting height h at an angle
of θ degrees and initial speed (velocity magnitude) v0, which gives:
r(0) = (0, h) and r˙(0) = v(0) = [v0 cos(θ), v0 sin(θ)].
2 Your tasks
Your program will be structured in two parts. In the first part, you will write a function to solve
the ODE system defined above. In the second part, you will use these functions to create a
simple targeting computer. In all of the below, you should assume that the projectile is a sphere
of diameter 0.1 m being fired through air of density ρ = 1.225 kg m−3
, with corresponding drag
coefficient cd = 0.479 and mass m = 5.5 kg.
2.1 Getting a trajectory
In the first part of this assignment, you will write a solver for the ODE in equation (1).
1. Using the techniques discussed in the lectures, first transform equation (1) as a system of
first-order equations for the two components of velocity and two components of position,
respectively. Mathematically, this system will be expressed as
y˙ = F(t, y(t)),
where y(t) = [vx(t), vy(t), rx(t), ry(t)]. Write a function drag ode, which takes as arguments
the time t and the vector y(t), and returns the function F(t, y) as a column
vector.
2.1 Getting a trajectory 2 YOUR TASKS
0 100 200 300 400 500 600 700 800 900
Downrange dustance [m]
Height [m]
Without drag
With drag
(a) Comparison of trajectories with and without drag
for initial conditions v0 = 150, h = 10 and α = 10.
0 200 400 600 800 1000
Downrange distance [m]
(b) Example of target being hit in target.m with
(tx, ty) = (1000, 50)
Figure 2: Graphs that your codes should produce.
2. Using the lecture notes, write a function called solve ode euler, which uses the forward
Euler time integration method and your drag ode function to solve the ODE in equation
(1). solve ode euler should take as arguments the initial speed v0, the angle of
attack θ, the initial height h and a timestep size ∆t. You should also include a maximum
integration time as a parameter. You should stop integrating once either the maximum
integration time has been reached or once the projectile reaches the ground.
The function should return: the solution time, the two components of v and the two
components of r. Include all outputs into a single variable with 5 columns. Each row
should contain these values at a given timestep (i.e. the first column is their values at
t = 0, the second at t = ∆t, and so on).
3. To validate your implementation, write a function called solve ode45, which instead
uses Matlab’s built-in function ode45 and your drag ode function to solve the ODE in
equation (1). solve ode45 should take the same arguments and return the same array as
solve ode euler.
4. You should include a script drag which will prompt the user for initial speed v0, the angle
of attack θ, the initial height h and will call the solve ode euler function to solve the
ODE system and produce a plot that compares this to the analytic solution of the particle
trajectory without drag. You may use the following formulae,
Distance travelled : d =v cos θgv sin θ +p(v sin θ)2 + 2gh
Height at position x : y = h + x tan θ −gx22(v cos θ)2
Travel time : t =dv cos θ
Think carefully about the choice of ∆t. Your script should output a graph that shows
the two trajectories with and without drag, as in Figure 2, and will print: initial speed,
angle of attack, initial height, travel time, distance travelled and maximum height of the
projectile.
2.2 A targeting computer 2 YOUR TASKS
2.2 A targeting computer
In your second part, you should use the Matlab’s optimisation package to build a simple targeting
system. The idea here is that you are given the position (tx, ty) of a target which is
downrange from your projectile. Given the constraints that:
❼ the projectile’s initial velocity can be between 100ms−1 and 800ms−1
;
❼ the projectile’s initial angle can be between 10◦ and 80◦
.
your program should determine what the initial velocity and angle of the projectile should be
in order to hit this target. In all of the below, you can assume that the initial height is h = 0.
1. Write a function called objective, which defines the objective function to be minimised.
objective should take as arguments: the two variables being optimised (i.e. the initial
velocity and angle), and the target coordinates tx and ty. You may combine the arguments
into two variables for convenience. It should call solve ode45 with the initial conditions
and return the shortest distance between the computed trajectory and (tx, ty).
2. Write a function called objective image. The purpose of this function will be to produce
an image, along the same lines as discussed in lectures, that shows the value of the
objective function across the range of constraints, as shown in Figure 3. objective image
should take as arguments the target (preferably in a single variable as in objective), and
two integers nx and ny, which are the number of pixels to be used to represent the range
of initial velocities (in the x-direction of the plot) and the initial angle (in the y-value of
the plot). Your function should then plot the objective function across these ranges and
produce an output file called objective.pdf.
3. Write a function called target, which takes as two arguments the target’s position tx and
ty. Use routines from Matlab’s optimisation toolbox to minimise objective, and return
True if the optimisation succeeds, along with v0, θ and the minimum distance; otherwise
return False and three zeros.
Some hints:
❼ This aspect is the most open-ended in the assignment. The function can be quite
hard to optimise and is sensitive to initial guesses.
❼ Use your images from the objective image as a guide: think about how the optimisation
(minimisation of the objective function) process works and how it might
work for some of these images to guide your choices of e.g. initial guesses.
❼ You can use any combination of functions from the optimisation toolbox, but you
may explore fmincon and ga in particular: note the former is much faster than the
latter, but is more prone to becoming stuck during optimization.
❼ You should consider in particular the bounds and args options in these functions
and you can consult Matlab’s documentation to determine appropriate arguments
to call.
Note that there is no unique solution to this problem: full marks will be awarded for
creativity and devising a relatively robust optimisation strategy that works for most cases,
which demonstrates knowledge of the programming and numerical techniques under the
hood.
4. Finally, write a script called test target that prompts the user for the coordinates (tx, ty),
calls your target function, and prints a message target was hit! along with some
3 MARKING SCHEME
Distance to target [m]
100 200 300 400 500 600 700 800
Initial velocity [m/s]
Figure 3: Image of the objective function, which shows the distance between computed trajectories
and a target (tx, ty) = (1000, 50) for each initial velocity and angle.
information if the target is hit: initial speed, and initial angle. The script should also
produce a plot which shows the target and the optimised trajectory, an example of which
can be seen in Figure 2(b).
5. On the ELE page, you will find a randomly-generated target to hit for your student ID.
You should run your code for this target and submit the objective function image from
objective image along with all your programs.
3 Marking scheme
Your submission will be assessed on the following criteria:
❼ For the drag part:
– Fully working implementation of the drag ode function. [5%]
– Fully working implementation of the solve ode euler method. [20%]
– Fully working implementation of the solve ode45 method. [10%]
– Fully working script
drag and correctly formatted graph. [10%]
❼ For the target part:
– Fully working implementation of the objective function. [5%]
– Fully working implementation of the objective image function. [10%]
– Fully working implementation of the target function. [15%]
– Fully working script test target which produces correctly formatted graphs of objective
function and trajectory. [10%]
4 SUBMISSION
❼ Discretionary marks for good programming technique, structure and commenting. [15%]
Details regarding the specific criteria that are used in the above can be found in the module
handbook on ELE. You should attempt all of the above routines; incomplete solutions will be
given partial credit. Note that these routines will all be marked independently, and you can still
receive full credit for complete solutions even if the preceding parts are not fully implemented.
4 Submission
You should package all your files along with the .pdf graph, into a single zip, rar or tar file,
and submit using the online Turnitin link on the ELE page.