Partnership for Integration of Computation into Undergraduate Physics

These exercises are designed to help you get started in using computers to solve physics problems. They use a language called "Python", which is an extremely powerful programming language useful for just about anything. This is not a programming course, though; to really learn Python will take a different course. This is just a set of exercises that will help you get started in using Python. The exercises build on each other, and should mostly be done in order. ## Exercise 1: Installation and Test Python is an open-source language that is extremely powerful yet still easy to use. It's also free, which is generally a good thing. Python by itself is a fairly basic computer language, but there are extra "libraries" that extend its capabilities to allow interesting things like plotting equations and data and 3D objects. The libraries you will need for these exercises are Numpy, Scipy, and MatPlotLib, which are often known collectively as "pylab". If you are already using Python, you may have these packages already; or you may need to do a more customized installation. If not, the simplest way to install what you need is to install "Anaconda". Anaconda is a full Python environment with all the tools you will need (and more!) It is available for free download from [continuum](//www.continuum.io/downloads). Download the appropriate version for your computer system, and follow the instructions to install it. Once Anaconda is installed, launch it and choose the "spyder-app" Scientific Python Development Environment. (Feel free to experiment with the others also, "qtconsole" is also a good choice.) A new window should appear, with several sub-windows. On the left will be a text editor, where you can write Python programs. On the bottom right will be a console window, in which you can type Python commands directly. In the console window, type the following commands, one per line: (Don't type the '>>>' symbols, they are just the Python prompt.) ``` >>> from pylab import * >>> t = linspace(0, 10*pi, 1000) >>> z = sin(t)*exp(-t/10) >>> plot(t,z,'b-') >>> title("Cool! It works!") >>> xlabel("time") >>> ylabel("displacement") >>> show() ``` After entering the last of these lines, a graph like the one below will appear in the console window.  Here's an explanation of what each of those lines do, for future reference. ``` >>> from pylab import * ``` This command tells Python to load up a big batch of mathematical goodies, and should be the first command you enter just about every time you use Python for graphing or numeric work. ``` >>> t = linspace(0, 10*pi, 1000) ``` $t$ is now an array of 1000 values, linearly spaced from 0 to 10$\pi$. ``` >>> z = sin(t)*exp(-t/10) ``` $z$ is calculated with the formula given. Since $t$ is an array, $z$ is also an array; one value for each value of $t$. ``` >>> plot(t,z,'b-') ``` Make a plot with the values contained in $t$ on the horizontal axis and the values contained in $z$ on the vertical axis. Show the data with a blue line ($b-$). ``` >>> title("Cool! It works!") >>> xlabel("time") >>> ylabel("displacement") ``` Make a title, a label for the horizontal axis, and a label for the vertical axis. (These are optional plot features, if you skip them the plot just won't have these useful details.) ``` >>> show() ``` Show the resulting plot. ### Assignment Make a change to the plot, print it out, and turn it in with a note calling attention to your change. ## Exercise 2: Writing a program In the previous exercise, you typed commands in one at a time. If you had to do things more than once, this would quickly become annoying. The usual way to do Python is to put all the commands into a text file, and then tell Python "go do those commands". Such a text file is called a "program". Instead of typing those commands in again and again each time we want to plot a mathematical function, we could make a program to do the same thing. We do this by creating a new text file containing those commands, in their proper order. We then save the file with a meaningful name, such as "plot1.py". (Do not put spaces in your file name, and end the name with ".py") In Spyder, we can then press the F5 key and Spyder will run the program. You can make changes to the program and run it again, as needed. Try it! Copy the Python code from exercise 1 and paste it into a new text editor window in Spyder, save the program, and press F5. Try changing things: change $\sin()$ to $\cos()$, change the title, whatever... just play with it a bit and see what happens. Any time you want to check how things are going, save and press F5 again to run the program. ### Functions One of the very useful tricks in Python is the ability to write functions. The $\sin(t)$ in the above program, for example, is a function. In the case of $\sin()$, that function is one of the things defined when you import pylab so you don't have to tell Python how to calculate sine. You can write your own functions with the 'def' (define) command. This allows you to simplify your programs into manageable chunks. Here's an example of how it works, using a rather silly and useless function. ```python def silly(x): y = sqrt(x**2 + 9) return y ``` There are several key things to note on that example. First, the colon and the indentation are part of the code! The colon tells Python to expect an indented chunk of code. All of the indented lines are part of the definition. Next, the $x$ in the definition is a place-holder. You can call silly(42), or silly(q), or silly(_anything_), as long as _anything_ is something that can be squared ($x$\*\*2) in the second line of the definition. Finally, the "return" value is what goes back to the rest of the program when you call this function. Once a function is defined in a program, we can use it just like we would use any more standard mathematical function. For example, the following code asks the user for a number and then prints "silly" of that number. ```python n = input("Enter a number: ") print("Silly of that number is", silly(n)) ``` (Oh, and by the way, input() and print() are simple ways to put numbers into a program and get numbers out. They are also Python functions, and are built-in.) This code works with arrays of numbers also: ```python a = linspace(0, 10, 11) print(silly(a)) plot(a, silly(a), 'go') plot(a, silly(a) + 1, 'rx') show() ``` ### Assignment If you drop an anvil from a great height, it accelerates downward at about 9.8 m/s$^2$. Write a program containing two functions: one that returns $v(t)$ and another that returns $y(t)$ for the anvil. Your program should then plot both these functions on one graph, for time ranging from 0 to 5 seconds. Have the program label the axes appropriately and put your name in the title. You could do this just using individual commands, as with the first exercise, but to increase your chance of surviving the rest of the exercises do it as a program, with functions. ## Exercise 3: Solving a hard problem If you launch a projectile off the top of a cliff, with initial height $h$, the range $R$ at which the projectile strikes the ground will depend on the initial speed $v_o$, the launch angle $\theta$, the initial height $h$, and the gravitational field $g$. The relationship between $R$, $v_o$, $g$, and $\theta$ in this case is *not* simply the "range equation" $$R = \frac{v_o^2\sin(2\theta)}{g}$$ The reason is that equation is derived for the case in which the launch point and impact point are at the same elevation, which is not the case here. ### Assignment part A 1. Derive an equation for $R$ as a function of $v_o$, $\theta$, $g$, and $h$. 2. Check your answer by letting $h$ go to zero: the result should simplify to the same-elevation answer. ### And now for some computer work The $h$-dependent range equation derived in part (1) above makes it easy to figure out how far something will go for a given launch angle. But what if you want to know what launch angle to use to hit a target at a known range $R$? It's possible to invert the equation, but it's quite difficult algebraically. Instead, let's use Python. What we have is $R = f(\theta)$, assuming that $v_o$ and $h$ are constants. The usual way to solve things computationally is to rearrange this like so: $$F(\theta) \equiv f(\theta) - R$$ Having done that, we can plot $F(\theta)$, and graphically determine the value of $\theta$ for which $F(\theta) = 0$. ### Assignment part B Using what you've learned so far, write a function $F(\theta)$ that returns the value of $f(\theta) - R$. Plot $F(\theta)$ versus $\theta$, and using the zoom tool on the resulting graph determine the value of $\theta$ for which $F(\theta)=0$. Report your answers to at least two decimal places. Be sure to save your work; you'll need it for the next exercise! Use these values: $R = 2.5$ m, $h = 1.2$ m, $v_o = 4.8$ m/s. There may be more than one correct answer. ## Exercise 4: Root-finding In the previous exercise, the computer problem involved using a graph to find the "roots" of an inconveniently-complex algebraic function. A root of function $f(x)$ is an $x$ value for which $f(x)=0$. Manually zooming into a graph works, but there are more convenient methods of numerical root-finding. One of the better alternative methods is the brentq() function. Brentq() requires three parameters: the function for which to find a root, and two starting guesses of where the answer might be. The two guesses must be on either side of a solution. Brentq() then returns the value of the solution between the two guesses. Let's use the problem from that set as an example: ```python from pylab import * # get mathematical tools from scipy.optimize import brentq # get the brentq() root-finder R = 2.5 h = 1.2 v = 4.8 g = 9.8 def F(theta): output = v*cos(theta)/g * (v*sin(theta) + sqrt(v**2*sin(theta)**2 + 2*g*h)) - R return output answer1 = brentq(F, 0.0, 0.4) answer2 = brentq(F, 0.8, 1.2) print(answer1, answer2) ``` The function definition shown is just the one from part A1 of the previous exercise. From a plot of that function (you know how to do this) you can see that the plot crosses zero somewhere between x=0 and x=0.4, then crosses again between x=0.8 and x=1.2. You can probably make better guesses than that, but there is no need: any reasonable guesses will work as long as there is one root between the two guesses you give brentq(). If you run this code, you'll probably recognize the results; they're answers to part B of the previous exercise. You can use this technique for solving many numeric problems: you just have to write a python function that returns zero at the desired solution. For example, to solve ``sin(x) = exp(-x)`` you would write the function like this: ```python def something(x): return sin(x) - exp(-x) print brentq(something, 0, pi/2.0) ``` Note that the value returned by the function is just the difference between the right and left sides of the equation, which is zero for a solution. Brentq() has more options: if you want to learn more google "Python brentq", but the simple options shown here work for 99% of problems. ### Assignment Redo exercise 3B, with the following parameters: 1. R = 1 m, other parameters unchanged. 2. R = 5 m, other parameters unchanged. 3. R = 5 m, h = 3.2 m, v = 6.8 m. ## Exercise 5: Euler's method You've probably noticed that you do a lot of "simplified" problems in your physics class. We assume air resistance is negligible in free-fall, for example, or pretend that there's no friction in the pulley. The reason we do this is that otherwise the problem is too hard to solve analytically, and it's worthwhile to find an approximate solution anyway. With a computer, though, we can often get numeric answers even if we can't get analytic answers. This exercise will show you an easy computational method to find approximate solutions to Ordinary Differential Equations. We'll start by finding a solution to a problem for which we already knonw the answer: free-fall with negligible air resistance. But first, there's another bit of Python to learn: How to do things more than once. ### For Loops If you know how many times you want to do something in Python, the "for loop" is the best way to make Python do that something that many times. (If you aren't sure how many times it'll take, then the "while loop" is a better choice, but let's just do the for loop for now.) Try this code in a new text file in Python: ```python for j in range(10): print("j =", j, "j**2 =", j**2) ``` Here's how that code works: It sets the value of $j$ to be each item in the list "range(10)", then does all the following indented lines once for each value of $j$. So "range(10)" produces a list that goes from 0 to 9 (ten elements total), and then goes through the indented print statements 10 times. The first time through, $j=0$. The second, $j=1$; and so on through $j=9$. One important thing to note is that Python always starts counting with 0. Keep this in mind, because it causes trouble occasionally! For example, if you wanted $j$ to range from 1 to 10 (instead of 0 to 9) you'd have to change the code a bit. This would work, and if you stare at it long enough it might make sense: ```python for j in range(1,11): print("j =", j, "j**2 =", j**2) ``` Another important thing to note is that there's nothing special about the name $j$. I could just as well call it $Rumplestiltskin$, it'd still work. A third important thing to note is that the for loop does't have to use range(): it can use any list of things. ```python listOthings=(1,2,3,42,17,"Hike!") for QB in listOthings: print(QB) ``` Of course, you can't print the square of 'Hike!', but you get the idea. Ok, got for loops? Now let's go back to doing physics. ### Euler's Method Let's start with the differential equations for free-fall: $$v = \frac{dx}{dt} \qquad a = \frac{dv}{dt}$$ where $a$ is a constant. We can rearrange these equations: $$dx = v \; dt \qquad dv = a \; dt$$ Remember what $dx$ and $dt$ and $dv$ mean: these are the infinitesimal changes in position, time, and velocity. This form of the equations suggests a possible way of approximating the solution to the differential equations: take little time-steps $dt$, and calculate the new $x$ and $v$ for each time-step. $$x_{new} = x_{old} + v \; dt \qquad v_{new} = v_{old} + a \; dt$$ This is called "Euler's method". We have to take really small steps to get a good solution, which means we have to take lots of steps to get anywhere. That's what the computer is for: doing all those boring calculations! Here's some example code. Don't try putting this into a Python console one line at a time: put it in a new text file, save it, and run it as a program. ```python from pylab import * g = 9.8 # We are on Earth. dt = 0.1 # 1/10 second time step N = 100 # I'd like to see 100 points in the answer array vi = 50.0 # initial velocity xi = 25.0 # initial position # first, set up variables and almost-empty arrays to hold our answers: t = 0 x = xi v = vi time = array([0]) # initial value of time is zero! height = array([xi]) # initial height is xi velocity = array([vi]) # initial velocity is vi # Note that 't', 'x' and 'v' are the current time, position and velocity, but # 'time', 'height' and 'velocity' are arrays that will contain all the positions # and velocities at all values of 'time'. # Now let's use Euler's method to find the position and velocity. for j in range(N): # here are the calculations: t = t + dt x = x + v * dt v = v - g * dt # And here we put those calculations in the arrays to plot later: time = append(time,t) height = append(height,x) velocity = append(velocity,v) # Now that the calculations are done, plot the position: plot(time, height, 'ro') xlabel("Time (s)") ylabel("Height (m)") # just for comparison, I'll also plot the known solution! plot(time, xi + vi*time - 0.5*g*time**2, 'b-') show() ``` You can see on the graph below that this actually works pretty well. It's a little bit off, but not bad; you can make it better by making dt smaller and N correspondingly larger.  ### Assignment Change that sample code to 'solve' Simple Harmonic Motion instead of free-fall. The equations for SHM, in Euler format, are $$dx = v \; dt \qquad dv = a \; dt$$ but instead of $a=-g$, use $a = -\omega x$. Use $x_i = 1$, $v_i=0$, and $\omega = 1$. Plot the position for 10 seconds. Make the code put your name in the title, of course, and '$A$' students will plot the known solution also! ## Exercise 6: Linear Curve Fits It is frequently the case in introductory physics labs (and in Advanced Physics Labs!) that we have a data set and a theoretical model for the data, and we would like to find the model parameters that make the model best match the data set. For example, a student in my lab was recently working on a sensor that was supposed to track the angular position of a rotating magnet, but it seemed as if the sensor was miscalibrated. It was difficult to tell for sure, though, since the the sensor (if it was working correctly) was reporting position in units of $\frac{1}{1024}^{ths}$ of a revolution, and it wasn't off by much if it was off at all. Moving the magnet by hand was not precise enough to determine whether the miscalibration was _real_, and if it was real whether it was _consistent_. To investigate the system further, he recorded the magnet position reported by the sensor each rotation, for 10 rotations. He could only "eyeball" the rotations, so his data was not exact but he was certainly within 2 degrees (5.7 sensor units) of 'zero' each time. His data is below. | Turns | position | | ----- | -------- | | 0 | 0 | | 1 | -2 | | 2 | -11 | | 3 | -15 | | 4 | -24 | | 5 | -32 | | 6 | -40 | | 7 | -50 | | 8 | -52 | | 9 | -62 | | 10 | -65 | Based on this data, could we say that the sensor was miscalibrated, or was it just random error? If it was miscalibrated, by how _much_ was it miscalibrated? In other words, if the sensor was supposed to put out 1024 pulses per rotation, how many pulses was it _really_ putting out per rotation? To answer this question, plot the data. Are the data points scattered around the average, or is there a definite trend? Is the trend _linear_? If the data is linear, what is the slope? How many pulses per revolution was the sensor really generating? Put the data into two Python _lists_: ```python turns = [0, 1, 2, 3, ... 10] position = [0, -2, -11, ... -65] ``` The ``polyfit()`` function in python[^2] does curve fitting of polynomials, including first-order (linear) polynomials. ### Assignment Look up "python polyfit()" (Google is your friend, here!) and figure out how to use ``polyfit()`` to find the best-fit line for this data. Report the slope and intercept, and whether the data seems to fit the theory that the sensor is miscalibrated. Answer the questions above. '_A_' students will plot the data, with the best-fit line, as part of their answer.