Printed and electronic copies of Modeling and Simulation in Python are available from No Starch Press and Bookshop.org and Amazon.
Modeling#
This chapter introduces the modeling framework we will use throughout the book, and works through our first example, using a simple model of physics to evaluate the claim that a penny falling from the height of the Empire State Building could kill someone if it hit them on the head. Also, you’ll see how to do computation in Python with units like meters and seconds.
This chapter is available as a Jupyter notebook where you can read the text, run the code, and work on the exercises. Click here to access the notebooks: https://allendowney.github.io/ModSimPy/.
The Modeling Framework#
This book is about modeling and simulating physical systems. The following diagram shows what I mean by modeling:
Starting in the lower left, the system is something in the real world we are interested in. To model the system, we have to decide which elements of the real world to include and which we can leave out. This process is called abstraction.
The result of abstraction is a model, which is a description of the system that includes only the features we think are essential. A model can be represented in the form of diagrams and equations, which can be used for mathematical analysis. It can also be implemented in the form of a computer program, which can run simulations.
The result of analysis and simulation might be a prediction about what the system will do, an explanation of why it behaves the way it does, or a design intended to achieve a purpose.
We can validate predictions and test designs by taking measurements from the real world and comparing the data we get with the results from analysis and simulation.
For any physical system, there are many possible models, each one including and excluding different features, or including different levels of detail. The goal of the modeling process is to find the model best suited to its purpose (prediction, explanation, or design).
Sometimes the best model is the most detailed. If we include more features, the model is more realistic, and we expect its predictions to be more accurate. But often a simpler model is better. If we include only the essential features and leave out the rest, we get models that are easier to work with, and the explanations they provide can be clearer and more compelling.
As an example, suppose someone asks you why the orbit of the Earth is elliptical. If you model the Earth and Sun as point masses (ignoring their actual size), compute the gravitational force between them using Newton’s law of universal gravitation, and compute the resulting orbit using Newton’s laws of motion, you can show that the result is an ellipse. Of course, the actual orbit of Earth is not a perfect ellipse, because of the gravitational forces of the Moon, Jupiter, and other objects in the solar system, and because Newton’s laws of motion are only approximately true (they don’t take into account relativistic effects). But adding these features to the model would not improve the explanation; more detail would only be a distraction from the fundamental cause. However, if the goal is to predict the position of the Earth with great precision, including more details might be necessary.
Choosing the best model depends on what the model is for. It is usually a good idea to start with a simple model, even if it is likely to be too simple, and test whether it is good enough for its purpose. Then you can add features gradually, starting with the ones you expect to be most essential. This process is called iterative modeling.
Comparing results of successive models provides a form of internal validation, so you can catch conceptual, mathematical, and software errors. And by adding and removing features, you can tell which ones have the biggest effect on the results, and which can be ignored.
Comparing results to data from the real world provides external validation, which is generally the strongest test.
The modeling framework is pretty abstract; the following example might make it clearer.
The Falling Penny Myth#
You might have heard that a penny dropped from the top of the Empire State Building would be going so fast when it hit the pavement that it would be embedded in the concrete; or if it hit a person, it would break their skull.
We can test this myth by making and analyzing two models. For the first model, we’ll assume that the effect of air resistance is small. In that case, the primary force acting on the penny is gravity, which causes the penny to accelerate downward.
If the initial velocity is 0 and the acceleration, \(a\), is constant, the velocity after \(t\) seconds is
and the distance the penny has dropped is
To find the time until the penny reaches the sidewalk, we can solve for \(t\):
Plugging in the acceleration of gravity, \(a = 9.8\) m/s\(^2\), and the height of the Empire State Building, \(x = 381\) m, we get \(t = 8.8\) s.
Then computing \(v = a t\) we get a velocity on impact of \(86\) m/s, which is about 190 miles per hour. That sounds like it could hurt.
Of course, these results are not exact because the model is based on simplifications. For example, we assume that gravity is constant. In fact, the force of gravity is different on different parts of the globe, and it gets weaker as you move away from the surface. But these differences are small, so ignoring them is probably a good choice for this problem.
On the other hand, ignoring air resistance is not a good choice, because in this scenario its effect is substantial. Once the penny gets to about 29 m/s, the upward force of air resistance equals the downward force of gravity, so the penny stops accelerating. This is the terminal velocity of the penny in air.
And that suggests a second model, where the penny accelerates until it reaches terminal velocity; after that, acceleration is 0 and velocity is constant. In this model, the penny hits the sidewalk at about 29 m/s. That’s much less than 86 m/s, which is what the first model predicts. Getting hit with a penny at that speed might hurt, but it would be unlikely to cause real harm. And it would not damage concrete.
The statistician George Box famously said “All models are wrong, but some are useful.” He was talking about statistical models, but his wise words apply to all kinds of models. Our first model, which ignores air resistance, is very wrong, and probably not useful. The second model, which takes air resistance into account, is still wrong, but it’s better, and it’s good enough to refute the myth.
The television show Mythbusters has tested the myth of the falling penny more carefully; you can view the results at https://modsimpy.com/myth/. Their work is based on a mathematical model of motion, measurements to determine the force of air resistance on a penny, and a physical model of a human head.
Computation In Python#
Let me show you how I computed the results from the previous section using Python. First we’ll create a variable to represent acceleration due to gravity in meters per second squared (m/s\(^2\)).
a = 9.8
A variable is a name that corresponds to a value. In this example, the name is a and the value is the number 9.8.
Suppose we let the penny drop for \(3.4\) seconds (s). I’ll create a variable to represent this time:
t = 3.4
Now we can compute the velocity of the penny after t seconds.
v = a * t
Python uses the symbol * for multiplication. The other arithmetic operators are + for addition, - for subtraction, / for division, and ** for exponentiation.
After you assign a value to a variable, you can display the value like this:
v
33.32
After \(3.4\) s, the velocity of the penny is about \(33\) m/s (ignoring air resistance). Now let’s see how far it would travel during that time:
x = a * t**2 / 2
x
56.644
It would travel about \(56\) m. Now, going in the other direction, let’s compute the time it takes to fall 381 m, the height of the Empire State Building.
h = 381
For this computation, we need the square root function, sqrt, which is provided by a library called NumPy.
Before we can use it, we have to import it like this:
from numpy import sqrt
Now we can use it like this:
t = sqrt(2 * h / a)
t
8.817885349720552
With no air resistance, it would take about \(8.8\) s for the penny to reach the sidewalk.
v = a * t
v
86.41527642726142
And its velocity on impact would be about \(86\) m/s.
False Precision#
Python displays results with about 16 digits, which gives the impression that they are very precise, but don’t be fooled. The numbers we get out are only as good as the numbers we put in.
For example, the “roof height” of the Empire State Building is \(380\) m (according to Wikipedia: https://en.wikipedia.org/wiki/Empire_State_Building). I chose \(h=381\) m for this example on the assumption that the observation deck is on the roof and you drop the penny from a 1 meter railing. But that’s probably not right, so we should treat this value as an approximation where only the first two digits are likely to be right.
If the precision of the inputs is two digits, the precision of the outputs is two digits, at most.
That’s why, if the output is 86.41527642726142, I report it as “about 86”.
Computation With Units#
The computations we just did use numbers without units.
When we set h=381, we left out the meters, and when we set a=9.8, we left out the meters per second squared.
And, when we got the result v=86, we added back the meters per second.
Leaving units out of computation is a common practice, but it tends to cause errors, including the very expensive failure of the Mars Climate Orbiter in 1999 (see https://en.wikipedia.org/wiki/Mars_Climate_Orbiter). When possible, it is better to include units in the computation.
To represent units, we’ll use a library called Pint (https://pint.readthedocs.io/).
To use it, we have to import a function called UnitRegistry and call it like this:
from pint import UnitRegistry
units = UnitRegistry()
The result is an object that contains variables representing pretty much every unit you’ve heard of. For example:
units.league
units.fortnight
But leagues and fortnights are not part of the International System of Units
(see https://en.wikipedia.org/wiki/International_System_of_Units); in the jargon, they are not “SI units”.
Instead, we will use meter and second.
meter = units.meter
meter
second = units.second
second
We can use meter and second to create a variable named a and give it the value of acceleration due to gravity.
a = 9.8 * meter / second**2
a
The result is a quantity with two parts, called magnitude and units, which we can access like this:
a.magnitude
9.8
a.units
We can also create a quantity that represents \(3.4\) s.
t = 3.4 * second
t
And use it to compute the distance a penny would fall after t seconds with constant acceleration a.
a * t**2 / 2
Notice that the units of the result are correct. If we create a quantity to represent the height of the Empire State Building:
h = 381 * meter
We can use it to compute the time the penny would take to reach the sidewalk.
t = sqrt(2 * h / a)
t
And the velocity of the penny on impact:
v = a * t
v
As in the previous section, the result is about \(86\), but now it has the correct units, m/s.
With Pint quantities, we can convert from one set of units to another like this:
mile = units.mile
hour = units.hour
v.to(mile/hour)
If you are more familiar with miles per hour, this result might be easier to interpret. And it might give you a sense that this model is not realistic.
Summary#
This chapter introduces a modeling framework that consists of three steps:
Abstraction is the process of defining a model by deciding which elements of the real world to include and which can be left out.
Analysis and simulation are ways to use a model to generate predictions, explain why things behave as they do, and design things that behave as we want.
Validation is how we test whether the model is right, often by comparing predictions with measurements from the real world.
As a first example, we modeled a penny dropped from the Empire State building, including gravity but ignoring air resistance. In the exercises, you’ll have a chance to try a better model, including air resistance.
This chapter also presents Pint, a library for doing computation with units, which is convenient for converting between different units and helpful for avoiding catastrophic errors.
Exercises#
Before you go on, you might want to work on the following exercises.
Exercise 1#
In mathematical notation, we can write an equation like \(v = a t\) and it’s understood that we are multiplying \(a\) and \(t\). But that doesn’t work in Python. If you put two variables side-by-side, like this:
v = a t
you’ll get a syntax error, which means that something is wrong with the structure of the program. Try it out so you see what the error message looks like.
Exercise 2#
In this chapter we used the sqrt function from the NumPy library. NumPy also provides a variable named pi that contains an approximation of the mathematical constant \(\pi\).
We can import it like this:
from numpy import pi
pi
3.141592653589793
NumPy provides other functions we’ll use, including log, exp, sin, and cos.
Import sin and cos from NumPy and compute
Note: A mathematical identity tells us that the answer should be \(1\).
Show code cell content
# Solution
from numpy import sin, cos
sin(pi/4)**2 + cos(pi/4)**2
1.0
Exercise 3#
Suppose you bring a 10 foot pole to the top of the Empire State Building and use it to drop the penny from h plus 10 feet.
Define a variable named foot that contains the unit foot provided by the UnitRegistry we called units. Define a variable named pole_height and give it the value 10 feet.
What happens if you add h, which is in units of meters, to pole_height, which is in units of feet? What happens if you write the addition the other way around?
h = 381 * meter
Show code cell content
# Solution
foot = units.foot
pole_height = 10 * foot
h + pole_height
Show code cell content
# Solution
pole_height + h
Exercise 4#
Why would it be nonsensical to add a and t? What happens if you try?
a = 9.8 * meter / second**2
t = 3.4 * second
In this example, you should get a DimensionalityError, which indicates that you have violated a rule of dimensional analysis: you cannot add quantities with different dimensions.
The error messages you get from Python are big and scary, but if you read them carefully, they contain a lot of useful information.
The last line usually tells you what type of error happened, and sometimes additional information, so you might want to start from the bottom and read up.
The previous lines are a traceback of what was happening when the error occurred. The first section of the traceback shows the code you wrote. The following sections are often from Python libraries.
Exercise 5#
Suppose instead of dropping the penny, you throw it downward at its terminal velocity, \(29\) m/s. How long would it take to fall \(381\) m?
Show code cell content
# Solution
h = 381 * meter
v = 29 * meter / second
t = h / v
t
Exercise 6:#
So far we have considered two models of a falling penny:
If we ignore air resistance, the penny falls with constant acceleration, and we can compute the time to reach the sidewalk and the velocity of the penny when it gets there.
If we take air resistance into account, and drop the penny at its terminal velocity, it falls with constant velocity.
Now let’s consider a third model that includes elements of the first two: let’s assume that the acceleration of the penny is a until the penny reaches \(29\) m/s, and then \(0\) m/s\(^2\) afterwards. What is the total time for the penny to fall \(381\) m?
You can break this question into three parts:
How long would the penny take to reach \(29\) m/s with constant acceleration
a.How far would it fall during that time?
How long would it take to fall the remaining distance with constant velocity \(29\) m/s?
Suggestion: Assign each intermediate result to a variable with a meaningful name. And assign units to all quantities!
a = 9.8 * meter / second**2
h = 381 * meter
Show code cell content
# Solution
v_terminal = 29 * meter / second
Show code cell content
# Solution
t1 = v_terminal / a
print('Time to reach terminal velocity', t1)
Time to reach terminal velocity 2.9591836734693877 second
Show code cell content
# Solution
h1 = a * t1**2 / 2
print('Height fallen in t1', h1)
Height fallen in t1 42.90816326530612 meter
Show code cell content
# Solution
t2 = (h - h1) / v_terminal
print('Time to fall remaining distance', t2)
Time to fall remaining distance 11.658339197748065 second
Show code cell content
# Solution
t_total = t1 + t2
print('Total falling time', t_total)
Total falling time 14.617522871217453 second
Exercise 7#
When I was in high school, the pitcher on the baseball team claimed that when he threw a fastball he was throwing the ball down; that is, the ball left his hand at a downward angle. I was skeptical; watching from the side, I thought the ball left his hand at an upward angle.
Can you think of a simple model you could use to settle the argument? What factors would you include and what could you ignore? What quantities would you have to look up or estimate?
I suggest you convert all quantities to SI units like meters and seconds (see https://en.wikipedia.org/wiki/International_System_of_Units).
Show code cell content
# Solution
# I suggest the following model:
# 1. Let's ignore the motion of the ball toward home plate and
# think about how far the ball would drop while it's in flight.
# 2. Let's ignore air resistance. Since we are only thinking
# about the relatively slow motion in the vertical direction,
# this is probably a good assumption.
# 3. Let's also ignore the effect of spin. This is probably a
# less good assumption.
# The distance from the pitcher's mound to home plate is about 60
# feet, but the point where the ball is released is a bit closer.
# An average pitcher in high school might be able to throw a ball
# at 80 mph.
Show code cell content
# Solution
v = (80 * mile / hour).to(meter/second)
x = (60 * foot).to(meter)
Show code cell content
# Solution
t = x / v
t
Show code cell content
# Solution
a = 9.8 * meter / second**2
h = a * t**2 / 2
h
Show code cell content
# Solution
# In the time it takes the ball to reach home plate, it drops
# about 1.3 m. If the release point is at 2 m, which is plausible,
# it would cross the plate at 0.7 m, which is in the strike zone.
# So I could be wrong -- it is plausible that the ball leaves
# the pitcher's hand at a downward angle, at least for some pitches.
Exercise 8#
Suppose I run a 10K race in 44:52 (44 minutes and 52 seconds). What is my average pace in minutes per mile?
mile = units.mile
kilometer = units.kilometer
minute = units.minute
Show code cell content
# Solution
t = 52 * second + 44 * minute
t
Show code cell content
# Solution
v = 10 * kilometer / t
v
Show code cell content
# Solution
pace_per_km = (1 / v)
pace_per_km
Show code cell content
# Solution
pace_per_mile = pace_per_km.to(minute / mile)
pace_per_mile
Show code cell content
# Solution
# To convert to minutes and seconds, we can use np.round.
# But we haven't covered that yet.
from numpy import round
round(pace_per_mile)
Show code cell content
# Solution
remainder = pace_per_mile - round(pace_per_mile)
remainder.to(second/mile)