CSC173: ODEs

Introduction

This assignment was originally for an introductory programming class. Feel free to ignore all the code-structuring instructions. If they're useful, fine. I've tried to mark this introductory material "Optional" below.

As usual, this is an open-ended project. Use your instincts and judgement to choose and ignore options with the goal of maximizing results while minimizing effort.

Do either:

  1. Ballistics: "Balls in the Air".
  2. "Balls to the Wall" is lots of fun but also lots of code.

The expectation is that you'll be using various sources for your investigations, so make sure you cite all references you used: the more (that you actually use) the better. Technique, results, and presentation (especially graphics) are all important. The elegance and organization (and correctness) of the code count, but not as much as the report, which should contain enough information to convince the reader that your code must be correct.

In all this work, you will use Matlab's ODE solvers: let's just stick to ode23. The ODE solution itself is a simple and well-defined problem: one 'derivative' function just encodes the problem's differential equations (often written out for you in mathematics). The serious computing is done by a 1-liner: calling ode23 with a few simple arguments. These ballistics problems call for various amounts of outside research (e.g. just how does atmospheric density vary with altitude?), and a little interesting programming (finding the range and range-maximizing elevation with no user interaction).

At least one good previous report on the Balls in the Air project is available on the 160 website .../u/brown/160/.

Going Ballistic: Writeup 60%, Code 40%

Ballistics Basics

There are plenty of problems you can solve in closed form concerning the trajectory of bodies in a vacuum: Like These. But we normally fight in the atmosphere, and it turns out that air friction at muzzle-velocity speeds goes as the square of the speed of the the projectile (unlike the linear friction we saw in the mass-spring-damper example). It's also proportional to the cross-sectional area of the ball and the density of air. We get two or three coupled non-linear (because of the square of the velocity) second-order equations (for the vertical and horizontal, or also sideways, component to the ball's flight.)

State: Our first problem is to write the derivative functions for our problem variants. Since we want range and height of the cannonball (and later, its right-to-left deviation), let's say x, y, and z are the range, height, and left-deviation coordinates of the ball, respectively. Then x' (called u below), y' (called v), and z' (w) are the x, y, and z velocity components, respectively. The physics (coming up) will give us x'', y'', and z''. That means our state vector will consist of (x, x', y, y', z, z') (in some order). The derivative function takes in the state vector and returns its derivative (it copies x', y', z' into the output and uses the physics to compute x'', y'', z'').

Physics: Let s be the ball's speed through the air; then
s2 = (u2 + v2 + w2).
A is the cross-sectional area of the ball and m its mass, ρ is the density of air (zero for a vacuum), g the force of gravity, and c is a dimensionless drag coefficient. For spheres, c tends to be about .5. Technically, c is the ratio between the actual drag force, and the theoretical drag force that would result from a hypothetical object that brought all the air it intercepted up to its own speed. Unstreamlined objects like cubes, flat surfaces normal to the airflow, and end-on cylinders tend to have drag coefficients pretty close to 1.0. Nice, streamlined, torpedo-shaped objects can have drag coefficients of 0.1 or even lower.

The force due to air resistance for a symmetrical object is directly against the direction of motion (that is, in the direction of the apparent "wind"), and has magnitude given approximately by
F = (c A ρ s2) / 2.
Since we are formulating the problem in terms of x, y, and z components, we have to determine the force components in those directions. In particular, we want the force in each direction as a function of the x, y, and z velocity components (x', y' and z').

In this section we assume the medium (air) is at rest, but below that assumption changes. For now let's stick with the "still air" assumption.

The unit vector in the direction of instantaneous velocity is
(α, β, γ) = (1/√s2) (u,v,w).
Its three components give the proportion of the force to be allocated to the (x,y,z) directions.

From " force equals mass times acceleration", or F = ma, we get a = F/m, thence
x'' = - (1/(2m)) (ρ c A s2 α) (drag alone)
y'' = -g - (1/(2m)) (ρ c A s2β) (gravity plus drag)

Express the area in terms of the diameter and pi:

A = π (d/2)2

Make all these substitutions and we have all we need to write the derivative function that computes x'' and y'' from x, x', y, y', and our physical constants. We'll deal with z later.

What Now? The cannon points up at some elevation angle θ. Let's say the cannon points out the x axis, y is up, z is left. Cannon has some muzzle (initial) velocity, which should be subsonic for these equations: let's standardize on 150 m/sec, which is just under 1/2 the speed of sound. Given the muzzle velocity and your chosen shooting elevation angle θ you can calculate the initial velocities x', y',z' at time zero (that is, u0, v0, w0). (hint: z' =0, and sin θ and cos θ are involved in the calculations for u0, v0). The initial [0 0 0] position and [u0 v0 0] velocity conditions are sent in as the two columns of the 3x2 initial-conditions matrix to ode23() (see below).

For the problems below unless otherwise directed, use these base specifications in SI units (used to be called MKS for meter-kg-sec): muzzle velocity = 150 m/sec, mass of cannonball = 5kg, ρ = 1.23 kg/m3, g = 9.81 m/sec/sec, c = .5. For an iron cannonball of mass 5kg, its diameter (we need that to compute A) follows from the density of iron: 7860 kg/m3.

We're going to use Matlab's ode23(), and it wants the handle of the derivative function, a time span, and an initial state vector (let's say the cannon is at x = y = z = 0). Some experimentation is needed on the timespan, since our equations don't know anything about the ground and we should stop our solution about when the ball hits the ground.

With constants close to those above and an elevation of 45 degrees I get a trajectory that looks like this: plotting height (in m) against range (in m). Notice the slight deviation from a pure parabola. In fact, you can just set air density (ρ) to 0 and run that case with
>> hold on
to superimpose the plots to visualize the difference. Soon you'll be asked to see how elevation, drag, and mass affect the trajectory. You'd expect a lighter, projectile to slow down faster (from intuition and also from the equations), so it would have a less symmetric trajectory.

.

Footnote:
Who's this kindly old gentleman? (from wikimedia commons, Scanned from German "Meyer's Encyclopedia", 1906.) Turns out for modern bullets there's a special drag coefficient called the The Ballistic Coefficient . It is a number based on the drag coefficient, only normalized to an "ideal bullet" that this benign duffer produced in 1881, with ballistic coefficient of 1 by definition. (So now you know: if he's dictating standards, he's not just somebody's beloved grandpa). Unlike the drag coefficient, a smaller ballistic coefficient means more drag. Regular bullets in use today have ballistic coefs varying from about .01 (.177 pellets) to about 1, but you can make a bullet with BC 1.1 out of special materials with special lathes...better than the first "ideal bullet".

Cannon at Crécy: Balls in the Air! Writeup 60%, Code 40%

.

Battle of Crécy between the English and French in the Hundred Years' War. From a 15th-century illuminated manuscript of Jean Froissart's Chronicles (BNF, FR 2643, fol. 165v). (Commons.Wikimedia.org)

  1. Start with the base specifications above and experiment with changing the elevation and the mass (drag coefficient too, if you want). You should be able to predict from the equations what these changes will produce, and hence see nothing crazy is happening. Make and comment on a few plots.
  2. Given the set of parameters, what elevation maximizes range for drag coefficient of .5? (we know it's 45 degrees in a vacuum, but...). You can find the range as the value of x where y hits zero: with linear interpolation, you want to know where the line between the two consecutive trajectory points (xp,yp) and (xq,yq) crosses y=0, where yp is the last positive value of y and yq is the first negative value of y.
  3. The wallahs in PSYOPS think if would be bad for the Froggies' morale if you lobbed disgusting biological matter (like dead cats, rancid maggoty horsemeat, POW pieces and parts, fresh organic rutabagas...) into their camp. About how far can your cannon with muzzle velocity of 150m/s shoot a 15 cm diameter meatball? Why don't we say the drag coefficient c=.7 for this material, and you need to estimate its mass and justify your number.
  4. That cocky West-Bog-graduate lieutenant thinks variation in air density with altitude can make a difference to your range calculations. What is that variation (reference please, he's picky...thank goodness for Google)? What difference do you get with the base specifications? What if your muzzle velocity was 300 m/s? (Hint: you'll need to extend the derivative equations for this part).
  5. Enter z: You're shooting south to north and a freakish Sirocco wind is rising, blowing from due east to west. Let's call the east-to-west direction the postive z direction. At maximum range, how far will an 80mph (remember SI (MKS) units, now!) Sirocco make your shot deviate, and in which direction? Make a plot of the trajectory seen from above. (Hint, you will need to set the state variables for z and z', and compute z' and z'' in the derivative routine.)

    We're now in 3-D, and we are measuring the forces on the cannonball in a medium that is moving at windvel, the wind velocity, here 80mph in the z direction. We thus relativize the z velocity to the moving atmosphere, and compute the squared velocity:
    s2 = x'2 + y'2 + (windvel- z') 2.
    The formula for z'' then looks like that for x'', only with a γ to project the force onto the (relative) z' direction, instead of an α (see the Physics section above for α, etc. defs.).

    For debugging, Should try it with no x or y velocity and a very light cannonball (like a balloon).

    FWIW, I get a pretty significant z deflection for the base specifications and 45 degree shooting elevation.

    Optional: Code Organization, Approaches, Opinions, Hints:

    function dy = cannon_der(t,y)
    script shoot.m
    function x0 = zero_cross(x,y)
    function x = range(angle, diameter, mass)
    function xmax = max_range(low_angle, high_angle)
    function dy = cannon_air_der(t,y)
    

    function dy = cannon_der(t,y) This is the always-needed derivative equations whose handle you pass to your ODE solver. It specifies the derivatives for the state vector that consists, in some order that makes sense to you, of the quantities (x,y,z,x',y',z'). The output is a COLUMN 6-vector specifying their derivatives (x', y', z', x'', y'', z''). The z (East-West) direction is used in the last part of the question. Note that since there is no explicit time dependence in this problem, the t value will be ignored. It is needed only to fit the form used by ode23().

    script shoot.m This is simply the call to ode23() and various plot commands to help get acquainted with the problem and for parts 1 and 5 of the assignment.

    function x = range(der_func_handle, angle, diameter, mass) This is the basic cannon-shooting simulator that we can use in all parts of the assignment. It uses your ODE solver (ode23() is fine) to compute the range of the cannonball (where its height (y) goes to 0). Thus this function needs the handle to cannon_der, it needs to set initial conditions (using the angle) and a time span (experiment until you get one that gets y (height) into negative values for the range of angles you want to consider) and it calls zero-cross. Its parameters are the most commonly-changed ones (and you have to change them in the second part of the assignment).

    IMPORTANT! cannon_der() needs to use these parameters (angle, diameter, mass) but we can't pass them into cannon_der() since its parameters have to look like (t,y). So declare these variables global in both cannon_der() and range().

    function x0 = zero_cross(x,y) processes the output of ode23() to estimate the value of x for which the ball hits the ground (i.e. for which y(x) = 0). ode23's output consists of the varying state vector Y (each row is a [x y z x' y' z'] vector) and the associated x value vector, where Y = f(x). x0 is the computed estimate of the value of x for which y=0. One coarse estimate is the first x for which y changes sign from positive to negative.

    Better is a linear interpolation: Let's say y starts out positive, as in this assignment. Given (x1,y1), with y1 the last positive y, and its neighboring point (x2,y2), with y2 the first negative y, you can draw a little picture with the two points, a line between them, and the line y=0, and notice you want to solve for x0, where a point on the line is (x0,0). A similar triangles argument and a little high school algebra gives you x0 in terms of x1, y1, x2, and y2.

    function xmax = max_range(low_angle, high_angle) This answers the second part of the assignment. It calls range() for various angles (maybe let them vary freom low to high by 2 degrees?), gets the ranges and returns the maximum one. I'd use a single for-loop that initializes a variable max_range to the first range and best_angle to the first angle, and subsequently updates the max_range and best_angle to the current range and angle if the current range is longer than the previous max.

    function dy = cannon_air_der(t,y) For the fourth part, we need a new derivative function that is just like cannon_der (globals and all) but instead of being a constant, the air-density is a function of ρ, the sea-level density, and y (the height). That change must be made in the three derivative functions for x', y', and z'. Come up with a close-enough function (e.g. google 'air density table'), justify it with a reference, and carry on!

    Baseball meets Forensic Engineering in Monte Carlo: Balls to the Wall! Writeup 60%, Code 40%

    .
    Courtesy Lancaster Jethawks at http://www.jethawks.com/
  6. You doubtless saw this WSJ Article from 23 June 2010 (you may need your NetID and password to get into this database if you don't save all your WSJ's like I do)...so here's a local copy. It's about a ballpark often featuring a tailwind that increases the number of home runs. Some statistics are provided.

    The NYTimes has hired us to quantify (and thereby they hope to invalidate) their right-wing rival's story. Here's what we (well, you) will do. (Note: NYTimes needs references for all numbers, functions, facts, etc.). About 10 minutes googling should be all you need.

    The strategy is to simulate the trajectories of hits, see if each hint is a homer, and to see if we can replicate the claimed statistical improvement in the number of homers given different velocities of following winds.

    Hint: don't even think about writing your final program and starting to play with it. Small routines, thoroughly understood and debugged: that's the thing. For instance, you could test your monte carlo routine by only varying one parameter (elevation or velocity), using a simple uniform distribution, and making sure the data make sense). Then the other parameter. Only then do both. Continue this conservative approach.


Last Change: 11/30/2012