Integration

Exercises taken from

5.3

Consider the integral

E(x) = \int_0^x e^{-t^2} dt.

Write a program to calculate E(x) for values of x from 0 to 3 in steps of 0.1. Choose for yourself what method you will use for performing the integral and a suitable number of slices.
When you are convinced your program is working, extend it further to make a graph of E(x) as a function of x. If you want to remind yourself of how to make a graph, you should consult Section 3.1, starting on page 88. Note that there is no known way to perform this particular integral analytically, so numerical approaches are the only way forward.

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 \big({T\over{\theta_D}}\big)^3 \int_0^{\theta_D/T} {x^4 e^4\over{(e^x - 1)^2}}

where V is the volume of the solid, \rho is the number density of atoms, k_B is 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. Write a Python function cv(T) that calculates C_V for a given value of the temperature, for a sample consisting of 1000 cubic centimeters of solid aluminum, which has a number density of \rho=6.022\times10^{28}\,\mathrm{m}^{-3} and a Debye temperature of \theta_D=428\,K. Use Gaussian quadrature to evaluate the integral, with N=50 sample points. Use your function to make a graph of the heat capacity as a function of temperature from T=5K to T=500K.

5.11 Plane  Wave

Suppose a plane wave, such as light or a sound wave, is blocked by an object with a straight edge, represented by the solid line at the bottom of this figure:
edge

The wave will be diffracted at the edge and the resulting intensity at the position (x,z) marked by the dot is given by near-field diffraction theory to be

I = \frac{I_0}{8} \big([2C(u) + 1]^2 + [2S(u) + 1]^2 \big)

where I_0 is the intensity of the wave before diffraction and

u = x \sqrt{2\over{\lambda z}},\,\,\,\,\,\, C(u) = \int_0^u cos \frac{1}{2} \pi t^2 dt,\,\,\,\,\,\, S(u) = \int_0^u sin \frac{1}{2} \pi t^2 dt .

Write a program to calculate I/I_0 and make a plot of it as a function of x in the range -5m to 5m for the case of a sound wave with wavelength \lambda=1\,m, measured z=3m past the straight edge. Calculate the integrals using Gaussian quadrature with N=50 points. You should find significant variation in the intensity of the diffracted sound—enough that you could easily hear the effect if sound were diffracted, say, at the edge of a tall building.

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 $\Iatex I(\omega) d\omega$, where

I(\omega) = {\hbar \over{4 \pi^2 c^2}} {\omega^3 \over{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.  Show that the total energy per unit area radiated by a black body is

W = {{k_B^4 T^4}\over{4\pi^2 c^2 \hbar^3}} \int_0^{\infty} {x^3\over{e^x -1}} dx .

Write a program to evaluate the integral in this expression. Explain what method you used, and how accurate you think your answer is. Even before Planck gave his theory of thermal radiation around the turn of the 20th century, it was known that the total energy W given
off by a black body per unit area per second followed Stefan’s law: W = \sigma T^4, where \sigma is the Stefan–Boltzmann constant. Use your value for the integral above to compute a value for the Stefan–Boltzmann constant (in SI units) to three significant figures.
Check your result against the known value, which you can find in books or on-line. You should get good agreement.

5.14 Gravitational pull of a uniform sheet

A uniform square sheet of metal is floating motionless in space:
plate

The sheet is 10m on a side and of negligible thickness, and it has a mass of 10 metric tonnes.
Consider the gravitational force due to the plate felt by a point mass of 1kg a distance z from the center of the square, in the direction perpendicular to the sheet, as shown above. Show that the component of the force along the z-axis is

F_z = G \sigma z \int \int_{-L/2}^{L/2} {dx dy \over{(x^2 + y^2 + z^2)^{3/2}}}

where G=6.674\times10^{-11}\,\mathrm{m}^3\,\mathrm{kg}^{-1}\,\mathrm{s}^{-2} is
Newton’s gravitational constant and \sigma is the mass per unit area of the sheet.
Write a program to calculate and plot the force as a function of z from z=0 to z=10m. For the double integral use (double) Gaussian quadrature, as in Eq. (5.82), with 100 sample points along each axis.
You should see a smooth curve, except at very small values of z, where the force should drop off suddenly to zero. This drop is not a real effect, but an artifact of the way we have done the calculation. Explain briefly where this artifact comes from and suggest a strategy to
remove it, or at least to decrease its size. This calculation can thought of as a model for the gravitational pull of a galaxy. Most of the mass in a spiral galaxy (such as our own Milky Way)
lies in a thin plane or disk passing through the galactic center, and the gravitational pull exerted by that plane on bodies outside the galaxy can be calculated by just the methods we have employed here.