PHY407-Lab03

PHY407-Lab03
1 Computational Background
Forward and Centred Differences The forward and centred approximations of a derivative are
forward difference: df
dx
 
 
 
 
xi
≈
f(xi+1) − f(xi)
xi+1 − xi
, (1)
centred difference: df
dx
 
 
 
 
xi
≈
f(xi+1) − f(xi−1)
xi+1 − xi−1
. (2)
Gaussian quadrature Section 5.6 of Newman provides a nice introduction to Gaussian quadrature. Example 5.2 on p.170 gives you the tools you need to get started. The files referenced are gaussint.py and
gaussxw.py. To use this code, you need to make sure that both files are in the same directory (so the import
will work), or to know how to fetch them in different directories.
The line x, w = gaussxw(N) returns the N sample points x[0], . . . , [N-1] and the N weights w[0], . . . ,
w[N-1]. These weights and sample points can be used to calculate any integral on the interval −1 < x < 1.
To translate this integral into one that approximates an integral on the interval a < x < b you need to
implement the change of variables formulas (5.61) and (5.62) of Newman, which are written in the code as
xp = 0.5*(b-a)*x + 0.5*(b+a)
wp = 0.5*(b-a)*w
The loop then sums things up into the summation variable s.
On p. 171, there’s a bit of code that lets you skip the change of variables by using gaussxwab.py, which
you can use if you wish.
In Section 5.6.3 there’s a discussion of errors in Gaussian quadrature, which are somewhat harder to
quantify than for the previous methods we’ve seen. Equation (5.66) suggests that by doubling N we can get
a pretty good error estimate:
ϵN = I2N − IN . (3)
Factorials There is a function numpy.math.factorial that calculates the factorial of an integer. Your
program might run out of memory if the numbers in expressions such as 2nn! get too large, so wrap them
in floats as float(2**n)*float(factorial(n)).
2 Physics Background
The relativistic particle on a spring The energy of a relativistic particle, which is conserved, is given
by
E =
mc2
p
1 − v
2/c2
+
1
2
kx2
, (4)
1
which can be rearranged to show that
v
2 = c
2
"
1 −

mc2
E − kx2/2
2
#
. (5)
Assume the particle started from rest, from an initial position x0. Then E = mc2 +
1
2
kx2
0 and after
rearranging terms we can write the following expression for the positive root:
v =
dx
dt = c
(
k

x
2
0 − x
2

2mc2 + k

x
2
0 − x
2


/2

2 [mc2 + k (x
2
0 − x
2) /2]2
)1/2
= g(x), (6)
where g(x) is a function of x.
Small-amplitude (“classical”) limit: Notice that for k

x
2
0 − x
2


/2 ≪ mc2 we find v ≈
p
k (x
2
0 − x
2), as
expected for an energy conserving linear harmonic oscillator.
Large-amplitude limit: Notice that for k

x
2
0 − x
2


/2 ≫ mc2
, v approaches c but remains less than c.
Given (6), the period for the oscillation is given by four times the time taken for the particle to travel
from x = x0 to x = 0. Using separation of variables (see Example 5.10 in the text for a somewhat similar
example),
T = 4 Z x0
0
dx′
g(x
′)
. (7)
In the small and large amplitude limits described above, we expect T to approach 2π
p
m/k and 4x0/c,
respectively. Because g(x) → 0 as x → x
−
0
, the integral diverges in this limit and a fairly large number of
points will be required for an accurate calculation.
Quantum uncertainty in the harmonic oscillator In units where all the constants are 1, the wavefunction of the nth energy level of the one-dimensional quantum harmonic oscillator —i.e., a spinless point
particle in a quadratic potential well— is given by
ψn(x) = 1
p
2
nn!
√
π
exp(−x
2
/2) Hn(x), (8)
for n = 0 . . . ∞, where Hn(x) is the nth Hermite polynomial. Hermite polynomials satisfy a relation somewhat similar to that for the Fibonacci numbers,
Hn+1(x) = 2xHn(x) − 2nHn−1(x). (9)
The first two Hermite polynomials are H0(x) = 1 and H1(x) = 2x.
The derivative of the Hermite polynomials satisfy
dHn(x)
dx = 2nHn−1(x) (10)
and as a result,
dψn(x)
dx =
1
p
2
nn!
√
π
exp(−x
2
/2) [−xHn(x) + 2nHn−1(x)] . (11)
The quantum uncertainty of a particle in the n
th level of a quantum harmonic oscillator can be quantified
by its root-mean-square position p
⟨x
2⟩, where


x
2
 
=
Z ∞
−∞
x
2
|ψn(x)|
2
dx. (12)
This is also a measure of twice the potential energy of the oscillator. A similar calculation tells us that in
these units the momentum uncertainty is


p
2
 
=
Z ∞
−∞
 
 
 
 
ψn(x)
dx
 
 
 
 
2
dx, (13)
Pa
which is a measure of twice the kinetic energy of the oscillator.
Summing potential and kinetic energy, the total energy of the oscillator is then
E =
1
2

x
2
 
+


p
2
 
. (14)
3 Questions
1. [20%] Numerical differentiation errors Read section 5.10.2 in the textbook. Demonstrate that
the optimum step size for forward difference differentiation schemes is indeed √
C ≈ 10−8 by using the
following example.
(a) Consider the function f(x) = e
−x
2
. Using a forward difference scheme, numerically find the
function’s derivative at x = 0.5 for a range of h’s from 10−16 → 100
increasing by a factor of 10
each step. (You should have 17 values for h). Compare the value of each numerical derivative
that you get to the analytic value: specifically, take the absolute value of the numerical error for
the derivative with each value of h.
Submit the pseudocode, code, and a table or list of the numerical derivative values
and their errors at each value of h.
(b) Plot the error as a function of step size on a log-log plot and show that the minimum is indeed
at ≈ 10−8
. Explain the shape of the curve by considering equation (5.91) in the text. When does
the truncation error dominate? When does the rounding error dominate?
Submit the plot and written answers.
(c) Now implement a central difference scheme for the derivative and calculate the error for the same
function. Plot the error in the central difference scheme on the same plot as your error for the
forward difference scheme from Q1b. Concisely explain the key features of your plot. Does the
central difference scheme always clearly beat the forward difference scheme in terms of accuracy?
Submit the pseudocode, code, plot with both curves on it, and written answers.
2. [40%] The period of a relativistic particle on a spring
Using Gaussian quadrature, we will numerically calculate the period of the spring with the period
given by eqn. (7), and see how it transitions from the classical to the relativistic case. The idea is to
calculate T multiple times for a range of initial positions x0, assuming a mass of m = 1 kg and a spring
constant k = 12 N/m.
(a) The first issue is accuracy of the solution.
As x0 gets smaller, the period should approach 2π
p
m/k, which is the “classical” value. Check
with the information given in the Physics background that x0 = 1 cm indeed corresponds to the
classical limit.
Numerically calculate the period for N = 8 and N = 16 for x0 = 1 cm using eqn. (7), and compare
this to the classical value.
Estimate the fractional error for these two N’s.
Submit the code and written answers.
(b) To get a better sense of what affects the calculation, plot the integrands 4/gk and the weighted
values 4wk/gk at the sampling points.
Describe how these quantities behave as the x0 limit (i.e. the upper limit) of integration is
approached. How do you think this behaviour might affect accuracy of the calculation?
Suggestion: for the integrands plot and for the weighted values plot, put the N = 8 and N = 16
cases on the same plot (using different colours), for better readability.
Submit plots and written answers.
Page
(c) For a classical particle on a spring
m
d
2x
dt = −kx, (15)
what initial displacement x0 = xc > 0 for a particle initially at rest would lead to a speed equal
to c at x = 0?
Submit written answer.
(d) For N = 200, what is your estimate of the percentage error for the small amplitude case?
Submit written answer.
(e) Plot T as a function of x0 for x0 in the range 1 m< x < 10xc, and compare it to the classical limit
and to the large-amplitude relativistic limit as suggested at the beginning of the problem.
Submit plot and written answers. Also submit pseudocode and code for the entire
question.
3. [40%] Calculating quantum mechanical observables (Adapted and expanded from Newman 5.13,
p. 182)
(a) Write a user-defined function H(n,x) that calculates Hn(x) for given x and any integer n ≥ 0.
Use your function to make a plot that shows the harmonic oscillator wavefunctions for n = 0, 1,
2, and 3, all on the same graph, in the range −4 ≤ x ≤ 4.
Submit figure and explanatory notes (code to be submitted later).
(b) Make a separate plot of the wavefunction for n = 30 from x = −10 to x = 10.
Note: the program should take only a second or so to run.
Submit figure.
(c) Write a program that
• evaluates 

x
2
 
,


p
2
 
and energy E, as described in the Physics Background, using Gaussian
quadrature on 100 points, and then
• calculates the position and momentum uncertainty (i.e., the root-mean-square position and
momentum of the particle) for a given value of n.
Use your program to calculate the uncertainty for n = (0, 1, 2, . . . , 15).
What is the relationship between the uncertainty in position and the uncertainty in momentum?
Do you notice a simple rule for the energy of the oscillator?
Note: you will need to use one of the evaluation techniques on p.179 of Newman to deal with the
improper integrals here. For n = 5,
p
⟨x
2⟩ ≈ 2.35.
Submit pseudocode, code, printed output, and written answers.
Page 4