Física Computacional

Semestre 1, 2024

David Cortés

1 Riemann integral
Compute the integral ˆ 4
sin (x/2) dx
0
using Riemann sums with the left hand, right hand and midpoint rules, for a number N of partitions. Compare
the solution to the analytical solution and plot the absolute error as a function of N , i.e. using different number
of partitions.

2 Velocities
In the on-line resources you will find a file called 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. The first few lines look like this: 0 0 1 0.069478 2 0.137694 3
0.204332 4 0.269083 5 0.331656 Write a program to do the following:
1. Read in the data and, using the trapezoidal rule, calculate from them the approximate distance traveled by
the particle in the x direction as a function of time. To read the data you can use different methods, such as
the open(...) Python function in its standard library, or Numpy’s loadtxt(...) function.
2. Extend your program to make a graph that shows, on the same plot, both the original velocity curve and
the distance traveled as a function of time.

3 Simpson’s Integration
´2
1. Write a program to calculate an approximate value for the integral 0 (x4 − 2x + 1)dx using Simpson’s rule
with 10 slices.
2. Run the program and compare your result to the known correct value of 4.4. What is the fractional error on
your calculation?
3. Modify the program to use an increasing number of slices. Note the improvement in the result. How do the
results compare with those from the trapezoidal rule with corresponding numbers of slices? Compute the
absolute error, and relative error (in a different plot), as a function of slices and plot the results for both the
trapezoidal and composite Simpson’s rule.

4 Simpson’s Integration 2
Consider the integral ˆ x
2
E(x) = e−t dt.
0
1. 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.
2. When you are convinced your program is working, extend it further to make a graph of E(x) as a function
of x.
Note that there is no known way to perform this particular integral analytically, so numerical approaches are the
only way forward.
5 The diffraction limit of a telescope
Our ability to resolve detail 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 λ, passes through the circular aperture of a telescope (which we’ll assume to have unit radius) and is
focused by the telescope in the focal plane, it produces not a single dot, but a circular 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
J1 (kr) 2
 
I(r) = ,
kr
where r is the distance in the focal plane from the center of the diffraction pattern, k = 2π/λ, and J1 (x) is a Bessel
function. The Bessel functions Jm (x) are given by
ˆ
1 π
Jm (x) = cos(mθ − x sin θ)dθ,
π 0

where m is a nonnegative integer and x ≥ 0.

1. Write a Python function J(m,x) that calculates the value of Jm (x) using Simpson’s rule with N = 1000
points. Use your function in a program to make a plot, on a single graph, of the Bessel functions J0 , J1 ,
and J2 as a function of x from x = 0 to x = 20.

2. Make a second program that makes a density plot of the intensity of the circular diffraction pattern of a
point light source with λ = 500 nm, in a square region of the focal plane, using the formula given above.
Your picture should cover values of r from zero up to about 1 µm.

Hint 1: You may find it useful to know that limx→0 J1 (x)/x = 1/2. Hint 2: The central spot in the diffraction
pattern is so bright that it may be difficult to see the rings around it on the computer screen. If you run into
this problem a simple way to deal with it is to use one of the other color schemes for density plots described in
Section 3.3. The “hot” scheme works well. For a more sophisticated solution to the problem, the imshow function
has an additional argument vmax that allows you to set the value that corresponds to the brightest point in the
plot. For instance, if you say “imshow(x,vmax=0.1)”, then elements in x with value 0.1, or any greater value, will
produce the brightest (most positive) color on the screen. By lowering the vmax value, you can reduce the total
range of values between the minimum and maximum brightness, and hence increase the sensitivity of the plot,
making subtle details visible. (There is also a vmin argument that can be used to set the value that corresponds
to the dimmest (most negative) color.) For this exercise a value of vmax=0.01 appears to work well.

6 Error estimation
´2
Write a program, or modify an earlier one, to once more calculate the value of the integral 0 (x4 − 2x + 1) dx from
Exercise (3), using the trapezoidal rule with 20 slices, but this time have the program also print an estimate of the
error on the result using the simplified calculation seen in class. To do this you will need to evaluate the integral
twice, once with N1 = 10 slices and then again with N2 = 20 slices. How does the error calculated in this manner
compare with a direct computation of the error as the difference between your value for the integral and the true
value of 4.4? Why do the two not agree perfectly?

7 Adaptive integration
Consider the integral ˆ 1 √
I= sin2 100x dx
0

1. Write a program that uses the adaptive trapezoidal rule method to calculate the value of this integral
to an approximate accuracy of ϵ = 10−6 (i.e., correct to six digits after the decimal point). Start with one
single integration slice and work up from there to two, four, eight, and so forth. Have your program print
out the number of slices, its estimate of the integral, and its estimate of the error on the integral, for each
value of the number of slices N , until the target accuracy is reached. (Hint: You should find the result is
around I = 0.45.)
 2. Now modify your program to evaluate the same integral using the Romberg integration technique
described in class. Have your program print out a triangular table of values, of all the Romberg estimates of
the integral. Calculate the error on your estimates using the expression for cm h2mi , and again continue the
calculation until you reach an accuracy of ϵ = 10−6 . You should find that the Romberg method reaches the
required accuracy considerably faster than the trapezoidal rule alone.

8 Adaptive integration 2
Write a program that uses the adaptive Simpson’s rule method to calculate the same integral as in problem 7,
again to an approximate accuracy of ϵ = 10−6 . Starting this time with two integration slices, work up from there
to four, eight, and so forth, printing out the results at each step until the required accuracy is reached. You should
find you reach that accuracy for a significantly smaller number of slices than with the trapezoidal rule calculation
in part (a) of Exercise 5.7, but a somewhat larger number than with the Romberg integration of part (b).

9 Legendre polynomials
´1
The Gaussian quadrature method to approximate the integral −1 W (x) f (x) dx, for the weight W (x) = 1,
determines the sample points as the roots (or zeroes) of the Legendre polynomials, depending on the number
of nodes. In addition, numerical weights depend on the Legendre polynomial derivatives. For this exercise,
1. Plot the first four Legendre polynomials Pi (x), i ∈ {1, 2, 3, 4}, and use a root finding function from Scipy or
Numpy function to find the zeroes of the polynomials. Compare your solutions with the values found in textbooks.
2. Find a method (for instance, interpolation) to find the first derivative of the Legendre polynomials for the
cases i = 2 and 3, and plot the results.

10 Heat capacity of a solid

Debye’s theory of solids gives the heat capacity of a solid at temperature T to be
 3 ˆ θD /T
T x4 e x
CV = 9V ρkB dx,
θD 0 (e − 1)2
x

where V is the volume of the solid, ρ is the number density of atoms, kB is Boltzmann’s constant, and θD is the
so-called Debye temperature, a property of solids that depends on their density and speed of sound.
1. Write a Python function cv(T) that calculates CV for a given value of the temperature, for a sample consisting
of 1000 cubic centimeters of solid aluminum, which has a number density of ρ = 6.022 × 1028 m−3 and a
Debye temperature of θD = 428 K. Use Gaussian quadrature to evaluate the integral, with N = 50 sample
points.

2. Use your function to make a graph of the heat capacity as a function of temperature from T = 5 K to
T = 500 K.

11 The Stefan–Boltzmann constant

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

ℏ ω3
I(ω) = .
4π 2 c2 (eℏω/kB T − 1)
This is the radiance expression given by Planck’s law. Here ℏ is Planck’s constant over 2π, c is the speed of light,
and kB is Boltzmann’s constant.
1. Show that the total energy per unit area radiated by a black body is given by Stefan-Boltzmann law,
4 T4 ˆ ∞
kB x3
W = 2 2 3 dx.
4π c ℏ 0 ex − 1

2. Write a program to evaluate the integral in this expression. Explain what method you used, and how accurate
you think your answer is.
 3. 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 = σT 4 ,
where σ 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. Compare the solution of the
integral to the analytical value of π 4 /15 and how well your method of choice approximates this integral.

12 Quantum uncertainty in the harmonic oscillator

In (dimensionless) 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
1 2
ψn (x) = p √ e−x /2 Hn (x),
2n n! π
for n = 0 . . . ∞, where Hn (x) is the nth Hermite polynomial. Hermite polynomials satisfy a relation somewhat
similar to that for the Fibonacci numbers, although more complex:

Hn+1 (x) = 2xHn (x) − 2nHn−1 (x).

The first two Hermite polynomials are H0 (x) = 1 and H1 (x) = 2x.
1. 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 x = −4 to x = 4. Hint: There is a function factorial in the math package that
calculates the factorial of an integer.
2. Make a separate plot of the wavefunction for n = 30 from x = −10 to x = 10. Hint: If your program takes
too long to run in this case, then you’re doing the calculation wrong—the program should take only a second
or so to run.
3. The quantum uncertainty of a particle in the nth level of a quantum harmonic oscillator can be quantified
by its root-mean-square position ⟨x2 ⟩, where
p

ˆ ∞
⟨x2 ⟩ = x2 |ψn (x)|2 dx.
−∞

Write a program that evaluates this integral using Gaussian quadrature on 50 and 100 points and then
calculates the uncertainty (i.e., the root-mean-square position of the particle) for a given value
p of n. Use
your program to calculate the uncertainty for n = 5. You should get an answer in the vicinity of ⟨x2 ⟩ = 2.3.
Note: The following code uses recurrence to compute the Fibonacci numbers in Python, which you might find
useful to compute Hermite polynomials:
def recur_fibo (n ) :
i f n <= 1 :
return n
else :
r e t u r n ( r e c u r _ f i b o ( n−1) + r e c u r _ f i b o ( n−2))

f o r i in range ( 1 0 ) :
print ( recur_fibo ( i ))

13 The gamma function

A commonly occurring function in physics calculations is the gamma function Γ(a), which is defined by the integral
ˆ ∞
Γ(a) = xa−1 e−x dx.
0

There is no closed-form expression for the gamma function, but one can calculate its value for given a by performing
the integral above numerically. You have to be careful how you do it, however, if you wish to get an accurate
answer.
1. Write a program to make a graph of the value of the integrand xa−1 e−x as a function of x from x = 0 to
x = 5, with three separate curves for a = 2, 3, and 4, all on the same axes. You should find that the integrand
starts at zero, rises to a maximum, and then decays again for each curve.

2. Show analytically that the maximum falls at x = a − 1.

3. Most of the area under the integrand falls near the maximum, so to get an accurate value of the gamma
function we need to do a good job of this part of the integral. We can change the integral from 0 to ∞ to
one over a finite range from 0 to 1 using the change of variables z = x/(1 + x), but this tends to squash the
peak towards the edge of the [0, 1] range and does a poor job of evaluating the integral accurately. We can
do a better job by making a different change of variables that puts the peak in the middle of the integration
range, around 1/2. We will use the change of variables:
x
z= .
c+x
For what value of x does this change of variables give z = 1/2? Hence what is the appropriate choice of
the parameter c that puts the peak of the integrand for the gamma function at z = 1/2? This is, find the
optimal c such that at z = 1/2 is the peak of the function.

4. Before we can calculate the gamma function, there is another detail we need to attend to. The integrand
xa−1 e−x can be difficult to evaluate because the factor xa−1 can become very large and the factor e−x very
small, causing numerical overflow or underflow, or both, for some values of x. Write xa−1 = e(a−1) ln x to
derive an alternative expression for the integrand that does not suffer from these problems (or at least not
so much). Explain why your new expression is better than the old one.

5. Now, using the change of variables above and the value of c you have chosen, write a user-defined function
gamma(a) to calculate the gamma function for arbitrary argument a. Use whatever integration method you
feel is appropriate. Test your function by using it to calculate and print the value of Γ( 32 ), which is known
√
to be equal to 12 π ≃ 0.886.

6. For integer values of a it can be shown that Γ(a) is equal to the factorial of a−1. Use your Python function to
calculate Γ(3), Γ(6), and Γ(10). You should get answers closely equal to 2! = 2, 5! = 120, and 9! = 362 880.