# 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.

## 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

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

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

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

## 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.

3. 6 January 2010
4. 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.

5. 8 January 2010
6. Q: How can you round an arbitrary float to the nearest int in python? (Provide examples)

7. 11 January 2010
8. 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?

9. 13 January 2010
10. 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!

11. 15 January 2010
12. 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.

13. 18 January 2010
14. 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.

15. 20 January 2010
16. 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.

17. 22 January 2010
18. 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).

19. 25 January 2010
20. 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.

21. 27 January 2010
22. 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)


23. 29 January 2010
24. 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.

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

27. 3 February 2010
28. 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.

29. 5 February 2010
30. 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?

31. 8 February 2010
32. 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.

33. 10 February 2010
34. 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.

35. 12 February 2010
36. 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?

37. 15 February 2010
38. 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.

39. 17 February 2010
40. 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.

41. 19 February 2010
42. 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).

43. 1 March 2010
44. 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.

45. 3 March 2010
46. 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.

47. 5 March 2010
48. 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.

49. 8 March 2010
50. 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.

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

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

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

57. 17 March 2010
58. Q: Create a contour plot for the equation $$z = x^2 - y^2$$ for -5<x<5 and -4<y<4

59. 19 March 2010
60. 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.