As a computational physicist (in training) you should be able to:
We are not doing any of this (but the principles remain the same):
Lectures:
Exercises
Project
Support
Grading
We also assume you know how to code to some degree, so we will focus on learning to use Python for computational physics rather than an introduction to Python in general.
In 1888 Johannes Rydberg published his famous equation describing the wavelengths of light emitted from a hydrogen atom:
$\frac{1}{\lambda}=R\left(\frac{1}{m^2}-\frac{1}{n^2}\right)$
where $R$ is the Rydberg constant $R=1.097 \times 10^{-2}$nm$^{-1}$, and $m$ and $n$ are positive integers. For a given value of $m$, the wavelengths $\lambda$ given by this formula for all $n \gt m$ form a series. The first three of these series, for $m=1,2,3$, are named after their discoverers Lyman, Balmer and Paschen respectively.
Let's make a simple Python code that solves this equation for $m=1,2,3$ and outputs the results.
R = 1.097e-2
for m in [1,2,3]:
print("Series for m =",m)
for k in [1,2,3,4,5]:
n = m + k
invlambda = R*(1/m**2-1/n**2)
print(" ",1/invlambda,"nm")
The agreement with experiment is actually remarkably good, with the differences due to small details of the electronic structure (spin, relativistic corrections, hyperfine splitting).
As a test that you can actually get Python code running on your system, try to re-run the Rydberg code yourself. You have several options for running Python:
python
at the command prompt and it should give you something like this:Python 3.7.4 (default, Aug 13 2019, 15:17:50)
[Clang 4.0.1 (tags/RELEASE_401/final)] :: Anaconda, Inc. on darwin
Type "help", "copyright", "credits" or "license" for more information.
>>>
then you can just type your code or paste it from a text file.
idle
,spyder
or any other Python flavoured development environment you prefer.jupyter-lab
, a recent update of the jupyter-notebook
environment.
for m in [1,2,3]
R = 1.097e-2
wavefunction = 1 + 2j
x = "Use sensible variable names!"
test_input=input("Enter the number of programming languages you know:")
Then the variable is read as a string whatever you enter. This is generally not ideal for a computational physics code:
print(test_input*2) # correct answer should be 2 x the number of languages obviously
print(test_input*1.5)
This can be handled easily by just enforcing the correct type after reading (this can be done directly with the input command as well):
type_input = float(test_input)
print(type_input*1.5)
There are some useful shortcuts for the common variable operations:
mod_var = 10
mod_var += 5
print(mod_var)
mod_var -= 10
mod_var *= 5
print(mod_var)
mod_var //= 7
print(mod_var)
mod_var = 3
eric,john,terry,michael,graham = mod_var,mod_var*10,mod_var**2,mod_var/4,mod_var**0.5
print(eric,john,terry,michael,graham)
Defining multiple variables at once makes the code more difficult to read, but can be useful at times. Note that everything on the right is resolved first.
mod_var = 3
eric,john,terry,michael,mod_var = mod_var,mod_var*10,mod_var**2,mod_var/4,mod_var**0.5 # notice the change in the 5th variable, demonstrating the timing of variable evaluation
print(eric,john,terry,michael,graham)
While many useful mathematical functions can usually be evaluated with a combination of standard Python commands, it is often easier and computationally more efficient to use imported functions from a package. Lists of common packages and the functions they contain can be easily obtained from Python documentation, but math
is an obviously useful beginning.
from math import sqrt
fun_var = 10**0.5/sqrt(10)
print(fun_var)
While you can import all functions from a given package:
from math import *
import math
x = math.sqrt(10)
this is generally dangerous, as you may be including duplicate functions you are unaware of. Much better to begin your code with an import of everything you will actually use. This is also generally easier to read than using modules.
from math import log,cos,sqrt,pi
bored_var = int(input("Enter your level of interest in this course (1-10): "))
if bored_var>10:
print("You are lying.")
bored_var = cos(pi)
print(bored_var)
elif bored_var<1:
print("You are dying.")
bored_var = sqrt(pi)
print(bored_var)
else:
while bored_var<10:
print("Be happy!")
bored_var = bored_var + log(10)
print(bored_var)
Once we move into handling real data from a physical problem, it is hardly practical to assign every value to a different variable, and the use of containers becomes important. Python has both lists and arrays, and they are suited to different tasks commonly encountered in computational physics.
int
,float
.Let's define a simple list:
simple_list = list(range(1,11))
print(simple_list)
This is just a list of integers, produced using the range
function and converted to a list. Then there a several functions that work directly on lists:
print(sum(simple_list)) # sum of the elements
print(len(simple_list)) # number of elements
print(max(simple_list)) # max of elements
print(min(simple_list)) # min of elements
print(sum(simple_list)/len(simple_list)) # mean of elements
Then you can also use the map
function to apply any function to a list:
from math import exp
exp_list = list(map(exp,simple_list))
print(exp_list)
We can also add elements to and remove elements from a list (note that the list must exist for this to work, but you can create an empty list empty_list = []
):
from math import exp # we don't need this everytime in the notebook, but as a reminder of good coding practice...
exp_list.append(exp(11))
print(exp_list)
exp_list.pop(0) # first element has an ID of zero.
print(exp_list)
We will see later how we can use list comprehensions as an elegant way to define and create lists based on existing lists.
For most problems in physics, arrays are likely to be more useful as we are commonly dealing with well-defined data sizes and physical quantities e.g. vectors. In Python the tools for manipulating arrays are contained within the numpy package.
from numpy import zeros,ones,empty
zero_mat = zeros([3,3],float)
one_mat = ones([3,3],float)
empty_mat = empty([3,3],float)
print(zero_mat)
print(one_mat)
print(empty_mat)
Note that empty
does not create an empty array, it just uses whatever numbers it has in memory. Both the zeros
and ones
functions actually use time to create, so if you just need an empty array ready to be filled, empty
is more efficient. But remember it will not actually be empty...
sub_mat = one_mat - empty_mat
print(sub_mat)
Also note that Python in general allows for the use of arrays as if they were normal variables, but there are some cases where it behaves differently - usually it is trying to help you, but you should be aware of it. Think about what the contents of these two arrays should be from reading the code:
from numpy import array
array_a = array([1,1],int)
array_b = array_a
array_a[0] = 2
Then let's look at what actually happens:
from numpy import array
array_a = array([1,1],int)
array_b = array_a
array_a[0] = 2
print(array_a)
print(array_b)
This is probably not what you were expecting. Common coding understanding (unless you come from C/C++) would assume that array_b = array_a
creates a new array array_b
with the same elements as array_a
, and then only one element of array_a
is changed. However, Python does not want to waste time copying potentially huge data arrays around and this actually just means array_b
also points to the same data storage as array_a
and they both are changed. If you really want a new array with the same values, then you need to do this:
from numpy import copy
array_a = array([1,1],int)
array_b = copy(array_a)
array_a[0] = 2
print(array_a)
print(array_b)
Most often you will want to read in the values of an array from a external source, such as a text file. Here we can combine some of the previous lists we generated into an array (assuming they have the same length) and use this as our data. We will then write it to a file and read it again...which obviously makes no sense, but demonstrates the functions needed.
from numpy import array,loadtxt,savetxt
print("List 1",simple_list)
print("List 2",exp_list)
output_array = array([simple_list,exp_list],float)
print("Combined array")
print(output_array)
savetxt("array_data.txt.gz",output_array) # save in gzipped format, since the data is so huge...
input_array = loadtxt("array_data.txt.gz",float)
print("Read array")
print(input_array)
Let's try and get our own data, read it in and analyze it using Python. You can use whatever you want, but an obvious source is the FMI. Just pick something sensible and download it in CSV format (comma separated values). If you do have a CSV file rather than a simple TXT, then you will likely need to edit it a little while reading it in. For the weather data, which has this format:
Year,m,d,Time,Time zone,Maximum temperature (degC)
2004,7,15,00:00,UTC,18.5
2004,7,16,00:00,UTC,21.3
2004,7,17,00:00,UTC,20.6
2004,7,18,00:00,UTC,21.9
2004,7,19,00:00,UTC,23.2
2004,7,20,00:00,UTC,20.3
....
something like this works:
from numpy import array,loadtxt
weather_array = loadtxt('weather.txt', skiprows=1, delimiter=',', usecols=(0,5), unpack=True, dtype=None)
Now calculate the number of values in your data, the sum of them and the arithmetic mean...or whatever you thinks makes sense for the data you are using.
#%%timeit # comment out print statements when testing speed. This is an iPython magic function.
from numpy import loadtxt,mean
weather_array = loadtxt('weather.txt', skiprows=1, delimiter=',', usecols=(0,5), unpack=True, dtype=None)
weather_mean = sum(weather_array[1])/len(weather_array[1])
weather_2004 = []
weather_2005 = []
rng = range(len(weather_array[1]))
for n in rng:
if (weather_array[0,n] == 2004):
weather_2004.append(weather_array[1,n])
elif (weather_array[0,n] == 2005):
weather_2005.append(weather_array[1,n])
weather2004_mean = sum(weather_2004)/len(weather_2004)
weather2005_mean = sum(weather_2005)/len(weather_2005)
print("Mean temperature:",weather_mean)
print("2004 mean:",weather2004_mean)
print("2005 mean:",weather2005_mean)
While this is an intuitive (old school perhaps) way to solve the problem, it is not a very pythonic approach to the problem...it can also be done as follows:
#%%timeit # comment out print statements when testing speed
from numpy import loadtxt,mean,where
weather_array = loadtxt('weather.txt', skiprows=1, delimiter=',', usecols=(0,5), unpack=True, dtype=None)
weather_mean = weather_array[1, :].mean()
weather2004_mean = weather_array[1, where(weather_array[0] == 2004)].mean()
weather2005_mean = weather_array[1, where(weather_array[0] == 2005)].mean()
print("Mean temperature:", weather_mean)
print("2004 mean:", weather2004_mean)
print("2005 mean:", weather2005_mean)
The approach you choose should really balance readability and efficiency...but for the latter, you should really only worry about functions that are used heavily.
While packages like numpy
offer a huge variety of optimized functions, in developing your own code it will often be useful to define custom functions and this is very straightforward in Python. Generally, as with importing, it is a good idea to define all the functions you need at the beginning of the code, as this guarantees they are defined before use, makes it easy to change them and is more readable for other users. User-defined functions can be used on lists and arrays as with any standard function, but check carefully they behave as you expect. Also remember that most packages use optimized code (math
functions are mostly written in C) and if someone has already implemented the function you want, theirs is probably faster.
Let's look at an example to find the factors of a given number. These can be found relatively easily by dividing repeatedly by all integers from 2 to n and checking to see if the remainder is zero:
def factors(n):
factorlist = []
for i in range(2,n+1): # note that range stops at n-1
if n%i==0: # % is the modulo function to find the remainder after division.
factorlist.append(i)
return factorlist
# using a list comprehension
def factors2(n):
factorlist = []
factorlist = [n//i for i in range(1,n) if n%i==0] # we are doing the sampling in the opposite direction
return factorlist
print(factors(99))
print(factors2(99))
Then if we want to know if a number is a prime, we can just check for all the numbers that only have one factor, themselves.
for n in range(2,100):
if len(factors(n))==1:
print(n)
# using a list comprehension
primelist = [n for n in range(2,100) if len(factors2(n))==1]
print(primelist)
You will see in the exercises that there are much more efficient ways to solve this problem...which is generally true for most coding problems.
It is already evident for several examples that we have looked at that it would be useful to have ways to visualize the data and for real computational physics problems this is normally the key to actually getting to the physics. Let's start by plotting perhaps the most classic of functions, sin
.
from numpy import sin,linspace # generally these are faster than the in-built Python functions
import matplotlib.pyplot as plt # useful for functions you will use a lot, but be careful that you don't lose track
x = linspace(0,10,100) # create a float array with values from 0 - 10 with 100 intervals
y = sin(x)
plt.plot(x,y,'g--') # plot a dashed line in green
plt.show() # you can do this with numpy as well, but is is clearer and safer to define the functions initially...but many people don't.
This is clearly very basic, but the options with the matplotlib library are very extensive (and there is a wealth on information online explaining them) and we can easily add a lot more information in the same manner.
from numpy import sin,cos,exp,linspace
import matplotlib.pyplot as plt
x = linspace(0,10,100)
y_sin = sin(x)
y_cos = cos(x)
y_exp = 1/exp(x)
plt.figure(figsize=(7.5,5)) # set the figsize
plt.rc('font',size=16) # set the font size
plt.plot(x,y_sin,'b--',label='sin(x)')
plt.plot(x,y_cos,'m--',label='cos(x)') # adds the second plot to the graph, then we plot the whole lot at the end with show()
plt.plot(x,y_exp,'r--',label=r'$\frac{1}{\exp(x)}$') # the 'r' forces this to be read as a raw string and the '\' is ignored by Python.
plt.title('The classics',fontsize=20)
plt.xlabel('x')
plt.ylabel('Function(x)')
plt.legend(fontsize=12)
plt.show() # notice that we don't need to include this in a Jupyter notebook, but it is good practice.
It is also easy to plot different styles of graph, depending on your data.
from numpy import random
from scipy import randn # another package, with many scientific focused functions
import matplotlib.pyplot as plt
rdist1 = 10.0*randn(1000)+20.0 # make a 1000 gaussian-distributed random numbers with a sigma=10 and mean=20
rdist2 = 10.0*random.randn(1000)+20.0
plt.hist(rdist1,color='blue')
plt.hist(rdist2,color='red',alpha=0.5)
plt.show()
plt.subplot(2,1,1) # now plot the two histograms as separate figures on 2 rows with 1 column, and this is row one.
plt.hist(rdist1,color='blue')
plt.subplot(2,1,2)
plt.hist(rdist2,color='red')
plt.show()
from numpy import loadtxt
from scipy import randn
import matplotlib.pyplot as plt
star_data = loadtxt("stars.txt",float)
x = star_data[:,0]
y = star_data[:,1]
rdist = 10.0*randn(len(star_data[:,0]))+20.0 # create a random dataset with the same number of values as the star data (easy version, more detail in lecture 6)
plt.xlabel("Temperature")
plt.ylabel("Magnitude")
plt.xlim(0,13000)
plt.ylim(-5,20)
plt.scatter(x,y,c=rdist,cmap='plasma_r') # the colourmap is unnecessary in this case, but just as an example
cbar = plt.colorbar()
cbar.set_label('Random')
plt.show()
It is very common when analyzing experimental data to include the errors in the measurement into the plot, as this gives more meaning to the values. A further common step is to guess a model that might fit the data and then to plot it to test the hypothesis. Let's use the star data again, and we can use a random dataset as the error in this case.
from numpy import loadtxt
from scipy import randn
import matplotlib.pyplot as plt
from sklearn import preprocessing # use scikit-learn to scale data
star_data = loadtxt("stars.txt",float)
scale_data = preprocessing.scale(star_data) # shifting the distribution of each attribute to have a mean of zero and a standard deviation of one (unit variance).
x = scale_data[:,0]
y = scale_data[:,1]
error = 0.2*randn(len(star_data[:,0]))
plt.figure(figsize=(7.5,5))
plt.rc("font", size=16)
plt.errorbar(x,y,abs(error),fmt="bo",capsize=5) # plot errorbars in blue with width 5 "caps"
plt.xlabel("Temperature")
plt.ylabel("Magnitude")
plt.show()
At first glance there appears to be two general trends - one seemingly linear and another almost a step function. We can guess some functional forms for these two trends and then plot them on the data to see how well they fit.
from numpy import loadtxt,linspace,exp
from scipy import randn
import matplotlib.pyplot as plt
from sklearn import preprocessing # use scikit-learn to scale data
# fitting functions
def linearF(x):
return a1*x + b1
def stepF(x):
return a2/(1+exp(x)) + b2 # try a sigmoid function
star_data = loadtxt("stars.txt",float)
# shifting the distribution of each attribute to have a mean of zero and a standard deviation of one (unit variance). Easier for guessing fits, and sometimes more efficient in machine learning.
scale_data = preprocessing.scale(star_data)
x = scale_data[:,0]
y = scale_data[:,1]
error = 0.2*randn(len(star_data[:,0]))
plt.figure(figsize=(7.5,5))
plt.rc("font", size=16)
plt.errorbar(x,y,abs(error),fmt="bo",capsize=5,zorder=1) # plot errobars in blue with witdth 5 "caps". Use zorder to make it at the "bottom".
plt.xlabel("Temperature")
plt.ylabel("Magnitude")
a1 = -0.3
b1 = 2.8
xmod1 = linspace(-2,6,100)
ymod1 = linearF(xmod1)
plt.plot(xmod1,ymod1,"r")
a2 = 5.0
b2 = -2.0
xmod2 = linspace(-4,6,100)
ymod2 = stepF(xmod2)
plt.plot(xmod2,ymod2,"y")
plt.show()
Obviously this is very inefficient way to actually fit data and there are a variety of powerful tools available for doing real fitting and calculating uncertainties.
In the previous example we introduced the colour scale for a third, random variable, but in many physical examples we actually will want to plot a two-dimensional function $F(x,y)$ e.g. temperature on a surface, particles incident on a detector. Here we make use of Python's capability for handling density plots. By default the package imshow
will automatically adjust the color-scale and axes for your input data - in this case the axes are just the rows and columns of the array. Note that the default is not how you would normally plot a graph, as the convention for matrices is that the row index is written first i.e. the vertical index rather than the horizontal.
from numpy import loadtxt
import matplotlib.pyplot as plt
circ_data = loadtxt("circular.txt",float)
print(circ_data) # note that it will not print all the data in this case, as it is pretty big, just the initial and final values.
plt.imshow(circ_data)
#plt.imshow(circ_data,origin="lower") # swap axes to that more conventional for plots
plt.show()
# rescale into region defined by [a,b,c,d] horizontal scale a-b and vertical scale c-d, but maintain an aspect ration of 2.0.
#plt.imshow(circ_data,origin="lower",extent=[0,10,0,5],aspect=2.0)
#plt.show()
# plt.xlim(0,10)
# plt.ylim(0,5)
# plt.imshow(circ_data,origin="lower",aspect=2.0) # showing that rescale is not the same as setting the axes limits
# plt.show()
The use of Pythonic methods to vectorize function calculations offers some very efficient methods for creating, handling and plotting density-like functions. Let's consider the waves generated from dropping a pebble into a pond, modelling them as the spreading out in the form of a sine
wave. If the centre is at $(x_1,y_1)$, then the distance $r_1$ from the centre to a point $(x,y)$ is:
$r_1 = \sqrt{(x-x_1)^2+(y-y_1)^2}$
where the height of the wave $h_1$ is defined by:
$h_1(x,y) = A_0 sin(kr_1)$
where $A_0$ is the amplitude of the waves and $k$ is the wavevector, $k=2\pi/\lambda$.
Then we can consider dropping a second pebble (exactly the same size and behaviour) at $(x_2,y_2)$ and assume that the waves will add linearly (a reasonable assumption for small waves), then the resultant height at $(x,y)$ will be:
$h(x,y) = A_0 sinkr_1 + A_0 sinkr_2$
Suppose the wavelength, $\lambda$, is 5 cm, the amplitude, $A_0$, is 1 cm and the centres of the impact circles are 20 cm apart, then we can create a programme to visualize the resulting interference pattern.
from numpy import pi,sqrt,sin
from numpy import linspace,meshgrid
import matplotlib.pyplot as plt
wavelength = 5.0
k = 2*pi/wavelength
amplitude = 1.0
separation = 20.0
side = 100.0 # side of the square we studying
points = 500 # number of grid points along each side
# calculate the positions of the centres of the circles
x1 = side/2 + separation
y1 = side/2
x2 = side/2 - separation
y2 = side/2
# calculate the height array based on the sum of the two sine waves
xvals = linspace(0,side,points) # create a float array with values from 0 - "side" with "points" intervals
yvals = xvals
x,y = meshgrid(xvals,yvals) # 2D arrays
r1 = sqrt((x-x1)**2+(y-y1)**2)
r2 = sqrt((x-x2)**2+(y-y2)**2)
h = amplitude*sin(k*r1) + amplitude*sin(k*r2)
# build the plot
plt.figure(figsize=(9,6))
plt.imshow(h,origin="lower",extent=[0,side,0,side])
cbar = plt.colorbar()
cbar.set_label("Amplitude")
plt.show()
Example - try playing around with the code for this and explore where the simple model breaks down.
Although density plots are a common approach for handling functions of two variables, another approach is to make a surface plot in 3D, with the height of the surface giving the value of the function at every position. This can give more insight for certain data or just be more visually appealing in a presentation. Let's try this by plotting the function $F(x,y) = cos(5x^2+y^2)$.
from numpy import cos
from numpy import linspace,meshgrid
import matplotlib.pyplot as plt
from matplotlib import cm # colourmaps for plotting
# generate the data
xyvals = linspace(-2,2,200)
x,y = meshgrid(xyvals,xyvals)
z = cos(5.0*x**2 + y**2)
# plot the function
plt.figure(figsize=(12,8))
plt.rc("font",size=12)
fig = plt.figure(figsize=(12,8))
ax = fig.add_subplot(projection="3d")
ax.plot_surface(x,y,z,cmap=cm.plasma)
# Note that this is an example of object-orientated programming (OOP). The general form is object.method(). Here ax is an instance of the class Axes3D, which we imported at the top of the program
# - by convention class names start with a capital letter. We created the object ax with the statement "ax = plt.gca(projection='3d')", which means "get current axes" (or create them if neccessary) of type "3d".
# If you look in the documentation for Axes3D, you will see that plot_surface() is one of its methods - things that Axes3D objects are able to do. You won't find set_xlabel() here, because that is a method that
# Axes3D inherited from it's base class Axes (which only knows how to do 2D plotting).
ax.set_zlim(-4.0,4.0)
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_zlabel("z")
plt.title(r"Surface plot of $F(x,y) = cos(5x^2+y^2)$")
plt.show()
As a practical example, read in the data file stm.txt
and plot it in density and surface styles. The file contains a grid of values from scanning tunneling microscope measurements of the (111) surface of silicon. A scanning tunneling microscope (STM) is a device that measures the shape of a surface at the atomic level by tracking a sharp tip over the surface and measuring quantum tunneling current as a function of position. Compare different plotting approaches and see if you understand what you are seeing.
Often in presentations and teaching environments, an animation of your results can be highly informative. Let's compare a static and dynamic version of a propagating wave packet defined by:
$f(x,t) = \frac{A_0 cos(k(x-v_p t))\exp(-(x-v_g t)^2}{2w^2}$
where $A_0$ is the amplitude, $w$ is the half-width, $k$ is the wavenumber, $v_g$ is the speed of the wavepacket (the group velocity) and $v_p$ is the speed of the waves (phase velocity).
from numpy import cos,exp
from numpy import linspace
import matplotlib.pyplot as plt
from numpy import cos,exp
# Parameters for the wave packet:
A0 = 3. # height
w = 0.5 # width
vg = 0.7 # group velocity
vp = 1.0 # phase velocity
k = 15. # wavenumber
def f(x,t):
"""Wave packet function."""
cpart = cos(k*(x-vp*t))
epart = exp(-(x-vg*t)**2/(2*w**2))
return A0*cpart*epart
xmax = 15.0 # x range of plot will be 0 to xmax
ymax = 4.0 # y range of plot will be -ymax to +ymax
x = linspace(0, xmax, 500) # array of x values for plotting
# Make the plot:
plt.rc('font', size=16)
plt.figure(figsize=(8,5))
t = 2.5 # chosen time
y = f(x,t) # array of y values at chosen time
plt.plot(x,y)
plt.text(xmax-3,-ymax+0.5,'t = {:.3f}'.format(t))
plt.xlim(0,xmax)
plt.ylim(-ymax,ymax)
plt.title('wave packet')
plt.xlabel('x')
plt.ylabel('f(x,t)')
plt.show()
In this case, we can vary the chosen time to explore how the wavepacket behaves as a function of time, but that is not very practical and does not really capture the relative differences in $v_g$ and $v_p$. Animating the plot would be more useful.
import matplotlib.pyplot as plt
from numpy import cos,exp
from numpy import linspace
from matplotlib.animation import FuncAnimation
# Parameters for the wave packet:
A0 = 3. # height
w = 0.5 # width
vg = 0.7 # group velocity
vp = 1.0 # phase velocity
k = 15. # wavenumber
def f(x,t):
"""Wave packet function."""
cpart = cos(k*(x-vp*t))
epart = exp(-(x-vg*t)**2/(2*w**2))
return A0*cpart*epart
xmax = 15.0 # x range of plot will be 0 to xmax
ymax = 4.0 # y range of plot will be -ymax to +ymax
x = linspace(0, xmax, 500) # array of x values for plotting
# Make the plot, but with no plotting data or label text:
plt.rc('font', size=16)
fig = plt.figure(figsize=(8,5))
curve = plt.plot([],[])[0] # this is the curves that we update in each frame
label = plt.text(xmax-3,-ymax+0.5,'') # the label on the legend also changes in each frame
plt.xlim(0,xmax)
plt.ylim(-ymax,ymax)
plt.title('wave packet')
plt.xlabel('x')
plt.ylabel('f(x,t)')
tmax = xmax/vg # pick so wave packet will reach x=xmax when t=tmax
nf = 200 # pick number of frames in animation
def frame0():
""""Initialization function,
used as background for all frames.
"""
curve.set_data([], [])
label.set_text('')
return curve,label
def frame(i):
""""Animation function,
called for each frame with i = frame number.
"""
t = tmax*i/nf
y = f(x,t)
curve.set_data(x, y)
label.set_text('t = {:.3f}'.format(t))
return curve,label
# Call the animator:
anim = FuncAnimation(fig, frame, init_func=frame0, frames=nf, interval=30, blit=True)
#anim.save('wavepacket.mp4') # optionally save the animation to a video file
plt.close() # supress automatic display of non-animated plot
plt.rc('animation', html='html5') # needed for Jupyter Notebook to display animation
anim # display the animation
Python also offers tools for creating simple 3D pictures and animations that can prove very useful for quick visualizations. In general, these cannot compete with dedicated software packages for features, but may be much more efficient for a lot of data. Most of these can be accessed by installing specific packages beyond MatPlotlib - a powerful example is Mayavi, which can be added to an Anaconda installation with conda install mayavi
. For those with skills in GPU programming, there is also Vispy. Here we introduce VPython which is a simpler, but still effective tool for 3D graphics (installation is more complex for Jupyter lab, but relatively painless for command line Python).
from vpython import sphere,vector,canvas
scene = canvas() # needed for Jupyter lab, but not in general
sphere()
There are a lot options for the creation of your 3D objects and they can be defined as variables that can later be manipulated e.g.
from vpython import sphere,vector,box,canvas
scene = canvas()
ball = sphere(radius=0.5,pos=vector(-5,0,0),color=color.blue)
wallR = box(pos=vector(6,0,0), size=vector(0.2,12,12), color=color.red)
This then opens the door to animations based on simple (or complex) physical ideas.
from vpython import sphere,vector,box,rate,canvas,color
scene = canvas()
ball = sphere(radius=0.5,pos=vector(-5,0,0),color=color.blue,make_trail=True)
wallR = box(pos=vector(6,0,0), size=vector(0.2,12,12), color=color.red)
ball.velocity = vector(25,0,0)
deltat = 0.005 # timestep of the ball position updates
t=0 # initial time
scene.autoscale = False # stop the window rescaling size
while t < 1: # run for 1 second
rate(100) # control the rate so that we see something
ball.pos = ball.pos + ball.velocity*deltat # physics!
t = t + deltat
Clearly we have very little of the physics included in this problem, and we leave it as an exercise if you want to make the ball bounce off the wall.
You can also harness the easy array commands in Python to rapidly create more complex structures in 3D form. Obvious examples can be found in atomic structures, where most systems are formed from repeating simple unit cells.
from vpython import sphere,canvas,color
from numpy import empty
scene = canvas()
size = 10
atoms = empty([size,3],float)
for x in range (0, size):
for y in range (0, size):
for z in range (0, size):
atoms = sphere(pos=vector(x, y, z), radius=1)
Example - Here the atoms are just placed on a grid, but as an exercise, visualize the crystal lattice of NaCl. It is simple cubic, with alternating Na and Cl ions at lattice sites, and a lattice constant of 0.564 nm. What does this distance represent?
As we have gone through the examples, we have commented on the balance between speed and simplicity, but the real war in computational physics is between speed and accuracy...they rarely have the inverse linear correlation one would wish for.
One the oldest, but still relevant errors in calculations are numerical errors. Normally these are due to a misunderstanding of how your particular code, programming language or computer architecture handles numbers. For example, by default Python has a precision of 16 significant figures for float
variables (integer
variables have arbitrary precision). As can be seen by printing the value of $\pi$.
from numpy import pi
print(pi)
This is generally enough for most applications, but physics is exactly one of those areas where precision really matters...consider this example.
x = 1e16
y = 10000000000000001.1823426537568
print(x-y)
While the error in representing $y$ is very small, the error in the difference between $x$ and $y$ is actually very significant, almost 50%. These types of errors can propagate throughout a code until all your results become meaningless. You can use numpy
packages to increase the precision, but in general it is just better to standardize or normalize your data so that these problems can be avoided and you don't have to work with numbers at the limit of precision.
In general, many physical questions solved computationally do not give absolute answers and we need to make a decision on how accurate we need the result versus the time it takes to provide it e.g. how many terms in the series? how many atoms in the crystal? As an example let's consider a quantum simple harmonic oscillator. It has energy levels $E_n = \hbar\omega(n+\frac{1}{2})$, where $n=0,1,2,\cdots\infty$. The average energy of a simple harmonic oscillator at temperature T is:
$$\langle E \rangle = \frac{1}{Z}\sum_{n=0}^{\infty}E_n e^{-\beta E_n}$$where $\beta=1/(k_B T)$ with $k_B$ being the Boltzmann constant, and $Z=\sum_{n=0}^{\infty}e^{-\beta E_n}{}$ is the canonical partition function. Let's consider that we want to evaluate the average energy $\langle E \rangle$ when $k_B T$ is 100. Since we are dealing with exponentials in $n$, several of the terms decrease rapidly and we can get a reasonable estimate with just the first 1000 terms in the sum. Using natural units ($\hbar=\omega=1$), we can solve this with the following code.
Note that in this example our choice of natural units means that $\frac{\hbar\omega}{k_{B}T}\ll 1$ and we are at the classical limit - the gap between quantum states is so small there is effectively a continuum and $\langle E \rangle$ should equal $k_{B}T$.
%%timeit -r3 -n10 # Remember timing will run the whole thing several times, if you don't specify it will evaluate what it needs for accuracy. Here we ask for 3 repeats of 10 loops.
from numpy import exp
terms = 1000
beta = 1/100
S = 0.0
Z = 0.0
for n in range(terms):
E = n + 0.5 # energy levels in natural units
weight = exp(-beta*E) # evaluate the exp only once as it is more expensive than standard arithmetic operations
S += weight*E
Z += weight
print(S/Z)
That was very fast, but we actually have no benchmark for how good this answer is without knowing the true solution (which you normally will not in computational physics). A common approach would be just to increase the number of terms in the sum to see how the answer changes.
%%timeit -r3 -n10 # remember timing will run the whole thing several times
from numpy import exp
terms = 10000
beta = 1/100
S = 0.0
Z = 0.0
for n in range(terms):
E = n + 0.5 # energy levels in natural units
weight = exp(-beta*E) # evaluate the exp only once as it is more expensive than standard arithmetic operations
S += weight*E
Z += weight
print(S/Z)
%%timeit -r3 -n10 # remember timing will run the whole thing several times
from numpy import exp
terms = 100000
beta = 1/100
S = 0.0
Z = 0.0
for n in range(terms):
E = n + 0.5 # energy levels in natural units
weight = exp(-beta*E) # evaluate the exp only once as it is more expensive than standard arithmetic operations
S += weight*E
Z += weight
print(S/Z)
At 100 000 terms the answer has completely converged to the precision we are using. This is likely much too accurate, but in general you need to decide for yourself what accuracy you need. Also notice that the computational cost was more or less linear with the number of terms, as you would expect - this is an Order N problem, and is as good as it gets. You will normally have to pay much more to increase your accuracy.
Obviously this benchmarking was a rather manual way of approaching the problem and could be done fully within the code itself. Let's do so...write a code that solves the problem (maybe using a different $k_b T$) as a function of the number of terms and then plots the convergence of the average energy. Plot also the time taken as a function of the number of terms included.
Finally, re-write the code for the quantum limit where $\frac{\hbar\omega}{k_{B}T}\gg 1$. What is this average energy that you get?
Methods for speed up of large calculations generally work by compiling Python code that would otherwise interpreted (assuming you have already vectorized your code using numpy arrays so the time-consuming parts are not interpreted). Two packages that do this are Cython and Numba - both come with Anaconda Distribution. A different avenue for speeding up large calculations is to use parallel or multi processing such as cluster computing, cloud and grid computing, etc.
jk1dict = {
"name": "Hrádecký",
"position": "Goalkeeper",
"year": 1989
}
print(jk1dict)
# assign variables using the key
x = jk1dict["year"]
print(x)
# change values using the key
jk1dict["year"] = 2020
print(jk1dict)
# print the length of the dictionary
print(len(jk1dict))
# add elements using dict, and then make a dictionary of dictionaries
jk2dict = dict(name="Toivio", position="Defender", year=1988)
jk3dict = dict(name="Sparv", position="Midfielder", year=1987)
jkdict = dict(player1=jk1dict,player2=jk2dict,player3=jk3dict)
print(jkdict)
A set
is an unordered collection of values and, if the format fits your data, are generally faster to handle for searching :
jkset = {"Hrádecký", "Toivio", "Sparv", "Pukki"}
print(jkset)
Most of the commands used for dict
are the same for controlling sets.
A tupl
is an ordered sequence of values, possibly of different types. It is written by separating the values with commas, and optionally or as necessary to prevent ambiguity surrounding the tuple by ordinary round parentheses:
a = 3,5
print(a)
print(type(a))
Tuples are very handy for making code simple and easy to read. Consider the following examples:
# use a tuple to swap values without using a temporary variable
x = 3
y = 4
x,y = y,x
print('x =',x,' y =',y)
# a function uses a tuple to return muliple values
def foofunc (x,y):
"""Function returns its two arguments squared."""
return x**2,y**2
a,b = foofunc(5,6)
print('a =',a,' b =',b)
# a list comprehension uses tuples to return a list of pairs of numbers
lst = [(n,m) for n in range(5) for m in range(5) if n > m]
print(lst)
In this last example (unlike the preceding ones) the parentheses surrounding the tuples are necessary to avoid ambiguity. Just like lists and numpy arrays, tuples can have zero or one element. A zero-element tuple is written (), while a one-element tuple is written (x,) or x,.
Object-oriented programming (OOP) is based on several useful ideas:
mylist
and supply the data by writing mylist = [2,4,6]
. Then, we can carry out the list-type method of appending the new element 8 to this list-type object by writing mylist.append(8)
.object.method(...)
object.item
Some languages like C do not have OOP features. Other languages like Java are completely based on OOP. Still other languages like C++ and Python can be used both for OOP and for non-OOP.
Here we define the class Bug, meant to represent a firefly (lightning bug) that can make flashes of light. By convention class names are capitalized to distinguish them from variables and functions, which should be spelled in lower case.
The definition of new class starts with class Name:
where Name
is the name of the class, and includes all of the indented code below this line. Within the class definition, methods are created by ordinary function definitions, except that the first argument of every method definition is self
. This "extra" argument is not supplied when the method is called (see examples below).
import numpy as np
from numpy.random import randn
class Bug:
"""Class to represent a firefly, 1st version.
Parameters
----------
xy : (optional) None or float,float
If None position of bug will be set randomly, otherwise
position will be supplied coordinates.
"""
def __init__(self,xy=None):
if xy:
self.x,self.y = xy # set position to supplied x,y
else:
self.x = randn()
self.y = randn() # pick x and y from normal dist
Our class Bug has a single method called __init__
(the word init
surrounded by double underscores, sometimes called "dunder init".)
Most classes will have an__init__
method, which will be called automatically when a new object in the class is created. To make a new Bug object called b
, you write b = Bug(xy)
. The arguments you supply to Bug
(xy
in this example) will be passed along to the the __init__
method.
Notice that you do not supply the self argument − this is supplied automatically whenever you call a class method, and will refer to the object itself. Looking at the code above, if we supply a value for xy
, then the two attributes self.x
and self.y
will be set equal to the two elements of xy
. Thus, if our bug object is called b
, then this code will make two variables that belong to b
, called b.x
and b.y
.
Finally, in the code above if we don't supply the argument xy
it will have a default value of None
which will count as False
in the if
statement, and the coordinates will be set randomly using randn
.
Here we make three Bug objects, two with random positions and one with specified position, and print out the coordinates of each bug:
# Make some bugs and print out their positions:
b1 = Bug() # b1 will have random position
b2 = Bug() # b2 will have random position
b3 = Bug((2,3)) # b3 will be at position x=2, y=3
print('b1: x = {:.3f}, y = {:.3f}'.format(b1.x,b1.y))
print('b2: x = {:.3f}, y = {:.3f}'.format(b2.x,b2.y))
print('b3: x = {:.3f}, y = {:.3f}'.format(b3.x,b3.y))
This simple example already shows some key features of OOP. Each Bug
object shares methods with other members of the class − in this case only __init__
. But each Bug
object has its own variables x
and y
. We don't need to give each Bug
object a separate name. Here we make a list of 1,000 bugs, print out the position of the 53rd bug, and make a plot showing the positions of all 1,000 bugs:
import matplotlib.pyplot as plt
# Make a list of 1,000 bugs at random positions:
bugs = [Bug() for b in range(1000)]
fav = 52 # starts at 0
# Print position of 53rd bug in the list:
print('b[{:d}]: x = {:.3f}, y = {:.3f}'.format(fav,bugs[fav].x,bugs[fav].y))
# Plot positions of all 1,000 bugs:
x = [b.x for b in bugs]
y = [b.y for b in bugs]
plt.rc('font',size=12)
plt.figure(figsize=(7,6))
plt.plot(x,y,'ro')
plt.title('1,000 bugs')
plt.xlabel('x')
plt.ylabel('y')
plt.show()
Here we put $N^2$ bugs at specified positions on an N×N
grid:
# Make N*N bugs on a regular grid:
nn = 20
bugs = [Bug((x,y)) for x in range(nn) for y in range(nn)]
# Plot their positions:
x = [b.x for b in bugs]
y = [b.y for b in bugs]
plt.figure(figsize=(7,6))
plt.plot(x,y,'ro')
plt.title('{} bugs'.format(nn*nn))
plt.xlabel('x')
plt.ylabel('y')
plt.show()
So far our bugs don't flash or do anything other than exist at location x,y
. To make them do something, we need to add more methods to the class beyond __init__
. To illustrate the syntax, here we add a (useless) method boo(n)
, which will print out the coordinates of the bug and print "Boo!" n
times:
class Bug:
"""Class to represent a firefly, 2nd version - include boo.
Parameters
----------
xy : (optional) None or float,float
If None position of bug will be set randomly, otherwise
position will be supplied coordinates.
"""
def __init__(self,xy=None):
if xy:
self.x,self.y = xy # set position to supplied x,y
else:
self.x = randn()
self.y = randn() # pick x and y from normal dist
def boo(self,n):
"""Prints coordinates of bug, then prints 'Boo' n times."""
print('Bug at x = {:.3}, y = {:.3} says'.format(self.x,self.y))
for j in range(n):
print(' Boo!')
b1 = Bug()
b2 = Bug() # make two bugs
b1.boo(3)
b1.boo(2)
b2.boo(5)
Things to note in this example:
boo
for the bug called b1
we write b1.boo(n)
(similar to the familiar lst.append('a')
which appends a
to the list lst
).boo
, we supply one less argument than in its definition − we do not supply self
.boo
, we can use self
to refer to "this object". Thus self.x
means "the x
variable for this bug".We want our bugs to flash on and off as time proceeds. To do this, we will add a boolean variable on
, which is True
when the bug is lit up, and False
when the bug is dark. We will also add a method step, taking no arguments, which will advance time by one unit for the bug.
Internally (in its little bug-brain), the bug must have some sort of clock that counts off the time steps. It is a Python convention that object variables that are intended for internal use only have names that start with an underscore, so we will call the clock variable _clk
. In addition to the clock, each bug needs to know how many time steps it will be lit up, and how many time steps are in a complete off+on period. Let's call these variables _osteps
and _per
, following the same underscore convention.
Finally, so each bug is a little different we will use the random-integer function randint
to chose the period and initial clock setting when the bug is initialized. Here then is a class definition for bugs that can flash:
from numpy.random import randn,randint
class Bug:
"""Class to represent a firefly, 3rd version - able to flash.
Parameters
----------
xy : (optional) None or float,float
If None position of bug will be set randomly, otherwise
position will be supplied coordinates.
"""
def __init__(self,xy=None):
if xy:
self.x,self.y = xy # set position to supplied x,y
else:
self.x = randn()
self.y = randn() # pick x and y from normal dist
self._osteps = 10 # flash lasts this many time steps
self._per = self._osteps + 50 + randint(20) # total steps in on+off
self._clk = randint(self._per) # start at a random place in cycle
self.on = self._clk<self._osteps # flash is last _osteps of the cycle
def step(self):
"""Steps time by one unit for this bug."""
self._clk -= 1 # count down clock by one step
if self._clk<0:
self._clk = self._per-1 # reset to start of cycle
self.on = self._clk<self._osteps # flash is last _osteps of the cycle
Note that the clock variable _clk
counts down and when it passes zero it is reset back to this bug's period (minus one), _per
- 1. The bug is "lit up" (on = True
) for the last _osteps
of this cycle. Here we make two bugs, and plot what they do over 500 time steps:
b1 = Bug()
b2 = Bug() # make two bugs
# Make arrays to plot on-off state of each bug:
steps = 500
b1on = np.zeros(steps) # for bug 1 start at zero
b2on = np.ones(steps) # for bug 2 start at one
for s in range(steps): # for each time step
# If a bug is on, bump up its array value by 0.8:
if b1.on:
b1on[s] += 0.8
if b2.on:
b2on[s] += 0.8
# Step time for both bugs:
b1.step()
b2.step()
# Plot the results:
plt.figure(figsize=(8,3))
plt.plot(b1on)
plt.plot(b2on)
plt.xlabel('time step')
plt.show()
We see that the two bugs start at different points in their on-off cycles, and that the relative time that they flash gradually changes because they have different periods. If we want to try this with a larger number of bugs, we can put the Bug object
in a list as was done above. Rather than making a separate list of plotting arrays for each bug, we can simply add a plotting array to each bug object. This illustrates a flexible aspect of OOP in Python - new attributes like this array can be added to an object at any time. In a less-freewheeling OOP language like Java, users of an object cannot even see the internal data structures of an object much less change them.
steps = 500 # number of time steps
bugs = [] # make a list of bugs
for j in range(15):
bugs.append(Bug())
bugs[j].y = j*np.ones(steps) # add a y-array offset by bug number j
# Track what they do over lots of steps:
for s in range(steps): # for each time step
for b in bugs: # for each bug
if b.on:
b.y[s] += 0.8 # if this bug on, bump up its line by 0.8
b.step() # step the bug
# Plot the results:
plt.figure(figsize=(8,5))
for b in bugs:
plt.plot(b.y)
plt.xlabel('time step')
plt.ylabel('bug number')
plt.show()
You can make a movie showing each flash at the x,y location of the bug that made it:
from matplotlib.animation import FuncAnimation
plt.rc('animation', html='html5') # needs this to write movies
bugs = [Bug() for b in range(200)] # make a list of bugs
def start_plot(xmin,xmax,ymin,ymax):
"""Make the plot, but with no plotting data or label text."""
fig = plt.figure(figsize=(7,7))
ax=plt.gca() # get current axes to change appearance
ax.set_facecolor('black') # set black background
ax.set_xticks([])
ax.set_yticks([]) # get rid of axis ticks and numbers
pts = plt.plot([],[],'yo')[0] # plot lit-up bugs as yellow dots
plt.xlim(xmin,xmax)
plt.ylim(ymin,ymax)
return fig,pts
def init():
"""Initialization function, used as background
for all frames."""
pts.set_data([], []) # background has no lit-up bugs
return pts,
def animate(i):
"""Animation function, called for each frame
with i = frame number."""
xon,yon = [],[] # lists of x,y of bugs that are on
for b in bugs:
if b.on:
xon.append(b.x)
yon.append(b.y)
b.step() # step the bug
pts.set_data(xon, yon) # plot points for the bugs that are on
return pts,
steps = 500 # set number of steps = animation frames
fig,pts = start_plot(-4,4,-4,4)
# Make and display the animation:
anim = FuncAnimation(fig,animate,init_func=init,\
frames=steps,interval=20,blit=True)
plt.close()
anim
This is supposed to show a tree full of fireflies. As there are 200 bugs, the flashes occur repeatedly at 200 specific locations in the plane.
Often we need to create a new class that is similar to an existing class but which has some additional features. Class inheritance makes it possible to build the new class starting with the existing class and only writing new code for the new things.
Here's an example: In Nature, groups of fireflies sometimes synchronize their flashing, by adjusting their flashing in response to the flashes produced by the other fireflies in the group. Bug
objects cannot synchronize since they have no information about what other Bug
objects are doing.
We can create a new class Bug2
in which each bug "sees" the flashes produced by other bugs. Here is one of many possible ways to do this: the step method gets a new argument, the total number bugs that are currently lit up. If the bug being time-stepped is not currently lit up, its clock is shifted by this number times a sensitivity parameter (set when the bug was created). The shift direction is chosen to speed the bug up if it is more than halfway through the cycle that ends in a flash, otherwise to slow the bug down. This will gradually shift the clock of the bug until it is lit up at the same time as most of the other bugs that it sees. Here is code that defines class Bug2
as a derived class with base class Bug
:
# code to define class Bug2
class Bug2(Bug):
"""Class to represent a firefly that can synchronize with
other fireflies.
Parameters
----------
xy : (optional) None or float,float
If None position of bug will be set randomly, otherwise
position will be supplied coordinates.
sens : (optional) float
Sensitivity of clock shift to on/off state of other bugs.
"""
def __init__(self,xy=None,sens=0):
Bug.__init__(self,xy) # initialize as for Bug class
self.sens = sens # save new sensitivity param
def step(self,non):
"""Steps time by one unit for this bug.
Parameters
----------
non : int
Number of other bugs that are on.
"""
if not self.on: # adjust clock only if not currently on
delta = self.sens*non # shift is sens times number of bugs on
if self._clk < 0.5*self._per: # more than halfway thru cycle,
self._clk -= delta # so speed up
else:
self._clk += delta # less than halfway, so slow down
Bug.step(self) # proceed with step as for Bug class
The statement class Bug2(Bug):
makes the variables and methods of the parent class Bug
available for Bug2
objects, unless overridden by the Bug2
class definition.
__init__
for Bug2
has a new argument sens
specifying the sensitivity of a bug to the flashes produced by other bugs. Note the Bug2
__init__
calls the initialization method of the parent class with the statement Bug.__init__(self,xy)
. Similarly the new step method uses Bug.step(self)
to call the parent time-step method. In this way, all of the code written for Bug
is re-used for Bug2
. If the parent class Bug
had any methods that were not re-defined in Bug2
(unlike this example), those methods would automatically be available for Bug2
objects.
We can update our program to use Bug2
objects by counting the number of bugs that are on and providing this as an argument to step
. With a suitable choice of the sensitivity parameter, these new bugs do indeed synchronize:
# Make a list of interacting bugs:
bugs = [Bug2(sens=0.007) for b in range(200)]
def animate2(i):
"""Animation function, called for each frame
with i = frame number."""
xon,yon = [],[] # make lists of x's and y's of bugs that are on
for b in bugs:
if b.on:
xon.append(b.x)
yon.append(b.y)
pts.set_data(xon,yon) # plot points for the bugs that are on
for b in bugs:
b.step(len(xon)) # step the bugs, giving number that are on
return pts,
steps = 2000 # set number of steps = animation frames
fig,pts = start_plot(-4,4,-4,4)
# Make and display the animation:
anim = FuncAnimation(fig,animate2,init_func=init,
frames=steps,interval=20,blit=True)
plt.close()
anim
The Bug2
bugs see all the other bugs equally much, no matter how far away the other bugs are. Since the distance between bugs has no effect on anything, position within the group has no real meaning and there is no possibility of spatial structure in the group behavior. In the language of statistical physics, this is called a mean-field model. A more realistic model with richer behavior has each bug "seeing" only the other bugs that are less than some fixed distance away. This is called a local model, since interactions between bugs have a limited spatial range.
To implement this local model we create a class Bug3
derived from Bug2
, with two new parameters giving the range of interactions and the list of other bugs already created. Each new Bug3
object builds a list of existing bugs that are within seeing range, and also updates the lists of these nearby bugs to include it.
To check that this is working correctly, we make 1000 bugs coloring a single chosen bug yellow, its in-range neighbours red, and all the other bugs blue. We expect to see a cloud of blue dots, with a single yellow dot somewhere in the cloud surrounded by a circular region of red dots.
class Bug3(Bug2):
"""Class to represent a firefly that can synchronize with
other nearby fireflies.
Parameters
----------
xy : (optional) None or float,float
If None position of bug will be set randomly, otherwise
position will be supplied coordinates.
sens : (optional) float
Sensitivity of clock shift to on/off state of other bugs.
rg : (optional) float
Distance range over with this bug can see other bugs.
buglist : (optional) list of Bug3
List of bugs already created.
"""
def __init__(self,xy=None,sens=0.0,rg=0.0,buglist=[]):
Bug2.__init__(self,xy,sens) # initialize as for Bug2 class
self.nearbugs = [] # construct list of other bugs that are in range
for b in buglist: # for other bugs created so far
d = ((self.x-b.x)**2 + (self.y-b.y)**2)**0.5
if d<rg: # if the other bug is in range,
self.nearbugs.append(b) # add the other bug to this bug's near-list,
b.nearbugs.append(self) # and add this bug to the other bug's near-list
def step(self):
"""Steps time by one unit for this bug."""
# Count number of nearby bugs that are on
nearon = [b.on for b in self.nearbugs].count(True)
Bug2.step(self,nearon) # use Bug2 step with nearon
bugs = [] # make 1000 limited-range interacting bugs
for j in range(1000):
bugs.append(Bug3(rg=1.0,buglist=bugs))
b1 = bugs[500] # pick an arbitrary bug as the central one
xnear,ynear = [],[]
xfar,yfar = [],[]
for b in bugs:
if b in b1.nearbugs:
xnear.append(b.x)
ynear.append(b.y)
else:
xfar.append(b.x)
yfar.append(b.y)
plt.figure(figsize=(6,6))
plt.plot(xnear,ynear,'ro')
plt.plot(xfar,yfar,'bo')
plt.plot([b1.x],[b1.y],'yo')
plt.xlabel('x')
plt.ylabel('y')
plt.show()
Here we create Bug3
bugs in two trees. The range that the bugs can see is large enough for them to see the other bugs in the same tree, but not in the other tree. As a result, the bugs in each tree synchronize, but the two trees of bugs are not synchronized with each other:
# Make a lot of limited-range interacting bugs
# living in two different trees:
bugs = []
for b in range(200):
bugs.append(Bug3(xy=(randn()-10,randn()),\
sens=0.01,rg=3.0,buglist=bugs))
bugs.append(Bug3(xy=(randn()+10,randn()),\
sens=0.01,rg=3.0,buglist=bugs))
steps = 2000 # set number of steps = animation frames
fig,pts = start_plot(-14,14,-4,4)
# Make and display the animation (use same animate function as
# used with Bug since Bug3 step does not need non):
anim = FuncAnimation(fig,animate,init_func=init,\
frames=steps,interval=20,blit=True)
plt.close()
anim
Bug3
locally-interacting bugs are capable of more interesting behaviors when placed in different geometries. Here we put them in a ring-shaped region. With suitable choices of the sensitivity and range parameters, propagating waves of flashes can be seen:
# Make a lot of limited-range interacting bugs
# living in a ring-shaped region:
from numpy.random import rand
from numpy import pi,cos,sin
bugs = []
for b in range(500):
r = 4 + rand() # ring from r=4 to r=5
theta = 2*pi*rand()
x = r*cos(theta)
y = r*sin(theta)
bugs.append(Bug3(xy=(x,y),sens=0.25,rg=.75,buglist=bugs))
steps = 2000 # set number of steps = animation frames
fig,pts = start_plot(-6,6,-6,6)
# Make and display the animation
anim = FuncAnimation(fig,animate,init_func=init,
frames=steps,interval=20,blit=True)
plt.close()
anim
Rather than thinking of fireflies, we might consider this as a model of a network of excitable elements like neurons which synchronize their firing to produce brain waves, or the muscle cells in the heart, which spontaneously synchronize their contractions.
For more (and almost endless options) you can begin with the official Python tutorial.