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 4×0/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

## Reviews

There are no reviews yet.