# PHYC 2050 - Computer Simulation in Science

In Computer Simulation in Science we learn the python programming language and use its computational capabilities to explore problems in science. As a key example, we develop a model for the projectile motion of a golf ball. Optimization techniques are used to find the best launch angles and spin rates for different conditions and drive speeds. Monte Carlo methods are introduced as a way of dealing with real variability in a system. The class ends with a projectile motion programming game where students pit their programming and problem-solving skills against each other.

## News and Upcoming Events

**Wednesday, 21 April 2010, 1-4 PM in Dunn 208:**Final exam. Note that this is a different time than was scheduled by the Registrar's Office.- 9 April 2010: Assignment 4 is due at least one hour before class
- 24 March 2010: Ch. 6 of the lecture notes released for download
- 19 March 2010: Second test, during class time in Room 208
- 17 March 2010: Ch. 5 of the lecture notes released for download
- 15 March 2010: Assignment 3 is due at the beginning of class
- 19 February 2010: Midterm exam, during class time in Room 208
- 17 February 2010: Ch. 3 and 4 of the lecture notes released for download
- 10 February 2010: Assignment 2 is due at the beginning of class
- 2 February 2010: Ch. 2 of the lecture notes released for download
- 25 January 2010: First test, during class time in Room 208
- 20 January 2010: Ch. 1 of the lecture notes released for download
- 20 January 2010: Assignment 1 is due at the beginning of class
- 7 January 2010: Jason Hopper will provide a second tutorial introduction to using linux in Room 208 at 3:30 PM
- 5 January 2010: Jason Hopper will provide a tutorial introduction to using linux in Room 208 at 3 PM

## Reference Documents

There is no required text for this class. All of the documentation that we will need is available online.

- Introduction to Unix - for those new to linux/unix
- Unix Command Summary - a list of the most used shell commands
- Python - the interactive, interpreted programming language
- Documentation index - we are using python 2.6
- The Python Tutorial - a supplement to what you learn in class
- The Python Standard Library - more modules than you can shake a stick at

- numpy - essential computational facilities for python
- Tutorial
- numpy functions by category - with examples

- scipy - advanced scientific computing for python
- matplotlib/pylab - matlab-style plotting for python
- A Guided Tour of Emacs - my favourite text editor (oblig. xkcd)

## Downloads

- Class syllabus
- Lecture notes: Only released as I am able, so be sure to take your own notes in class.
- Assignments
- Assignment 1
- Assignment 2
- Assignment 3 (updated Mar. 11)
- Assignment 4

- Programs
The following appeared in class, in the lecture notes, or both

- Hello, World!
- Fitting a straight line to data
- linefit1.py -- points on a line with noise (plot)
- linefit2.py -- fit to the points (plot)
- linefit3.py -- includes a function for fitting
- linefit4.py -- calls fit function in tools module

- Projectile motion: Cannon shells
These programs were only used in class. There are analogous programs for golf ball trajectories in the lecture notes.

- vacuum_1d.py -- vertical projectile motion in a vacuum
- projectile.py -- 2D projectile motion with drag

- Projectile motion: Golf
- golf_vacuum_1d.py -- vertical projectile motion in a vacuum (plot)
- golf_vacuum_2d.py -- 2D projectile motion in a vacuum (plot)
- golf_drag.py -- golf ball trajectory with the air resistance (plot)
- golf_magnus.py -- golf ball trajectory including the magnus force (plot)
- golf.py -- the golf module (updated March 26)
- cgolf.py -- the golf module with embedded C code for speed (updated March 26)
- golf_trajectories.py -- test program used to change parameters and see the results (plot)
- golf_angles.py -- plots the travel distance against the launch angle (plot)
- golf_parameter_space.py -- total drive distance versus angle and spin rate (plot)
- golf_angles_findmax.py -- optimizes the launch angle for maximum distance given the launch speed and rotation rate
- golf_angle_optimization.py -- optimizes the launch angle against speed and spin rate (plot)
- golf_angle_optimization_wind.py -- optimizes the launch angle against speed and spin rate for different wind speeds (plot)
- golf_monte_carlo.py -- statistical golf simulation with random variability (plot)

- Random numbers and PDFs
- gauss.py -- some sample Gauss distributions (plot)
- gauss_accumulated.py -- corresponding accumulated distributions (plot)
- pseudorandom.py -- output from the Linear Congruential Generator (plot)
- gauss_random.py -- Gauss-distributed noise (plot)
- gamma_random.py -- Gamma-distributed noise (plot)

- Combat Golf
The most recent update was April 8 @ 4:45 PM. Please make sure you are using the most up-to-date version of each file.

- README -- important notes
- LICENSE -- distribution licences
- combat_golf.py -- main program
- game_engine.py -- the game engine and plotter
- golf.py -- physics module
- cgolf.py -- C extensions to golf.py
- random_golfer.py -- randomly firing golfer control module
- target.py -- a target golfer that doesn't fire back
- test.py -- unit tests
- ks.py, jl.py, ap.py, bt.py, cr.py, gh.py lb.py -- student-written player modules

## Interesting Links

- The Hello World Collection - in over 400 different programming languages
- 99 Bottles of Beer - programmed in 1314 variations (c, fortran, java, perl, python)

## Class Log

- Introduction
- Basic mathematics
- Functions and modules
- Variables; Arrays
- Function definition and conditional statements; Plotting
- Hello World; Fitting a line to data
- Review
- Taylor polynomials and Euler's method
- Test 1
- Euler's method (continued); Array indexing; while loops
- For loops; Array slicing; Difference formula implementations
- Numbers with different bases; Representation of integers
- Representation of floating-point numbers
- Machine epsilon
- 1 and 2D Projectile motion
- Projectile motion with drag
- The Magnus force
- Snow day -- class cancelled
- Midterm
- STUDY BREAK
- The Magnus force revisited
- Parameterizations for lift and drag and the flight of a golf ball
- Spin decay, carry and roll
- Implementation and Code Optimization
- Mixed language programming; Exploring the parameter space
- More exploring; The Brent method; Optimizing the launch angle
- Optimizing the angle against launch speed and spin rate; Optimizing with a tail wind
- Intro to Monte Carlo methods
- Test 2
- Probability density functions (PDFs) and the accumulated probability; The Linear Congruential Generator
- Noise distributed according to a PDF
- Golf Monte Carlo; Introduction to Combat Golf
- Object-oriented programming
- More object-oriented programming
- Dictionaries; Error handling; Sorting
- Finding and extracting things from lists; The in, continue and break keywords
- Combat Golf Day

## FAQ (Frequently Asked Questions)

(... although the frequency may only be once.)

- Can you provide technical support for my home computer?
- Is it essential that I use linux? I can run python in Windows.
- How can I set up my home computer to run python with the numpy, scipy and matplotlib extensions?

Sorry, no. However, I am happy to provide advice where I can.
*Disclaimer: you are responsible for what you do to your own
computer, not me!!*

[Answer: +]

[Answer: +]

## Exercises

The following problems are left as exercises (no mark value). Solutions will be posted as they are submitted to me by students in the class.- 4 January 2010
- 6 January 2010
- 8 January 2010
- 11 January 2010
- 13 January 2010
- 15 January 2010
- 18 January 2010
- 20 January 2010
- 22 January 2010
- 25 January 2010
- 27 January 2010
- 29 January 2010
- 1 February 2010
- 3 February 2010
- 5 February 2010
- 8 February 2010
- 10 February 2010
- 12 February 2010
- 15 February 2010
- 17 February 2010
- 19 February 2010
- 1 March 2010
- 3 March 2010
- 5 March 2010
- A ball launched with a 5 iron at 189.5 ft/s, 16.6 deg and 91.1 rev/s which carried 187.6 yards.
- A ball launched with a sand wedge at 130.1 ft/s, 25.7 deg and 178.8 rev/s which carried 110.7 yards.
- 8 March 2010
- 10 March 2010
- 12 March 2010
- 15 March 2010
- 17 March 2010
- 19 March 2010
- 22 March 2010
- 24 March 2010

Q: Determine the types returned by operations using the following binary arithmetic operators: +, /, %, **, <, and ==. Use both ints, floats, and combinations of the two.

[Answer: +]

Q: Explain the following result from the python interpreter, which would be 1 if typed into a standard calculator.

>>> (7+2)/7 * 7/(7+2.) 0.77777777777777779

You can ignore the problem with the last digit of the result, which we will deal with later.

[Answer: +]

Q: How can you round an arbitrary *float* to the nearest *int*
in python? (Provide examples)

[Answer: +]

Q: Numpy has a variety of functions for creating arrays. How can we create an array of integer fives? How can we convert this array to an array of floats instead?

[Answer: +]

Q: Functions that call themselves are called recursive. Write a recursive function f(n) that evaluates the arithmetic series 1+2+3+...+n. To do this recursively, your function f(n) should return n + f(n-1). Be careful: you will need to provide a way for the recursion to end!

[Answer: +]

Q: It is possible to pass a function as an argument in python. For example:

>>> def foo(func,x): ... return func(x) ... >>> def square(x): ... return x**2 ... >>> foo(square,3) 9

Functions in python are said to be "first class" because they can be treated just like objects.

The built-in map() function can be used to apply a function of one variable to all of the elements in an array individually. Use the map() function to apply the math.sin() function to the range of integers between 0 and 10.

[Answer: +]

Q: Lists, tuples, strings, and numpy arrays are three examples of sequence
types in python. They are, however, distinct types, and so can be
expected to behave differently in different circumstances.

a) What is the primary difference between a list and an array?

b) Show that the multiplication of list by an integer produces
different results than for an array multiplied by an integer.
Describe the behaviour in each case.

c) Show that adding two lists together produces different results than
for adding two arrays. Describe the behaviour in each case.

[No answer yet]

Q: A function may explicitly return a single value, a comma-separate sequence
of values, or nothing at all. What *actually* gets returned may
surprise you.

a) What type is returned when a comma-separated sequence is returned? How are objects of this type different than the list and array sequence types?

b) What type is returned when nothing is explicitly returned in the function body? Are there any values associated with this type? Note: If you want to do nothing in the body of a function, use "pass", which is one of the reserved words in python. You can use "pass" any time you want to leave a block of code empty.

[Answer: +]

Q: For the exercise on January 13 you wrote a function that evaluated the arithmetic series (1+2+3+...+n) using recursion.

a) Write the same function using only functions in numpy (no recursion).

b) Write the same function using a loop (no recursion, no numpy).

[Answer: +]

Q: The python *dictionary* type allows the mapping of arbitrary keys
to values. Look up the dictionary type in the
Python Tutorial. Create
a dictionary that maps the names of each month to its corresponding numerical
value. Show how values in the dictionary can be accessed.

[Answer: +]

Q: Derive the central difference formula for the *second* derivative:
$$ f_i'' = \frac{f_{i+1} - 2f_i + f_{i-1}}{h^2} + O(h^2) $$
If you want your answer to look pretty, format it
LaTeX-style (see this
short math tutorial for LaTeX). For example, the source for the
above equation is

f_i'' = \frac{f_{i+1} - 2f_i + f_{i-1}}{h^2} + O(h^2)

[No answer yet]

How can you slice an array [0,1,2,3,4,5,6,7,8,9] such that the
result is [8,6,4,2,0]? Recall that the slice notation for array *A* is
*A[start:stop:step]*, where *start*, *stop* and *step* are integers.

[Answer: +]

Q: Can 10.1 be represented exactly as a "double precision" floating-point number? Explain using the formalism provided in class.

[Answer: +]

Q: List comprehensions in python can achieve the same ends as the map() function, but are more powerful and generally considered better form. Rewrite the factorial function to replace map() with a list comprehension when the factorial argument is an array.

[Answer: +]

Q: Single precision floats represent numbers in 32 bits.

a) How are the bits partitioned into the different parts of a float?

b) What is the largest and smallest positive single-precision float?

c) What is machine epsilon for single precision floats?

[Answer: +]

Q: What is the largest argument to the sine function we can choose such that the maximum error due to machine precision is on the order of 0.01 ? Note: This problem requires a theoretical answer, not a computation.

[No answer yet]

Q: Modify the 1D projectile motion in a vacuum program so that it determines and plots the trajectory in 2D for a launch angle of 45 degrees.

[Answer: +]

Q: Many physical problems (transmission of light in a medium, radioactive decay, etc) obey the following differential equation: $$ \frac{df}{dt} = -\alpha f $$ for decay constant alpha.

a) Solve the equation analytically.

b) Choose $$ \alpha = 1\ s $$ and solve the equation numerically. How do the computed and analytic solutions compare?

[No answer yet]

Q: The differential equation describing the variation of pressure with altitude in the atmosphere is $$ \frac{d\ln p}{dz} = -\frac{g}{RT} $$ where the temperatures (K) approximately decrease with altitude according to $$ T(z) = 290 - C z $$ for gravitational acceleration g, gas constant R = 287 J/K/km and lapse rate C = 6.5 K/km.

Determine the atmospheric pressure profile numerically (be careful with units!). Plot T and p in two different panels on the same Figure using pylab's subplot() functionality.

[No answer yet]

Q: Normally we provide a print statement with a comma-separated list of arguments.

from numpy import pi print 'Pi is approximately', pi

This form, however, does not provide very good control over the output.
We may, for example, want to limit the number of decimal places printed out.
Such control is possible through the use of *string formatting*.
Read up on "old-style" (very C-like) and
new style
(using a string's format method) string formatting
in the python documentation. Rewrite the example above in both styles so
that pi is printed with only two decimal places.

[Answer: +]

Q: Read up on how to write data to a file in the python documentation. Create a program that writes the first N values of the geometric progression $$ 1,\ \frac{1}{2},\ \frac{1}{4},\ \frac{1}{8},\ \ldots $$ to a file, arranged in a single column (N should be defined in the program).

[Answer: +]

Q: The optparse module provides a command-line option parser for use with your programs. Modify the program in the previous exercise so that it reads in the variable N from the command line. For example, executing

$ python geometric.py 10

should set N as 10 and write the geometric progression as before.

[Answer: +]

Q: Numerically solve the spin decay equation, $$ \dot\omega = -C_\omega\frac{\omega V}{r} $$ assuming the initial spin rate is 3000 rpm, the speed is a constant 70 m/s, r is the radius of a golf ball, and the decay constant is 2.0e-5. Plot the result.

[Answer: +]

Q: Using the new golf.py module (available in the downloads section), compare the model data against the following measured test data from the USGA "Second Report on the Study of Spin Generation":

[Answer: +]

Q: Write a program that plots the spin rate versus time for a golf ball launched at 70 m/s, 15 degrees, and 3000 rpm. Compare the result to that obtained from the exercise on March 3.

[No answer yet]

Q: Write a program that computes and plots the optimum launch angle versus spin rate for a golf ball with V0=45 m/s.

[No answer yet]

Q: Numerically compute the zeros and minimum of the equation $$ y = \frac{x^3}{10} - x^2 + 10 $$ for 0<x<10.

[Answer: +]

Q: Use the Brent Method to compute where the curves $$ y = \sin x $$ and $$ y = x/3 $$ cross.

[Answer: +]

Q: Create a contour plot for the equation $$ z = x^2 - y^2 $$ for -5<x<5 and -4<y<4

[Answer: +]

Q: scipy.optimize.fmin() locates the minimum of a function
of one *or more* variables by starting at a given point and moving
downhill. Use the function to locate the minimum of
$$ z = x^2 + y^2 $$
Although the location of the minimum is obvious, you should choose a
starting point somewhere away from that to prove to yourself the
optimization works.

[No answer yet]

Q: Is the sum of two Gauss-distributed random arrays also Gauss-distributed? What about the product? Write a program that generates two Gauss-distributed random arrays with different means and standard deviations. Test each case by comparing histograms against a Gauss distribution using the measured mean and standard deviation. Hint: One case is Gauss-distributed, and one isn't.

[No answer yet]

Q: Write a program that creates a straight line modulated by Gaussian-distributed noise. Fit a straight-line to the data and compare the statistical parameters of the fit to the known inputs.

[No answer yet]