Introduction to Computational Physics Lab

name: ____________

This is a Jupyter Notebook. It is organized in "cells", some containing Python code (those labelled with []:) and others containing text (such as this one).

When you start this file with Jupyter, you can edit and execute the code blocks. Click in a code block, type or edit some code, then type Shift-Enter to execute the code. You can also edit text blocks like this one; type Shift-Enter to render the text after editing. If the address of this page in your browser ends in .html then you are viewing a static rendering of the notebook and you will not be able to edit it.

To complete this lab, you must

1. Get Jupyter running. See How to Get Started with the Jupyter Notebook.
2. Download the code files and save them in a folder where you plan to do your computational work this semester. The files may be found on the lab web page.
3. Type your name in the text block at the top of this document. (Double-click the block to open it for editing.)
4. Read through this entire notebook to see how it's organized.
5. Execute the import lines without error.
6. Answer every "Question" in B. Exercises by writing and executing the appropriate code. That section starts with "Example 1" and "Question 1". You can answer the "Question 1" by following the model in the "Example 1". Repeat for all Questions.
7. Section C. Importing Packages is optional.
8. Study the methods illustrated in D. Vector Sum Example.
9. Answer all the E. Vector Problems as directed by writing and executing the code required.
10. Save this notebook as a .ipynb file and submit it on the course web page when the is lab due.

---

A. Set Up

To be able to execute much of the code we use, you will need to first execute these two import statements.

from numpy import *
import scipy.constants

The NumPy package (short for numerical python) will allow you to use vectors, trig functions and much more.

The scipy.constants package contains a lot of constants and unit conversion values. The SciPy package (short for scientific python) contains many advanced computational physics libraries.

B. Exercises

Below are a series of examples, each followed by a similar question. Copy the example code and paste it into the question's answer cell. Then edit as needed to answer the question.

Example 1

Mt. Everest is 29 000 ft high. How many meters is that?

h_ft = 29000
h_m = h_ft * scipy.constants.foot
print(f'height is {h_m} meters')

height is 8839.199999999999 meters

Note: The scipy.constants conversion factors all convert to SI units.

Question 1

The Capitol is 22 800 ft from Truax. How many meters is that?

Example 2

The SI unit of atmospheric pressure is the Pascal (Pa). One Pascal is one Newton of force per square meter of area:

$$1\;\rm{Pa} = 1 \frac{\rm{N}}{\rm{m}^2}$$

Standard Atmospheric Pressure is $1.013 \times 10^5\;\rm{Pa}$. Express this in pounds per square inch.

Notes:

• to convert from (not to) SI units, you have to divide by the conversion factor.
• to square a quantity, python uses a double-asterisk, as in x**2

p = 101300                                   # p is atmospheric pressure in N per m²
p = p / scipy.constants.lbf                  # p is now in lb per m²
p = p * scipy.constants.foot**2              #             lb per ft²
p = p / 12**2                                #             lb per in²
print(f'Standard Atmospheric Pressure is {p} lb/in²')

Standard Atmospheric Pressure is 14.692322832070191 lb/in²

Question 2

A tire needs to be inflated to a pressure of $2.07 \times 10^5$ Pa. What is this pressure in lb/in²?

Hint: enter numbers in scientific notation with an e: $2.07 \times 10^5$ = 2.07e5.

Example 3

You measure the diameter of a sphere to be $4.02 \pm 0.02$ cm.

Find the sphere's volume in cm$^3$ along with its absolute uncertainty.

d = 4.02                # diameter, cm
unc_d = 0.02/d          # relative uncertainty in d
r = 0.5*d               # find radius
unc_r = unc_d
V = 4/3*pi * r**3       # use ** to raise to a power
unc_V = 3*unc_r
unc_V = unc_V * V       # convert to absolute uncertainty 
print(f'volume is {V} ± {unc_V} cm^3')

volume is 34.01549392577856 ± 0.5076939391907248 cm^3

You can control the number of digits after the decimal that print out with string formatting as follows:

print(f'volume is {V:.2f} ± {unc_V:.2f} cm^3')

volume is 34.02 ± 0.51 cm^3

The :.2f directs python to show 2 decimal places of a "floating point" (ie, decimal) value. If you would like to use scientific notation, use an e instead of f as follows:

print(f'volume is {V:.2e} ± {unc_V:.2e} cm^3')

volume is 3.40e+01 ± 5.08e-01 cm^3

You can learn more about fancier printing in Python by a web search.

Question 3

Astronomers report the diameter of a planet to be $(6.37 \pm 0.02) \times 10^6$ m.

Find the planet's volume in m$^3$ along with its absolute uncertainty.

Example 4

The coordinates of point $P$ are measured to be $x = 2.55 \pm 0.02$ m, $y = 1.24 \pm 0.04$ m. What is the angle $\theta$ (in degrees) and its uncertainty?

Finding $\theta$ requires the inverse tangent function (also known as the arctangent) as follows.

x  = 2.55     # m
dx = 0.02     # m, absolute uncertainty
y  = 1.24     # m
dy = 0.04     # m, absolute uncertainty
Q = arctan(y/x)
print(f'theta = {Q*180/pi} degrees')

theta = 25.93247363696144 degrees

Computers always work in radians, so the last line converts the radians to degrees (*180/pi).

We have no rule about propagating the uncertainty through an arctan function. But we can find the range of values we would get when $x$ and $y$ are at the edge of their uncertainty.

print(f'the angle is theta = {Q*180/pi} degrees')
print(f'other theta values within our uncertainty are...')
print( arctan((y+dy)/(x+dx)) *180/pi )
print( arctan((y+dy)/(x-dx)) *180/pi )
print( arctan((y-dy)/(x+dx)) *180/pi )
print( arctan((y-dy)/(x-dx)) *180/pi )

the angle is theta = 25.93247363696144 degrees
other theta values within our uncertainty are...
26.47580548677303
26.836164829561802
25.02909705182668
25.375373720723683

So a good measure of the uncertainty in $\theta$ is 1/2 of the difference $\max (\theta) - \min (\theta)$. Looking at the four values printed above, you need to find the maximum and minimum.

Q1 = arctan((y+dy)/(x-dx))*180/pi 
Q2 = arctan((y-dy)/(x+dx))*180/pi
dQ = 0.5*( Q1 - Q2 )
print(f'theta = {Q*180/pi:.2f} ± {dQ:.2f} degrees')

theta = 25.93 ± 0.90 degrees

We've answered the question, but in case you're curious, a better way to find the min and max would be to put the results in a Python list (comma separated between square brackets) and then use the min and max functions:

Q_list = [ arctan((y+dy)/(x+dx)) *180/pi,
           arctan((y+dy)/(x-dx)) *180/pi,
           arctan((y-dy)/(x+dx)) *180/pi,
           arctan((y-dy)/(x-dx)) *180/pi ]
dQ  =   0.5*( max(Q_list) - min(Q_list) )
Q_ave = 0.5*( max(Q_list) + min(Q_list) )
print(f'theta = {Q*180/pi:.2f} ± {dQ:.2f} degrees')

theta = 25.93 ± 0.90 degrees

Question 4

A distance along a slope is measured to be $d=2.00 \pm 0.02$ m. The slope lies at an angle of $\theta = 30.0 \pm 0.1$ degrees above the horizontal. The vertical height of this slope is given by

$$y = d \sin \theta$$

Find the value and uncertainty of $y$ (ie, find $y \pm dy$).

C. (Optional) Importing packages

If you have a computational problem to solve, there's a good chance someone has already solved it. Python's FLOSS (free-libre open source software) philosophy lends itself to community-shared solutions.

For example, you might wonder if anyone has written code to automatically handle uncertainty calculations. Answering this kind of question is best done with an internet search for, in this case, "python uncertainty". You would quickly learn that the uncertainties package is popular.

However, you cannot import it since it does not come with the standard Python/Jupyter installation. You will receive an import error if you execute the import line below before installing the package. In this case (and for many Python packages), installation is easy. Execute the following at a command line:

pip install uncertainties

Note: You can obtain a command line in Windows by running app cmd, in a Mac by launching a terminal. This is not available in jupyter.org/try.

Now you may import and use the package. Check its documentation for how to use it.

import uncertainties

Now uncertainty calculations can be done automatically (and more accurately)!

Consider Example 3, above, using this package.

d = uncertainties.ufloat(3.14, 0.02)  # from Example 3, above
r = 0.5*d
V = 4/3*pi * r**3
print(f'volume is {V} cm^3')

volume is 16.21+/-0.31 cm^3

D. Vectors in Python

Suppose we know the vector sum $\vec{F_1} +\vec{F_2} + \vec{F_3} = \vec{0}$, and we are given the following information.

• vector $\vec{F_1}$ lies along the positve $x$ axis, with magnitude 200 N
• vector $\vec{F_2}$ makes an angle of $107^{\circ}$ with the $+x$ axis, and has magnitude 300 N

Answer the following questions:

1. What is the direction and magnitude of vector $\vec{F_{4}} = \vec{F_1} + \vec{F_2}$ ?
2. What must be the direction and magnitude of vector $\vec{F_3}$ ?

---

Solution

First we need to define vectors $\vec{F_1}$ and $\vec{F_2}$. In Python, this is done with the array statement. Here are two equivalent ways we could define $\vec{F_1}$ :

F1 = array(( 200, 0 ))       # specify the components of F1
F1 = 200 * array(( 1, 0 ))   # write F1 as a scalar times a unit vector

For $\vec{F_2}$ it's easier to use the unit vector form. All angles must be in radians.

Q = 107*pi/180     # angle in radians 
F2 = 300 * array(( cos(Q), sin(Q) ))

If you add the variables F1 and F2, it will (correctly) add the components of the vectors. Let's add them and call the sum $\vec{F_{4}}$:

F4 = F1 + F2
print('F4 = ',F4)

F4 =  [112.28848858 286.89142679]

This means that $\vec{F_{4}} = 112 \hat{\imath} + 287 \hat{\jmath}$ N.

Let's find the magnitude of this vector: $F = \sqrt{ F_x^2 + F_y^2}$.

• to square a Python array use F4**2; this will square each component individually
• sum will add up all the (squared) array components
• sqrt is the square root

mag = sqrt(sum(F4**2))
print(f"magnitude of F4 is {mag} N")

magnitude of F4 is 308.08342284724057 N

To find the direction of $\vec{F}$, we use $\theta = \arctan \left (\frac{F_y}{F_x} \right )$.

The components of a python array are accessed by indexing them, starting from zero, like this:

$F_x =$ F4[0]

$F_y =$ F4[1]

angle = arctan( F4[1] / F4[0] )
print(f'angle is {angle*180/pi} degrees')  # convert rad to degrees

angle is 68.62476644579117 degrees

The vector $\vec{F_3}$ should be in the opposite direction to $\vec{F_{4}}$. So it should have the same magnitude with an angle $\theta+180^{\circ}$ or $\theta+\pi$ radians.

F3 = mag * array(( cos(angle+pi), sin(angle+pi) ))
print("F3 =",F3)

F3 = [-112.28848858 -286.89142679]

We know that $\vec{F}_3 = -\left (\vec{F_1} +\vec{F_2} \right ) = -\vec{F_{4}}$. Look above and note that the components of $\vec{F_3}$ are the negative of the components of $\vec{F_{4}}$.

That's it, we're done. The following code block summarizes the operations required to answer the question.

F1 = 200 * array(( 1, 0 ))
F2 = 300 * array(( cos(107*pi/180), sin(107*pi/180) ))
F4 = F1 + F2
mag = sqrt(sum(F4**2))
angle = arctan(F4[1]/F4[0])
print(f'magnitude of F4 is {mag:0.2f} N and its direction is {angle*180/pi:0.2f} degrees')

magnitude of F4 is 308.08 N and its direction is 68.62 degrees

A Note about arctan

The arctan function above has a problem: it can give only give angles between $-\pi/2$ and $+\pi/2$, in the first and fourth quadrants. This is because the computer cannot distinguish between, for example, arctan(-1/1), which should be $- \pi/4$ and arctan(1/-1), which should be $3 \pi/4$.

The solution is to use the arctan2(y,x) function, giving it the $y$ and $x$ values individually (not their ratio) so it knows which angle to return.

print('arctan(1/-1)=',arctan(1/-1)*180/pi)
print('    wrong quatrant!')
print('arctan2(1,-1)=',arctan2(1,-1)*180/pi)
print('    correct!')

arctan(1/-1)= -45.0
    wrong quatrant!
arctan2(1,-1)= 135.0
    correct!

Note the arctan2 function requires y,x, not y/x.

E. Vector Problems

Now it's your turn.

Find the magnitude and direction of the following vectors. Give your answers as follows:

• all values are decimals expressed with three significant figures
• directions are angles in degrees off the $+x$ axis (postive angles for counterclockwise and negative angles for clockwise)

(a) $\vec{A} = - \frac{1}{\sqrt{2}} \hat{\imath} - \frac{1}{\sqrt{2}} \hat{\jmath}$

(b) $\vec{B} = - 10 \hat{\imath} + 10 \hat{\jmath}$

(c) $\vec{C} = 1000 \hat{\imath} - \hat{\jmath}$

(d) $10\vec{A} - \frac{\vec{B}}{2}$

Hint: for question (d), you can write 10*A - B/2 (assuming you defined A and B correctly).