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.

Downloads

Interesting Links

Class Log

  1. Introduction
  2. Basic mathematics
  3. Functions and modules
  4. Variables; Arrays
  5. Function definition and conditional statements; Plotting
  6. Hello World; Fitting a line to data
  7. Review
  8. Taylor polynomials and Euler's method
  9. Test 1
  10. Euler's method (continued); Array indexing; while loops
  11. For loops; Array slicing; Difference formula implementations
  12. Numbers with different bases; Representation of integers
  13. Representation of floating-point numbers
  14. Machine epsilon
  15. 1 and 2D Projectile motion
  16. Projectile motion with drag
  17. The Magnus force
  18. Snow day -- class cancelled
  19. Midterm
  20. STUDY BREAK
  21. The Magnus force revisited
  22. Parameterizations for lift and drag and the flight of a golf ball
  23. Spin decay, carry and roll
  24. Implementation and Code Optimization
  25. Mixed language programming; Exploring the parameter space
  26. More exploring; The Brent method; Optimizing the launch angle
  27. Optimizing the angle against launch speed and spin rate; Optimizing with a tail wind
  28. Intro to Monte Carlo methods
  29. Test 2
  30. Probability density functions (PDFs) and the accumulated probability; The Linear Congruential Generator
  31. Noise distributed according to a PDF
  32. Golf Monte Carlo; Introduction to Combat Golf
  33. Object-oriented programming
  34. More object-oriented programming
  35. Dictionaries; Error handling; Sorting
  36. Finding and extracting things from lists; The in, continue and break keywords
  37. Combat Golf Day

FAQ (Frequently Asked Questions)

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

  1. Can you provide technical support for my home computer?

    2. 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!!

  2. Is it essential that I use linux? I can run python in Windows.

    3. [Answer: +]

  3. How can I set up my home computer to run python with the numpy, scipy and matplotlib extensions?

    4. [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.
  1. 4 January 2010

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

    [Answer: +]


  2. 6 January 2010

    3. 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: +]


  3. 8 January 2010

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

    [Answer: +]


  4. 11 January 2010

    5. 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: +]


  5. 13 January 2010

    6. 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: +]


  6. 15 January 2010

    7. 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: +]


  7. 18 January 2010

    8. 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]


  8. 20 January 2010

    9. 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: +]


  9. 22 January 2010

    10. 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: +]


  10. 25 January 2010

    11. 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: +]


  11. 27 January 2010

    12. 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]


  12. 29 January 2010

    13. 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: +]


  13. 1 February 2010

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

    [Answer: +]


  14. 3 February 2010

    15. 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: +]


  15. 5 February 2010

    16. 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: +]


  16. 8 February 2010

    17. 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]


  17. 10 February 2010

    18. 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: +]


  18. 12 February 2010

    19. 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]


  19. 15 February 2010

    20. 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]


  20. 17 February 2010

    21. 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: +]


  21. 19 February 2010

    22. 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: +]


  22. 1 March 2010

    23. 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: +]


  23. 3 March 2010

    24. 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: +]


  24. 5 March 2010

    25. 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":

    1. 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.
    2. 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.

    [Answer: +]


  25. 8 March 2010

    26. 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]


  26. 10 March 2010

    27. 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]


  27. 12 March 2010

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

    [Answer: +]


  28. 15 March 2010

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

    [Answer: +]


  29. 17 March 2010

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

    [Answer: +]


  30. 19 March 2010

    31. 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]


  31. 22 March 2010

    32. 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]


  32. 24 March 2010

    33. 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]