Computational Physics Lab 4

Basic Numerical Analysis Techniques

The goal of this lab is to explore basic numerical analysis techniques: function minimization and numerical integration.

Part 1. Optimizing Projectile Firing Angle

Note that numerical "function minimization" algorithms are applied to function maximization problems as well. Obviously, the problem of maximizing function f(x) is equivalent to the problem of minimizing −f(x). Historically, early optimization problems (brachistochrone, least action principle, cost management) involved minimization of the objective function, so the word "minimization" is commonly used in describing this type of problems, and most optimization software packages are designed to locate function minima.

In this part of the lab you will use the optimize package from SciPy to determine the optimal projectile firing angle in the presence of quadratic air drag. The goal is to shoot the projectile as far as possible. To determine the shooting distance, calculate the projectile trajectory using the 2^nd or the 4^th order Runge-Kutta ODE solving algorithm. Bracket the expected location of the maximum horizontal distance (as a function of the firing angle) and then apply Brent's parabolic interpolation minimization scheme to pinpoint the location of the maximum. Note that you will need to determine the projectile impact point with sufficient precision so that the projectile flight distance as a function of firing angle has a continuous second derivative (within the requested tolerance of Brent's method), otherwise parabolic interpolation becomes of no use.

It makes sense to determine the point where the projectile hits the ground in such a way that the precision of this point is at least as good as the precision with which the trajectory itself is calculated. For example, the global precision of the 4^th order Runge-Kutta metod is O(Δt⁴). In order to maintain a similar precision for the shooting distance, the last 4 points of the projectile trajectory should be interpolated using a 3^rd order polynomial which gives you y(x) up to O(Δx⁴) (this assumes that the ODE integration procedure is stopped as soon as the projectile's y coordinate gets below the ground level). The interpolation polynomial can be built using the polyfit function of the polynomial module from numpy. The root of this polynomial can then be determined using the roots function (you will need to pick the right root from the returned array).

You will need to make an appropriate choice of tolerance for the Brent's method. If the objective function is calculated with relative precision ε (you might want to estimate this precision by studying the projectile flight distance as a function of the ODE solver step size for a few different firing angles), Brent's method tolerance should not, as a rule of thumb, be less than √ε. This recommendation comes from noticing that, near the minimim located at some value b, the linear term in the Taylor's expansion of the objective function vanishes:

The second term will be negligible compared to the first (that is, will be factor of ε smaller and comparable with the noise) whenever

The reason for writing the right-hand side in this way is that, for most functions, the final square root is a number of order unity. For such functions, it is hopeless to ask for a bracketing interval of width less than √ε times its central value.

Use the following system parameters and initial conditions (they describe, roughly, a shot fired from the M16 rifle using the standard NATO 5.56 mm cartridge):

Projectile mass: 4.0 g
Projectile cross section area: 2.43 × 10^-5 m²
Initial projectile speed: 930 m/s
Initial elevation: 0.3 m
Air density: 1.2 kg/m³
The drag force expression is standard (see, e.g., Eq. 2.9 in Giordano & Nakanishi), but the dependence of the air drag coefficient on the projectile speed is rather complicated:

where M is the Mach number of the projectile: M = v/v_s, and v_s is the speed of sound in the air: 340 m/s. Note the strong characteristic jump in the value of C near M = 1. This jump limits the usefulness of high-order ODE integration schemes; in fact, a scheme with an adaptive step size would work best for this application. As before, you can use modules v3.py, cpode.py, and cpforces.py to assist you with building the force model for the projectile and with running the ODE solver.

We will encounter numerous problems of multivariate function minimization in the future, when we consider fitting a model with multiple parameters to a set of data points.

Part 2. Calculating the Muon Lifetime

Muon is an elementary particle with properties very similar to those of the electron but with much higher mass: it is about 200 times heavier than the electron. To produce a muon, one needs roughly 106 MeV of energy. This energy is far above typical energies released in chemical and nuclear reactions. On Earth, muons are produced in the reactions between energetic galactic cosmic rays and atmospheric molecules. The primary cosmic rays collide with the nuclei of air molecules and produce a shower of particles that include protons, neutrons, pions (both charged and neutral), kaons, photons, electrons and positrons. These secondary particles then undergo electromagnetic and nuclear interactions to produce yet additional particles in a cascade process. Figure below illustrates the general idea:

Muons are produced in the decays of charged pions in the cascade. About 10,000 muons hit every square meter of the Earth's surface per minute. Muons themselves decay in flight into muon neutrino, electron, and electron antineutrino:

Antimuons (positively charged muons) decay into muon antineutrino, positron, and electron neutrino (that is, all decay products are changed into their respective antiparticles). Measurement of the muon lifetime using scintillator slabs and photomultipliers is one of the classic experiments in the TTU advanced physics labs course.

In this lab, however, we want to calculate the muon lifetime using numerical analysis techniques. According to the theory of particle interactions (which is beyond the scope of this note), the rate for a muon to decay into the electron and two neutrinos with particular values of initial and final momenta consistent with energy-momentum conservation is

The rate is called |M|² (matrix element squared) according to the particle physics convention. G_F is the Fermi constant of weak interaction: G_F = 1.16637 × 10^-5 GeV^-2. Terms (P_a ⋅ P_b) stand for the relativistic scalar products of the particle 4-momenta:

Energies are measured here in units of GeV and 3-momenta are measured in units GeV/c.

This matrix element is not exact: it is given here within first-order perturbation theory. However, this should be good enough already to estimate the muon lifetime with 1% precision.

To get the total reaction rate (usually called decay width, Γ) one has to integrate |M|² over all possible states of the final particles in the reaction, as well as over those states of the initial particle which are relevant to the question at hand. Effectively, one assumes that the density of states in the 4-dimensional relativistic momentum space is constant (this assumption is in good agreement with experimental results). The exact value of this constant depends on the conventions chosen for the units of action and for the wavefunction normalization. In the next formula, this constant is assumed to be (2 π)^-1 for every integration dimension (this corresponds to action measured in units of h/(2 π) where h is the Plank's constant):

Only non-negative energy values are used in this integration. The function T(...) embodies out notion of which values of final 4-momenta are possible and which initial values are relevant. In particular, every particle in the reaction has a definite mass. This means T(...) contains the mass constraint term

for each particle in the reaction. This is just a mathematically consistent way to assign definite masses to relativistic particles inside such an integral. The coefficient 2π appears here (and in all other delta functions in the momentum space) due to the normalization convention adopted for the relativistic wavefunctions. Also, the function T(...) includes the term

which imposes the energy-momentum conservation on the reaction.

Due to the relativistic time dilation effect (time passes slower for a muon in flight), the observed muon decay rate depends on its speed. This is clearly undesirable for our theoretical calculation because we would like to study the intrinsic property of the muon rather than how it behaves in a particular experimental situation. Because of this, we will calculate the decay rate in the system of coordinates in which the muon is at rest. This rest system is selected by introducing the term

into the T(...) function. The muon decay rate in any other coordinate system (in the lab system in particular) can be obtained from the rest frame decay rate by using the standard relativistic time dilation formula. The lifetime is related to the decay rate by τ = 1/Γ. The lifetime of the particle in its own rest frame is called proper lifetime.

Finally, we have the 16-dimentional integral for the total decay rate Γ in which the function T(...) contains 11 delta functions. Brute-force calculation of this integral is hopeless -- there is simply no chance that a point chosen randomly or on some grid in 16 dimensions will satisfy all the constraints imposed by the delta functions. We must reduce the integral to 5 dimensions and get rid of all delta functions in the process. Fortunately, the procedure for doing this is relatively straightforward. It is based on the standard variable transformation formula:

where J is the Jacobian of the transformation:

One chooses the new set of variables u₀, ..., u_n in such a way that one or more of these variables look exactly like the arguments of the delta functions inside the integrand. Then this new variable can be integrated over. For example, if the integrand contains δ(u₀), u₀can be integrated over. After this, the only trace of δ(u₀)remaining in the integral is the requirement that the rest of the integrand must be evaluated assuming u₀= 0.

To begin with, let's change the variables so that only one integration variable changes: instead of E_e- we introduce

The Jacobian of this transformation is just

After this transformation, one can integrate away the mass constraint term for e^-. The net effect of the e^- mass constraint on the Γ formula is the following:

where E_e- is no longer an independent variable and must be calculated using

In this formula we assume that the electron mass is measured in units GeV/c², momentum in GeV/c, and energy in GeV. Obviously, all other mass constraint terms can be integrated out in a similar manner. The delta function which restricts our consideration to the rest system of the decaying muon is in the form which can be integrated out without any special change of variables. Thus, the only delta function which remains is the 4-dimensional delta function responsible for the energy-momentum conservation:

Note that, since E_μ- has to be calculated for the muon at rest, E_μ- (in units of GeV) is just m_μ- (in units of GeV/c²).

In principle, the remaining delta function can be removed in a similar manner, using variable transformations. There is, however, a more efficient technique for doing this called phase space splitting. In this technique, we assume that the decay proceeds in two stages: first, the muon decays into an electron and a fictitious compound particle, and then the compound particle decays into two neutrinos (it does not matter which two of the three decay products are used to build the compound particle). Since we do not know the compound particle mass, we will integrate over it. Formally, this program can be carried out by inserting the compound particle with 4-momentum P_c into the formula for Γ like this:

Here comes the phase space splitting trick: let's use the square of the mass of the compound particle as another independent variable. Then, under the integral, the following is formally true:

The delta function in this formula looks just like the delta function in the mass constraints, and we already know how to eliminate those. We can combine this formula with the previous one by noticing that

so that finally

Now, the integral for Γ contains (among other things) two very similar terms:

and

Each of these terms describes integration over 4-momentum space in a decay of one particle into two. Moreover, each of these terms is written in a relativistically invariant form -- it does not matter which system of coordinates is used to evaluate these integrals. This is because in the expression

the δ(P² - m²) term is manifestly invariant because it involves only invariant quantities, and d⁴P term is invariant because the Jacobian of any Lorentz transformation is 1. Because these terms are invariant, one can choose any coordinate system to evaluate these integrals. In particular, the choice of the system in which the decaying particle is at rest results in an especially simple treatment of the energy-momentum conservation delta function. In that system, integration over the 3-momentum of one of the particles (let say, ν_μ) results in elimination of 3 out of 4 delta functions (those which correspond to the conservation of 3-momenta). To eliminate the energy conservation delta function, one performs the change of variables using the magnitude of the other particle's momentum. For example, one introduces the variable

Since the 3-momentum of the muon neutrino is no longer an independent quantity (it has been integrated over, and the momentum conservation requires that it is equal in magnitude but opposite in direction to the 3-momentum of the electron antineitrino), the Jacobian of this transformation is

where we are using simplified notation for the momentum magnitude:

We can now write

Integration over u₁ together with the remaining delta function gives 1. Finally,

where the right side (including dΩ) must be evaluated in the rest system of the compound particle. In that system E_c = m_c. It is also customary to solve the energy-momentum conservation equations and express the momentum of the decay products in the rest system of the decaying particle using invariant quantities:

where the function λ is defined by

The treatment of the muon decay into the electron and the compound particle is analogous. Collecting everything together, the muon decay rate becomes

As expected, this is a 5-dimensional integral. Note that integrations over solid angles are performed in different coordinate systems which move w.r.t. each other. In practice, this does not result in any serious problem -- one builds the integration grid (or generates random numbers in case Monte Carlo method is used) in one system and then transforms that grid into another system using standard Lorentz transformation formulae.

When decays into more than three particles are considered, the phase space splitting trick can be invoked recursively, thus representing the whole reaction as a chain in which one has to deal with two-body particle decays only. This technique is widely used in calculation of rates of various processes in particle physics by Monte Carlo method. Each step in the chain consists of simulating a random decay of a fictitious compound particle (or, for the first step in the chain, the real initial particle) into two daughter particles in the rest system of the parent. Then the two daughters are boosted back into the lab system, and the procedure is repeated for the decays of the daughters. The masses of the compounds are generated randomly within their respective limits.

So far, the discussion of the integration method has been completely generic. It can be applied to an arbitrary decay of any particle into three daughter particles -- the only difference would be in the particle masses and in the expression for |M|². Even though it is possible to simplify the integral further for the muon decay (using the fact that neutrino masses are extremely small and that |M|²for this decay is very simple), we will stop algebraic reasoning here and solve the general problem numerically. The advantage of this approach is that it can be easily reused in the future in some other similar theoretical calculation.

5-d integration is somewhere on the border where Gaussian quadratures or repeated 1-d integration are still feasible. However, I suggest applying Monte Carlo and Quasi-Monte Carlo methods to this problem in order to get a feeling of how these techniques work. Please write your program in such a manner that switching from Monte Carlo to Quasi-Monte Carlo does not require any effort. Probably, the simplest way to do this is to assume that you have a sequence of 5-d (quasi) random vectors. Each vector has components between 0 and 1. You can pick the random numbers from that sequence, one vector per one integral evaluation point. You should be able to generate the sequence before running your main integration program by employing different generators of pseudo-random and quasi-random numbers.

A good way to generate the compound mass squared is to pick it randomly and uniformly between its smallest and largest possible values:

The integral must be multiplied by the range of the random variable (if you are not sure why this is so, please review the Monte Carlo part of the numerical integration lecture). To transform random numbers on the (0, 1) range to an arbitrary range (a, b) use the following formula: r_{(a, b)} = a + r_{(0, 1)}*(b - a). This transformation maps 0 into a and 1 into b. All numbers in between are mapped in a linear fashion.

A good way to generate random solid angles is to generate points uniformly on the surface of a sphere. This can be achieved by generating cosine of the polar angle, cos(θ), uniformly between −1 and 1, and the azimuthal angle, φ, uniformly between 0 and 2π (this works because dΩ = −dcos(θ) dφ). Each random solid angle contributes the range factor 4π into the integral. The angles generated randomly in the rest system of the compound particle should be transformed into the rest system of the muon.

With all the range factors included, the muon decay rate becomes

where the averaging is performed over the sequence of random values of compound mass squared and solid angles generated as described above.

Using this formula, calculate the muon lifetime with 1% precision using the pseudo-random numbers provided by the standard random module. Use the following values for the particle masses:

Mass of the muon: 105.66 × 10^-3 GeV/c².
Mass of the electron: 0.511 × 10^-3 GeV/c².
Masses of the neutrinos: 0.0 GeV/c².

The lifetime 1/Γ which you will obtain in this calculation will be expressed in units GeV^-1. To convert this time into seconds, multiply your result by the conversion coefficient: h/(2 π) expressed in units GeV⋅s. The numerical value of the conversion coefficient is 6.582119 × 10^-25. The experimentally measured value of the muon lifetime, which is pretty close but not exactly equal to the number you expect to obtain, is (2.19703 ± 0.00004) × 10^-6 s.

The module rk.py provides basic support for calculations in relativistic kinematics on which you can base your integration program. It includes 4-vector and Lorentz transformation (boost) classes, as well as an implementation of the λ function.

The sequence of quasi-random numbers can be obtained in the following way: download the files quasi_random.f, quasi_random.pyf, exor.c, sobol.py, and niederreiter.py and then run

f2py -c quasi_random.pyf quasi_random.f exor.c

This will create a compiled python module "quasi_random" in your directory (the module file will be named "quasi_random.so"). After this, you should be able to generate a sequence of N 5-dimensional Sobol pseudo-random vectors in the following manner:

import sobol
sobol.init(5)
sequence = []
for i in xrange(N):
sequence.append(sobol.next())

Note that it only makes sense to use a set of N = 2^k Sobol numbers where k is an integer. Run your code using quasi-random instead of pseudo-random numbers for k = 6, 7, 8, 9, ..., 16, and plot the result as a function of k. Can you comment on the algorithm convergence? Note that you do not need to reinitialize the generator for all these different k values. Just use k = 16 to generate the quasi-random sequence of N = 65536 vectors, but along the way store the values of the integral calculated using 64, 128, 256, etc. vectors.

Please submit your lab reports before 2 pm 02/19/2008. I suggest including the scripts which solve the problems described in parts 1 and 2 of the lab, computation results, and some description of the important choices which you had to make in the course of writing your programs and in selecting algorithm parameters.