Fall 2013

Office hours

Introduction

I assume in this course that you have taken Physics 265 or its equivalent, and have some experience programming in python. I will not assume, however, that you remember everything you did in that course, so please do not panic! I will begin with a very brief review of simple python, and pretty quickly we will get down to numerics and then physics.

You can run python from a shell (interactive terminal) by typing

ipython

or

python

either of which will give you an interactive python interpreter. Using python interactively can be helpful when you're debugging, or if you want to quickly visualize something. At the python interpreter, you can do simple math

>>> 1 + 3
4
>>> 3*4 - 1
11
>>> x = 5
>>> y = 3
>>> x + y
8

The numpy and matplotlib modules

Will be extensively using the numpy_ and matplotlib_ modules. Python has several ways to import modules. From Physics 265, you should already be familiar with

from visual import *

which allows you to use visual python. In this class, we will instead use

import numpy # for matrix stuff
import pylab # for plotting

which imports the numpy and matplotlib modules. The import style import numpy differs from from numpy import * in that it requires you to prefix any "numpy" functions with numpy.. This makes it somewhat more cumbersome to type, but ensures that you always know which module provides the functions that you are using. When we used only one module, of course, this wasn't much of an issue. Henceforth, I will always assume that you have imported these two packages.

# the following creates a 1-D array of real numbers going from 0 to 10
# by steps of 0.1
x = numpy.arange(0,10,0.1)

# The following defines each element of y to be the cosine of the
# corresponding element of x
y = numpy.cos(x)

# This tells pylab what you want to plot
pylab.plot(x, y)

# Always label your axes!
pylab.xlabel('time (s)')
pylab.ylabel('voltage (mV)')

# with savefig you can save the plot as a file
pylab.savefig('sin-wave.pdf')
# and with show you can view the plot interactively
#pylab.show()

There is a lot more that you can do with matplotlib and numpy, but when you need to do more, you can figure it out using their webpages. Of particular interest is the matplotlib gallery, which can be a handy resource for figuring out how to do a new kind of plot.

Pair programming

In this course, you will be doing pair programming, which means you will be working in pairs at a single computer. One member of the pair will have the keyboard and will be the driver while the other is the navigator. The role of the navigator is to watch for mistakes, think ahead, give advice, and plan strategically. I will periodically ask you to switch roles. From time to time, your pairs will also be switched up, possibly in the middle of a project, so you need to be careful to always understand how your program works, and to write it in such a way as it will be easily understood by a new partner!

Learning goals

The goal of this course is to teach you to use computational methods to solve real-world physics problems. At the same time, we aim to help improve and solidify the material you are learning in the paradigms. The following is

Numerical methods:

  1. Verlet's method of integration
  2. Quadrature methods such as trapezoidal method or midpoint method
  3. Fast Fouriet transform (may be skipped)
  4. Quadrature in multiple dimensions and curvilinear coordinates

Programming in python:

  1. Plotting 1D plots with matplotlib
  2. Using numpy arrays
  3. Writing functions
  4. Performing FFTs (time permitting)

Electrostatics

Plotting functions in three dimensions (four point charges)

To begin with electrostatics, we are going to simply work on plotting the electrostatic potential for a collection of point charges. As you have seen in the paradigm, there are various ways you can visualize the potential. You will write a function that computes the potential given $x$ and $y$, and then plot this function in various ways that seem good to you.

The picture to the right is to give you an idea of what a plot of a dipole might look like. You can do better than this.

Integrating to find the electrostatic potential

In the paradigms, you've been learning to integrate in order to find the electrostatic potential due to a charge distribution. At some point (past or future) you will come across integrals that you can't perform, or that lead to special functions that are unfamiliar to you. Doing integrals analytically is useful because it allows you to examine limiting cases analytically. But often you only want to know the answer, or obtain a plot, and for those purposes a numerical integral is perfectly adequate.

In the first project you learned to perform quadrature. In this project we will apply those techniques to multidimensional integrals in different coordinate systems.

Electrostatic potential of a square of charge

Let's begin by computing the electrostatic potential due to a square sheet of charge with side $d$ and charge density $\sigma$. Write a function to compute the electrostatic potential at an arbitrary point in space.

Once you have written the above function, use it to plot the electrostatic potential versus position in the three cartesian directions. If you have extra time, try plotting along an arbitrary third direction. What happens when you change the number of grid points used to perform the integral?

Electrostatic potential of a solid cylinder of charge

Now let's compute the electrostatic potential due to a solid cylinder of charge with length $L$, radius $R$ and charge density $\rho$. Once again, write a function to compute the electrostatic potential at an arbitrary point in space. When you have written this function, plot the potential.

Electrostatic potential of a solid sphere of charge

To make use of the third coordinate system, compute the electrostatic potential due to a solid cylinder of charge with radius $R$ and charge density $\rho$. Write a function to compute the electrostatic potential at an arbitrary point in space, and use it to plot the potential. Does it behave as you expect?

To the right is the potential of a solid sphere, so your plot ought to look something like it.

Electric field due to a large circular sheet of charge

Write a function to compute the electric field due to a large circular sheet of charge. Plotting the electric field is a bit more cumbersome, since it is a vector. Try plotting the field in some symmetry directions where the field is pointing in a known direction. You can also plot multiple components of the electric field on a single plot. Try plotting $\mathbf{E}$ versus $z$ for several different radii, and see how they compare when $z \ll R$.

Electric field due to solid hemisphere of charge

Write a function to compute the electric field due to a solid hemisphere of charge. Different groups will be asked to write their functions using either cylindrical or spherical coordinates. As usual, plot your result.

A pendulum with finite amplitude

Pendulum in the time domain

The first topic we will study will be a pendulum. $\begin{align} V = mgh \end{align}$

In groups of two or four, I'd like you to work out the torque on a pendulum as a function of its angle. This will allow us to compute the angle as a function of time using the differential equation

$ \begin{align} ML^2 \frac{d^2 \theta}{d t^2} = \tau \end{align} $

As you may recall, Verlet's algorithm consists of taking the finite difference formula for the second derivative and solve for the value at a time in the future. In this case we're looking at the angle: $\begin{align} \frac{d^2 \theta}{d t^2} &= \frac{\tau}{ML^2 } \\ &= \lim_{\Delta t \rightarrow 0} \frac{\theta(t+\Delta t) + \theta(t-\Delta t) - 2\theta(t)}{\Delta t^2} \\ &\approx \frac{\theta(t+\Delta t) + \theta(t-\Delta t) - 2\theta(t)}{\Delta t^2} \end{align}$ Solving for $\theta(t+\Delta t)$, we find that \begin{align} \theta(t+\Delta t) &\approx 2\theta(t) - \theta(t-\Delta t) + \frac{\Delta t^2}{ML^2 }\tau \end{align} In Physics 265, we solved this by keeping track of the ball.pos and ball.oldpos. This is convenient, in that we are able to compute the position at any given time while only storing a couple of positions, but has the disadvantage that we forget where the ball was, and thus lose the opportunity to plot its position versus time.

In this first project, you will store the angle as a function of time for a simple harmonic oscillator, and then plot it. You should also work out the prediction from the small angle approximation, and plot the two together, to see the effect of the finite angle.

Pendulum in the angular domain

In the previous activity, we computed the angle of a pendulum as a function of time. We will now use the first integral of the motion to compute the time as a function of angle. This is a fancy way of saying that we will use energy conservation to find the angular velocity as a function of angle and then integrate that.

The kinetic energy of a pendulum is just $\frac12 I\omega^2$, so if we know the energy of a pendulum, we can easily solve for $\omega$, which is of course the same thing as $\frac{d\theta}{dt}$. You might think that we could integrate $\frac{d\theta}{dt}$ to find $\theta(t)$, but the problem is that we know $\frac{d\theta}{dt}$ as a function of $\theta$, not as a function of $t$. Fortunately, we can integrate the inverse of this derivative, which will give us $t(\theta)$.

First, solve for $\omega(\theta)$ for a grid of $\theta$ values. Then you should invert this, so you have $\frac{dt}{d\theta}$ on a grid. Finally, you will need to integrate this function numerically, to find $t(\theta)$.

Once you have finished computing $t(\theta)$ (which may be after the next section), you will want to plot your Verlet's solution computed in the time domain against the solution computed using the first integral of the motion to be sure the answers agree. As long as you're at it, you may as well also plot the solution found using the small-angle approximation. It will be interesting to compare these solutions both for small maximum angle (e.g. $\theta = 0.1$) and for large maximum angle (e.g. $\theta = 3$).

Quadrature

There are a number of ways we can integrate numerically. In the common case, you have a function whose value you know on a regularly spaced grid of points, and want to know its integrall over that range. The simplest approach is to just add everything up $ \begin{equation*} \int_a^b f(x)dx \approx \sum_x f(x)\Delta x \end{equation*} $

This is a reasonable approach, and I recommend that you start with it. However, there are some improvements we can make, based on the locations of the starting and ending points. Often your grid points will start and end on $a$ and $b$. In that case, the above approach over-estimates the importance of those two points. Imagine if you had just those two endpoints, and your function was constant: you'd get the integral off by a factor of two! In fact, with very little work, we can construct an approach that gives the exact answer for any straight line, just by working out the area under

Above is a graphical construction for the derivation of the trapezoidal rule. The integral is approximated by the area under all the trapezoids constructed by connecting the sampled points with straight lines.

You might think that if the trapezoidal rule provides exact answers for any straight line, we could do better with a rule that gives exact answers for any parabola. The approach that works for this is called Simpson's rule, and is based on the idea of interpolating every pair of intervals with a parabola.

Above is a graphical construction for the derivation of the Simpson's rule. Each pair of intervals is interpolated with a parabola passing through three points.

The trapezoidal method is what is called a first-order method, and for functions that are well-approximated by a polynomial the error is proportional to $\Delta x$ and linear functions are treated exactly. Simpson's method, in contrast, is a second-order method, and thus errors are proportional to $\Delta x^2$, and parabolas are handled exactly.

Write two functions to integrate a function over a given range using the trapezoidal method or Simpson's rule (both of which you will need to work out geometrically). Your function should have the following type:

def trapezoidal(f, a, b, Nsteps):
    """This function integrates the function f from a to b, using N
    steps in between."""
    # insert actual code here.
# You should be able to call your function with something like:
trapezoidal(math.sin, 0, math.pi, 5)

Once you have written your trapezoidal and Simpson's quadrature routines, you can test them by integrating a function whose integral you know analytically, and seeing how the error changes as you increase the number of points (or equivalently decrease $\Delta x$).

Fast Fourier Transform (FFT)

Finally, let's analyze your output using a fast Fourier transform (FFT). The FFT is an incredibly useful algorithm that performs a Fourier transform on discretly sampled data very quickly. The scaling of an FFT is $O(N\log N)$, where the scaling for a naive Fourier transform is $O(N^2)$. The difference between $O(N^2)$ and $O(N\log N)$ can be huge as $N$ becomes large.

In numpy, we can perform an FFT by simply

t = numpy.arange(0,10,0.125)
y = numpy.cos(5*numpy.pi*t) + numpy.sin(3*numpy.pi*t)
ytransformed = numpy.fft.fft(y)
omega = numpy.zeros_like(t) # there are the same number of frequencies as times
omega_min = 2*numpy.pi/10 # lowest possible frequency observable
for i in range(len(omega)):
    if i < len(omega)/2:
        omega[i] = i*omega_min
    else:
        omega[i] = (i - len(omega))*omega_min
pylab.plot(omega, ytransformed.real, 'o-r', label='real')
pylab.plot(omega, ytransformed.imag, 'o-b', label='imaginary')
pylab.legend() # use a legend!
pylab.xlabel('angular frequency')
pylab.ylabel('fourier transform of two sinusoids')
pylab.savefig('fft.svg')
pylab.show()

To the right is the fast Fourier transform computed in the example.

As you can see from the code, the frequency is returned in a somewhat unusual order: negative frequencies come in the second half of the array. This is the standard behavior of any FFT function, so it's worth getting used to. Note that the output of an FFT is a complex array.

Take your data for $\theta(t)$ and FFT it, and see what frequencies show up.