Chapter 5: Integrals

Exercises are from "Computational Physics" by Mark Newman.
Programmed by Sasha Bakker.

In [1]:
from numpy import *
import matplotlib.pyplot as plt
import math
from scipy import integrate

Exercise 5.1: Trapezoidal rule

Download velocities.txt, which contains two columns of numbers, the first representing time t in seconds and the second the x-velocity in meters per second of a particle, measured once every second from time $t=0$ to $t=100$.

Our approximation for the area under the whole curve is the sum of the areas of the trapezoids for all $N$ slices:

$I(a,b) \simeq \Sigma_{k=1}^{N} A_{k} = h \left[\frac{1}{2}f(a)+\frac{1}{2}f(b) + \Sigma_{k=1}^{N-1} f(a+kh) \right]$
where the interval is divided from $a$ to $b$ into $N$ slices each with width $h$.

In [24]:
data = loadtxt("velocities.txt",float)
t = data[:,0]
v = data[:,1]
N = len(t)-1 # number of slices
a = t[0]
b = t[-1]

def trap(v,N):
    h = (b-a)/(N)
    
    weight = ones(N+1) # number of points is number of slices + 1
    weight[0] = 0.5
    weight[-1] = 0.5
    
    return h*sum(weight*v)

plt.figure(figsize=(6,3))
plt.title("Particle velocity (m/s) vs. time (s)")
plt.plot(t,v,color='purple')
plt.xlim(a,b)
plt.show()

area = trap(v,N)
print(f"The approximate distance travelled in the x-direction as a function of time is {round(area,2)} meters.")

The approximate distance travelled in the x-direction as a function of time is 8.22 meters.

Exercise 5.2: Simpson's rule

Approximate the integral $\int_{0}^{2} (x^{4}-2x +1) dx$ with Simpson's rule given by

$I(a,b) \simeq \frac{1}{3} h \left[f(a) + f(b) + 4 \Sigma_{k=1}^{N/2}f(a+(2k-1)h) + 2 \Sigma_{k=1}^{N/2-1}f(a+2kh)\right]$.
In [13]:
N = 999 # number of slices
a = 0
b = 2

def f(x):
    return x**4-2*x+1

def simpson(f,N):
    if N%2 !=0: 
        N += 1 # must have an even number of slices

    h = (b-a)/N
    x = linspace(a,b,N+1) # number of points is number of slices + 1

    weight = 4*ones(N+1) 
    weight[::2] = 2
    weight[0] = 1
    weight[-1] = 1

    return ((h/3)*sum(weight*f(x)),N)

def trapezoid(f,N):
    h = (b-a)/(N)
    x = linspace(a,b,N+1) # number of points is number of slices + 1
    
    weight = ones(N+1) # number of points is number of slices + 1
    weight[0] = 0.5
    weight[-1] = 0.5
    
    return h*sum(weight*f(x))

values = simpson(f,N)
area = values[0]
N = values[1]
print("Simpson's Rule:")
print(f"The approximate area under the curve is {area}.")
print(f"The functional error on the calculation is {abs(area-4.4)} for {N} slices.\n")

area = trapezoid(f,N)
print("Trapezoidal Rule:")
print(f"The approximate area under the curve is {area}.")
print(f"The functional error on the calculation is {abs(area-4.4)} for {N} slices.")

Simpson's Rule:
The approximate area under the curve is 4.400000000004266.
The functional error on the calculation is 4.2659209498197015e-12 for 1000 slices.

Trapezoidal Rule:
The approximate area under the curve is 4.400010666665601.
The functional error on the calculation is 1.0666665600567171e-05 for 1000 slices.

Exercise 5.3: More integral approximation

Consider the integral $E(x) = \int_{0}^{x} e^{-t^{2}} dt$. Note there is no known way to perform this particular integral analytically, so numerical approaches are the only way forward.

In [25]:
a = 0
N = 1000
x = arange(a,3,0.1)
area = list()

def f(x):
    return math.e**(-x**2)

for b in x:
    area.append(simpson(f,N)[0])

plt.figure(figsize=(6,3))
plt.title("E(x) vs. x")
plt.plot(x,area,color='purple')
plt.xlim(x[0],x[-1])
plt.show()

Exercise 5.4: The diffraction limit of a telescope

Our ability to resolve deatil in astronomical observations is limited by the diffraction of light in our telescopes. Light from stars can be treated effectively as coming from a point source at infinity. When such light, with wavelength $\lambda$ passes through the circular aperture of a telescope and is focused by the telescope in the focal plane, it produces not a single dot, but a circulat diffraction pattern consisting of central spot surrounded by a series of concentric rings. The intensity of the light in this diffraction pattern is given by

$I(r) = \left( \frac{J_{1}(kr)}{kr} \right)^2 $
where $r$ is the distance in the focal plane from the center of the diffration pattern, $k=2\pi / \lambda$, and $J_{1}(x)$ is a Bessel function. The Bessel functions $J_{m}(x)$ are given by

$ J_{m}(x) = \frac{1}{\pi} \int_{0}^{\pi} \cos (m\theta - x \sin \theta) d\theta$
where $m$ is a nonnegative integer and $x\geq 0$.

In [26]:
N = 1000
a = 0
b = pi
x1 = 0
x2 = 20
xvals = linspace(x1,x2,N+1)
M = [0,1,2]

def f(theta):
    return cos(m*theta-xval*sin(theta))

def J(m,xval,f,N):
    return (1/pi)*simpson(f,N)[0]

plt.figure(figsize=(6,3))
plt.title("Bessel functions")

for m in M:
    Jvals = list()
    for xval in xvals:
        Jvals.append(J(m,xval,f,N))
    plt.plot(xvals,Jvals,label='m='+str(m))

plt.xlim(x1,x2)
plt.legend()
plt.show()
In [21]:
a = 0
b = pi
N = 1000

lam = 500 # nm
k = (2*pi)/lam
l = 500
nn = 100
x = linspace(-l,l,nn)
y = x
x,y = meshgrid(x,y)
r = sqrt(x**2+y**2)
I = zeros([nn,nn])

j1 = lambda t: cos(t-(m)*sin(t))

for j in range(len(y)):
    for i in range(len(x)):
        
        rval = r[i,j]
        m = rval*k
        number = (1/pi)*integrate.quad(j1,0,pi)[0]
        number = (number/m)**2
        I[i,j] = number

print("Via scipy.integrate.quad()")
plt.figure(figsize=(6,6))
plt.title("Diffraction of a telescope")
plt.imshow(I,origin='lower',vmax=0.01,extent=(-N,N,-N,N))
plt.set_cmap("twilight")
plt.colorbar()
plt.show()

Via scipy.integrate.quad()

Exercise 5.5: Error on Simpson's Rule

For an integral over the interval from $a$ to $b$ the approximation error of Simpson's rule is given to leading order by

$\epsilon =\frac{1}{90}h^{4}\left[f^{'''}(a)-f^{'''}(b)\right]$.
This allows us to calculate the error on our integrals provided we have a known closed-form exression for the integrand f(x), so that we can calculate the derivatives that appear in the formulas. Suppose we are evaluating an integral over the interval from $x=a$ to $x=b$ using Simpson's rule. Let us perform the integral with some number of steps $N_{1}$, so that the step size is $h_{1}=(b-a)/N_{1}$.

Then here's the trick: we now double the number of steps and perform the integral again. That is we define the new number of steps $N_{2}=2N_{1}$ and a new step size $h_{2}=(b-a)/N_{2}=\frac{1}{2}h_{1}$ and we reevaluate the integral using Simpson's rule, giving a new answer, which will normally be more accurate than the previous one. As we have seen, Simpson's rule introduces an error of order $O(h^{4})$, which means when we half the value of $h$ we are dividing the size of our error by $16$. Knowing this fact allows us to estimate how big the error is.

Suppose that the true value of our integral is $I$ and let us denote our first estimate using Simpson's rule with $N_{1}$ steps by $I_{1}$. The difference between the true value and the estimate, which is the error on the estimate, is proportional to $h^{4}$, so let us wrtie it as $ch^{4}$, where $c$ is a constant. Then $I$ and $I_{1}$ are related by $I = I_{1}+ch_{1}^{4}$.

We can also write a similar formula for our second estimate $I_{2}$ of the integral, with $N_{2}$ steps: $I=I_{2}+ch_{2}^{4}$. Equating the expressions for $I$ we then get

$I_{1}+ch_{1}^{4}=I_{2}+ch_{2}^{4}$,
or
$I_{2}-I{1}=ch_{1}^{4}-ch_{2}^{4}=15ch_{2}^{4}$,
where we have made use of the fact that $h_{1}=2h_{2}$. Rearranging this expression then gives the error $\epsilon_{2}$ on the second estimate of the integral to be
$\epsilon_{2}=\frac{1}{15} (I_{2}-I_{1})$.

Exercise 5.6: Error on the trapezoidal rule

Calculate the value of the integral $\int_{0}^{2}(x^{4}-2x+1)dx$ with $N_{1}=10$ slices, and estimate the error of the result, given by

$\epsilon_{2}=\frac{1}{3} (I_{2}-I_{1})$.
In [17]:
N1 = 10 # number of slices
a = 0
b = 2

def f(x):
    return x**4-2*x+1

I1 = trapezoid(f,N1)
I2 = trapezoid(f,2*N1)
error = (1/3)*(I2-I1)
print(f"The estimated error on the integral for {2*N1} slices is {error}")
print(f"Directly calculating the error via difference between actual and calculated result produces {4.4-I2}")
print(f"The percent difference between estimated and true error is {round(100*abs(4.4-I2-error)/abs(4.4-I2),2)}%")

The estimated error on the integral for 20 slices is -0.026633333333333137
Directly calculating the error via difference between actual and calculated result produces -0.026660000000000572
The percent difference between estimated and true error is 0.1%

Exercise 5.7a: Adaptive trapezoidal rule

Consider the integral $I=\int_{0}^{1}\sin ^{2}\sqrt{100x}dx$.

In [3]:
a = 0 # integration lower bound
b = 1 # integration upper bound
N = 1 # intial number of integration slices
e = 1e-6 # accuracy bound

def f(x):
    return (sin(sqrt(100*x)))**2

def trapezoid_adapt(f,N,I0):
    h = (b-a)/(N)
    x = linspace(a,b,N+1) # number of points is number of slices + 1
    
    weight = ones(N+1) # number of points is number of slices + 1
    weight[0] = 0
    weight[-1] = 0
    weight[::2] = 0 # summing over odd only
    
    I = 0.5*I0 + h*sum(weight*f(x))
    
    return I, (1/3)*(I-I0)


Ivals, Evals, Nvals = list(), list(), list()
I0, E0 = trapezoid_adapt(f,N,0)[0], 1

while abs(E0) > e:
    
    N *= 2
    I0, E0 = trapezoid_adapt(f,N,I0)
    Ivals.append(I0)
    Evals.append(E0)
    Nvals.append(N)

hz = ones(len(Nvals))*1e-6
plt.figure(figsize=(7,4))
plt.title("Estimated error of the adaptive trapezoidal integration")
plt.ylabel("Estimtated error")
plt.xlabel("Number of slices")
plt.plot(Nvals,Evals,'b--')
plt.plot(Nvals,hz,'k')
plt.xlim(0,200)
plt.show()

print(f"An accuracy of 1e-6 is reached at the number of slices N = {N}")

An accuracy of 1e-6 is reached at the number of slices N = 65536

Exercise 5.7b: Romberg integration

In [4]:
def trapezoid(f,N):
    """
    Regular Trapezoidal Rule
    """
    h = (b-a)/(N)
    x = linspace(a,b,N+1) # number of points is number of slices + 1
    weight = ones(N+1) # number of points is number of slices + 1
    weight[0] = 0.5
    weight[-1] = 0.5
    
    return h*sum(weight*f(x))

I1 = trapezoid(f,1)
R0 = list()
R0.append(I1)
error = 1
h = (b-a)/(N)

def romberg(R0,m_max):
    
    R = list()
    error = 1.0
    
    for m in range(m_max):
        
        m = int(m)
        if m == 0: 
            val = trapezoid(f,m_max)
            error = (1/3)*(val-R0[0])
        else:
            num = 1/(4**m-1)*(R[m-1]-R0[m-1])
            val = R[m-1]+num
            error = num + h**(2*m+2)
            
        R.append(val)
    
    return R, error

i = 1
print(f"Slice {i}: {R0}\n")
while abs(error) > e:
    
    i += 1 
    R0, error = romberg(R0,i)
    print(f"Slice {i}: {R0}\n")

print(f"The integration achieves an accuracy of 1e-6 for the number of slices N = {i}")

Slice 1: [0.147979484546652]

Slice 2: [0.3252319078064746, 0.38431604889308213]

Slice 3: [0.43079757183944845, 0.4659861265171064, 0.47143079835870805]

Slice 4: [0.5122828507233315, 0.5394446103512924, 0.5443418426069049, 0.5454991607695747]

Slice 5: [0.45902066640404526, 0.4412666049642832, 0.43472140460514924, 0.4329813976527404, 0.43254015152287045]

Slice 6: [0.42216668286887116, 0.40988202169047977, 0.40778971613889287, 0.4073622290203809, 0.40726176169241085, 0.40723705163392554]

Slice 7: [0.4067588542747065, 0.4016229114099849, 0.4010723040579519, 0.40096567846936554, 0.4009405939574008, 0.40093441490780646, 0.40093287580237885]

Slice 8: [0.40299744847824825, 0.4017436465460955, 0.40175169555516954, 0.40176247954718886, 0.40176560425729796, 0.40176641071897917, 0.4017666138925594, 0.4017666647830001]

The integration achieves an accuracy of 1e-6 for the number of slices N = 8

Exercise 5.8: Adaptive Simpson's rule

In [5]:
N = 1 # intial number of integration slices
f_a = f(a)
f_b = f(b)


def simpsons_adapt(f,N,I0):
    h = (b-a)/(N)
    x = linspace(a,b,N+1) # number of points is number of slices + 1
    fx = f(x)
    
    weight = ones(N+1) # number of points is number of slices + 1
    weight[0] = 0
    weight[-1] = 0
    weight_S, weight_T = weight, weight
    weight_S[1::2] = 0 # summing over even only
    weight_T[0::2] = 0 # summing over odd only
    
    S = (1/3)*(f_a + f_b + 2*sum(weight_S*fx))
    T = (2/3)*sum(weight_T*fx)
    I = h*(S+2*T)
    return I, (1/15)*(I-I0) # integral, estimated error

Ivals, Evals, Nvals = list(), list(), list()
I0, E0 = simpsons_adapt(f,N,0)

while abs(E0) > e:
    
    N *= 2
    I0, E0 = simpsons_adapt(f,N,I0)
    Ivals.append(I0)
    Evals.append(E0)
    Nvals.append(N)

hz = ones(len(Nvals))*1e-6
plt.figure(figsize=(7,4))
plt.title("Estimated error of the adaptive Simpson's integration")
plt.ylabel("Estimtated error")
plt.xlabel("Number of slices")
plt.plot(Nvals,Evals,'b--')
plt.plot(Nvals,hz,'k')
plt.xlim(0,200)
plt.show()

print(f"An accuracy of 1e-6 is reached at the number of slices N = {N}")

An accuracy of 1e-6 is reached at the number of slices N = 8192

Gaussian quadrature

The python file guassxw.py written by Mark Newman can be found here. We define the Guassion quadrature method of integration, which will be some exercises, in the following cell.

In [2]:
from gaussxw import gaussxwab
def gquad(function, b):

    xp, wp = gaussxwab(N,a,b)

    s = 0.0

    for k in range(N):
        s += wp[k]*function(xp[k])

    return s

Exercise 5.9: Heat capacity of a solid

Debye's theory of solids gives the heat capacity of a solid at temperature $T$ to be

$C_V= 9V\rho k_{B} \int_{0}^{\theta_{D} /T} \frac{x^{4}e^{x}}{(e^{x}-1)^{2}}dx$

where $V$ is the volume of the solid, $\rho$ is the number density of atoms, $k_{B}$ is the Boltzmann's constant, and $\theta_{D}$ is the so-called Debye temperature, a property of solids that depends on their density and speed of sound.

In [27]:
temps = linspace(5,500,500)

V = 1e-6 # cubic meters
rho = 6.022e28 # m**-3
debye = 428 # K
kB = 1.38064852e-23 # m**2 kg s**-2 K**-1
N = 50 # sample points
a = 0


def dcdt(x):
    return (9*V*rho*kB*(temp/debye)**3)*(x**4*exp(x))/(exp(x)-1)**2


S = list()

for temp in temps:
    b = debye/temp
    S.append(gquad(dcdt,b))

plt.figure(figsize=(7,4))
plt.plot(temps,S,color='purple')
plt.title("Heat capacity of solid aluminum")
plt.xlabel(r"Temperature ($K$)")
plt.ylabel(r"$C_V$ ( $J \cdot K^{-1}$)")
plt.xlim(temps[0],temps[-1])
plt.show()

Exercise 5.10: Period of an anharmonic oscillator

The harmonic oscillator corresponds to a quadratic potential $V(x) \propto x^{2}$. Any other form gives an anharmoic oscillator. One way to calculate the motion of an oscillator is to write down the equation for the conservation of energy in the system. If the particle has mass $m$ and position $x$, then the total energy is equal to the sum of the kinetic and potential energies thus:

$E = \frac{1}{2}\left(\frac{\text{dx}}{\text{dt}}\right)^2 + V(x)$

Since the energy must be constant over time, this equation is effectively a (nonlinear) differential equation linking $x$ and $t$.


Lwr us assume that the potential $V(x)$ is symmetric about $x=0$ and let us set our anharmonic oscillator going with amplitude $a$. That is, at $t=0$ we release it from rest at position $x=a$ and it swings back towards the origin. Then at $t=0$ we have $\text{dx}/ \text{dt}=0$ and the equation above reads $E=V(a)$, which gives us the total energy of the particle in terms of the amplitude.


When the particle reaches the origin for the first time, it has gone through one quarter of a period of the oscillator. By rearranging the equation above for $\text{dx}/ \text{dt}$ and then integrating with respect to $t$ from $0$ to $\frac{1}{4}T$, we show that the period $T$ is given by:

$T = \sqrt{8m} \int_{0}^{a} \frac{dx}{V(a)-V(x)}$

Let us suppose the potential is $V(x) = x^{4}$ and the mass of the particle is $m=1$.

In [10]:
from gaussxw import gaussxwab
m = 1 # kg
N = 20 # points
a = 0
amplitudes = linspace(0.1,2,N)


def dTdt(x):
    
    return sqrt(8*m)/sqrt(amp**4-x**4)
    
S = list()
for amp in amplitudes:
    S.append(gquad(dTdt,amp))

plt.figure(figsize=(7,4))
plt.title("Period of an anharmonic oscillator")
plt.xlabel("Amplitude (m)")
plt.ylabel("Period (s)")
plt.plot(amplitudes,S,'k.')
plt.show()

Notice that the oscillator gets faster as the amplitude increases, even though the particle has further to travel for larger amplitude. This is due to the stronger restoring force the larger amplitude supplies.

Exercise 5.11: Plane wave diffraction

Suppose a plane wave of wavelength $\lambda$, such as light or a sound wave, is blocked by an object with a straight edge. The wave will be diffracted at the edge and the resulting intensity at the position $(x,z)$ is given by near-field diffraction theory to be:

$I=\frac{I_{0}}{8}\left([2C(u)+1]^2+[2S(u)+1]^2\right)$

where $I_{0}$ is the intensity of the wave before diffraction and

$u = x\sqrt{\frac{2}{\lambda z}}$, $\hspace{1cm} C(u)=\int_{0}^{u}\cos{\frac{1}{2}\pi t^2}dt$ $\hspace{1cm}S(u)= \int_{0}^{u} \sin{\frac{1}{2}\pi t^2}dt$

In [11]:
a = 0
lam = 1 # m
z = 3 # m
N = 50 # points

def u(x):
    return x*sqrt(2/(lam*z))

def dCdt(t):
    return cos(0.5*pi*t**2)

def dSdt(t):
    return sin(0.5*pi*t**2)

def ratio(x):
    val_u = u(x)
    val_S = gquad(dSdt,val_u)
    val_C = gquad(dCdt,val_u)
    return (1/8)*(2*val_C+1)**2 + (1/8)*(2*val_S+1)**2
    
x_vals = linspace(-5,5,N) # (-5,5) m
r_vals = list()
for x in x_vals:
    r_vals.append(ratio(x))

plt.figure(figsize=(7,4))
plt.xlabel(r"$x$ (m)")
plt.ylabel(r"$\frac{I}{I_{0}}$")
plt.plot(x_vals,r_vals,'k--')
plt.title(r"$I \div I_{0}$ as a function of $x$")
plt.show()

Exercise 5.12: The Stefan-Boltzmann constant

The Planck theory of thermal radiation tells us that in the (angular) frequency interval $\omega$ to $\omega + d\omega$, a black body of unit area radiates electromagnetically an amount of thermal energy per second equal to $I(\omega) d\omega$, where

$I(\omega) = \frac{\hbar}{4 \pi^2 c^2} \frac{\omega^3}{\text{e}^{\hbar\omega / k_{B} T}-1}$.

Here $\hbar$ is Planck's constant over $2\pi$, $c$ is the speed of light, and $k_{B}$ is Boltzmann's constant. From this we derive the total rate at which energy is radiated by a black body per unit area over all frequencies:

$W = \frac{k_{B}^{4}T^{4}}{4 \pi^{2}c^{2}\hbar^{3}} \int_{0}^{\infty} \frac{x^{3}}{e^{x}-1}dx$.
In [62]:
"""
## Warning filter
import sys
if sys.path[0] == '':
    del sys.path[0]
import warnings
warnings.filterwarnings("ignore", category=RuntimeWarning) 
"""

kB = 8.6173e-5 # eV / K
hbar = 4.1357e-15/(2*pi) # ev*s 
c = 3e8 # m/s
N=50 # steps
T = 293.15 # room temp
a = 0

def dWdt(z):
    cst = (kB**4*T**4)/(4*pi**2*c**2*hbar**3)
    return cst*(z**3) / ((1-z)**5*(exp(z/(1-z))-1))

W1 = gquad(dWdt,1)
N = 2*N
W2 = gquad(dWdt,1)
error = (W2-W1)


stefan = W1/T**4
stefan_error = (W1+error)/T**4-stefan
stefan, stefan_error = stefan*1.60218e-19, stefan_error*1.60218e-19

print(f"The Stefan-Boltzmann constant is computed with N = {int(N/2)} slices at RT to be:\n {stefan} J / (m^2 s K^4)")
print(f"\nThe error due to the number of slices for the Gaussian quadrature "
       f"integration is:\n{stefan_error} J / (m^2 s K^4)")

The Stefan-Boltzmann constant is computed with N = 50 slices at RT to be:
 5.6623233113355344e-08 J / (m^2 s K^4)

The error due to the number of slices for the Gaussian quadrature integration is:
1.5275960958251952e-17 J / (m^2 s K^4)

The accpeted value for the Stefan-Boltzmann constant to three significant digits is $5.67\cdot10^{-8} \frac{\text{J}}{\text{m}^2\cdot\text{s}\cdot\text{K}^{4}}$. Our computation, to three significant digits is $5.66\cdot10^{-8} \frac{\text{J}}{\text{m}^2\cdot\text{s}\cdot\text{K}^{4}}$. Because the error on the calculation is negligible, the error is due to the method of integration chosen.

Exercise 5.13: Quantum uncertainty in the harmonic oscillator

In units where all the constants are $1$, the wavefunction of the $n$th energy level of the one-dimensional quantum oscillator is given by

$\psi_{n}(x)=\frac{1}{\sqrt{2^{n}n!\sqrt{\pi}}}e^{-x^{2}/2H_{n}(x)}$.

for $n=0...\inf$, where $H_{n}(x)$ is the $n$th Hermite polynomial. Hermite polynomials satisfy a relation somewhat similar to that for the Fibonacci numbers, although more complex:

$H_{n+1}(x)=2xH_{n}(x)-2nH_{n-1}(x)$
In [29]:
def H(n,x):
    if n < 0:
        print("Error: n cannot be less than 0")
    elif n == 0:
        return 1
    elif n == 1:
        return 2*x
    else:
        return 2*x*H(n-1,x)-2*(n-1)*H(n-2,x)
    
N = 500
xs = linspace(-4,4,N)
ns = [0,1,2,3]
Hvals = zeros((4,N))


for n in ns:
    Hvals[n,:] = H(n,xs)
   
plt.figure(figsize=(7,4))

for n in ns:
    plt.plot(xs,Hvals[n,:],label="n = "+str(n))

plt.title("Hermite polynomials vs. x")
plt.ylim(-30,30)
plt.xlim(-4,4)
plt.legend()
plt.show()
In [30]:
from math import factorial

def psi(n,x):
    
    n = int(n)
    coeff = 1/sqrt( 2**n*factorial(n)*sqrt(pi) )
    ex = exp(-x**2/2)
    herm = H(n,x)
    
    return coeff*ex*herm

Pvals = zeros(N)
xs = linspace(-10,10,N)
plt.figure(figsize=(7,4))
plt.title("Wavefunction vs. x")
n = 30
Pvals = psi(n,xs)
plt.plot(xs,Pvals,label="n = "+str(n),color='purple')
plt.xlim(-10,10)
plt.legend()
plt.show()

The quantum uncertainty in the position of a particle in the $n$th level of a harmonic oscillator can be quantified by its root-mean-square position $\sqrt{\langle x^{2} \rangle}$, where

$\langle x^{2} \rangle = \int_{-\infty}^{\infty} x^{2} |\psi_{n}(x)|^{2}dx$
In [160]:
N = 100 # points
n = 5
a = -1

def dfdz(z):
    coeff = ((z/(1-z**2))**2)*(1+z**2)/((1-z**2)**2)
    wf = psi(n,(z/(1-z**2)))
    return coeff*abs(wf)**2

val = gquad(dfdz,1)
print(f"The uncertainty for n = {n} is {round(sqrt(val),3)}")

The uncertainty for n = 5 is 2.345

Exercise 5.14: Gravitational pull of a uniform sheet

A uniform square sheet of metal is floating motionless in space (it lies in the XY plane) with side length L. The force of a small piece of sheet on a 1kg point mass a distance $z$ above the origin:

$dF_{z} = \frac{G(\vec{dM}\cdot\vec{m})}{|r|^{2}}$

where $(\vec{dM}\cdot\vec{m}) = |dM||m|cos{\theta} $ and $|dM| = \frac{MdS}{A} = \sigma dxdy$ such that

$dF_{z} = \frac{G |dM||m|\cos{\theta}}{|r|^{2}}$
. We substitute $\cos{\theta}=\frac{z}{|r|}$ and $|r| = (x^2+y^2+z^2)^{1/2}$ then integrate to get:
$F_{z} = G m \sigma z \int \int_{-L/2}^{L/2} \frac{dx dy}{(x^2+y^2+z^2)^{3/2}}$

where $G=6.674\cdot10^{-11} \text{m}^3 \text{kg}^{-1} \text{s}^{-2}$ is Newton's gravitational constant and $\sigma$ is the mass per unit area of the sheet.

In [32]:
G = 6.674e-11
zi = 0 # m 
zf = 10 # m
L = 10 # m
a = -L/2
b = L/2
N = 100 # points
sigma = 5 # kg / m^2
m = 1 # kg

def dFdS(x,y,z):
    return (G*sigma*z*m)/(x**2+y**2+z**2)**(3/2)


def gquaddouble(function):
  
    xp, wp = gaussxwab(N,a,b)

    s = 0.0

    for i in range(N):
        for j in range(N):
            s += wp[i]*wp[j]*function(xp[i],xp[j],z)

    return s
 
Zvals = linspace(zi,zf,N)
Svals = list()
for z in Zvals:
    Svals.append(gquaddouble(dFdS))
plt.figure()
plt.title(r"$F_{z}$ due to plate (N) vs. z (m)")
plt.plot(Zvals,Svals,'b--')
plt.show()

There is a smooth curve, except at very small values of $z$, where the force drops off suddenly to zero. This drop is not a real effect, but an artifact of the way we have done the calculation. Since, having small values of $z$ does not constitute having small values of $|r|$ at every small piece of sheet with mass $dM$.

In [ ]: