import sys
print(sys.executable)

/home/docs/checkouts/readthedocs.org/user_builds/physics-555-quantum-technologies/conda/latest/bin/python3

%matplotlib inline
import numpy as np, matplotlib.pyplot as plt
import mmf_setup;mmf_setup.nbinit()

This cell adds /home/docs/checkouts/readthedocs.org/user_builds/physics-555-quantum-technologies/checkouts/latest/src to your path, and contains some definitions for equations and some CSS for styling the notebook. If things look a bit strange, please try the following:

  • Choose "Trust Notebook" from the "File" menu.
  • Re-execute this cell.
  • Reload the notebook.

Assignment 0: Introduction#

Assignment Preparation#

Complete the following in one of two ways:

  1. A simple brute force solution that you completely understand but which may not be fully accurate.

  2. A more sophisticated solution that is accurate to close to machine precision.

For the second method, please feel free to use any tools in the NumPy or SciPy libraries. (You may use mpmath to check, but don’t use it as part of your solution.)

Series#

Numerically check the following formula:

\[ \sum_{n=1}^{M} n = \frac{M(M+1)}{2} \]
\[ \sum_{n=1}^{\infty} \frac{1}{n^{2}} = \zeta(2) = \frac{\pi^2}{6} \]

where \(\zeta(s)\) is the Riemann zeta function.

Integrals#

Numerically check the following integrals for all viable values of \(p\) (both positive and negative):

\[ \int_0^1 x^p \d{x} = \frac{1}{p+1} \]

What about infinite limits?

\[ \int_{-\infty}^{\infty} e^{-x^2} \d{x} = \sqrt{\pi} \]

This one is probably challenging to do accurately:

\[ \int_{-\infty}^{\infty} e^{-x^2}\sin\frac{1}{x^2}\d{x} = \sqrt{\pi}e^{-\sqrt{2}}\sin(\sqrt{2}) \]

This should be easy if you could do the previous one:

\[\begin{split} \int_0^{\infty} xe^{-x^2}\sin^2\frac{1}{x}\d{x} = \frac{\sqrt{\pi}}{4} G_{14}^{13}\left( \left. \begin{matrix} \tfrac{1}{2}\\ \tfrac{1}{2} & \tfrac{1}{2} & 0 & -\tfrac{1}{2} \end{matrix} \right| z=1 \right) = 0.32006330909018418888\cdots, %37810082287661243051390948375855688739690521501945667936210716620524245639515708 \end{split}\]

where \(G_{mn}^{pq}\bigl(\begin{smallmatrix}\vect{a}_p\\ \vect{b}_q\end{smallmatrix}\big|z\bigr)\) is the Meijer G-function.

Roots#

Find all solutions \(x\) to the following equations:

Polynomials#

\[ x^2 - (1+\epsilon)x + \epsilon = 0, \qquad x = \{1, \epsilon\} \]

for \(\epsilon \in \{1, 10^{-10}, 10^{-20}\}\).

\[ x^4 + 3x^2 - 6x + 10 = 0, \qquad x = \{1\pm\I, -1\pm 2\I\} \]

Lambert W function#

\[ xe^x = w, \qquad x = L_0(1) \]

for \(w = 1\).

The function \(w = W_k(x)\) is the Lambert W function.

More Details#

Write a function phys_555.assignment_0.quadratic_equation() which returns the roots of the equation

\[ ax^2 + bx + c = 0. \]

For example:

\[ x^2 - 3x + 2 = (x-2)(x-1) \]

so we expect

from phys_555.assignment_0 import quadratic_equation
np.allclose(quadratic_equation(a=1, b=-3, c=2), [1, 2])

---------------------------------------------------------------------------
ModuleNotFoundError                       Traceback (most recent call last)
Cell In[3], line 1
----> 1 from phys_555.assignment_0 import quadratic_equation
      2 np.allclose(quadratic_equation(a=1, b=-3, c=2), [1, 2])

ModuleNotFoundError: No module named 'phys_555.assignment_0'

Note: if you attempt to blindly use the quadratic formula as is done in the default version, you will encounter errors when \(b \approx \pm \sqrt{b^2 - 4ac}\) because the two terms can cancel. This is the main source of error associated with floating point comptations.

  1. When will this become a problem?

  2. How can you overcome this issue?

  3. How will you test your code to make sure it works well?

The goal should be for every function to return an answer that has a relative error comparable to the machine precision of the computer: sometimes called \(\epsilon = \)eps. This is not always possible if the problem is ill-conditioned.

for epsilon in [1, 1e-10, 1e-20]:
    roots = quadratic_equation(a=1, b=-(1+epsilon), c=epsilon)
    exact_roots = [1, epsilon]
    print(f"epsilon = {epsilon}")
    print(f"roots =       {roots}")
    print(f"exact roots = {exact_roots}")
    print(f"Roots allclose? {np.allclose(roots, exact_roots)}")
    print()

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[4], line 2
      1 for epsilon in [1, 1e-10, 1e-20]:
----> 2     roots = quadratic_equation(a=1, b=-(1+epsilon), c=epsilon)
      3     exact_roots = [1, epsilon]
      4     print(f"epsilon = {epsilon}")

NameError: name 'quadratic_equation' is not defined

Floating Point Numbers#

Readings:

  • [Gezerlis, 2020] Chapter 2. Try some of the “experiments” Alex suggests.

  • [Goldberg, 1991] “What Every Computer Scientist Should Know About Floating-Point Numbers.” This is a rather technical, but complete account about floating point numbers, and contains a detailed analysis of this assignment.

We can use NumPy to see the properties of the floating point numbers. The defult float in python is the IEEE double-precision floating point number which has \(64 = 52 + 11 + 1\) bits. 52 of these are used for the mantissa, 11 for the exponent and 1 for the sign. Roughly, the relative error in a floating point number is eps=\(\epsilon = 1/2^{52}\approx 2.22\times 10^{-16}\) (16 digits of precision in decimal), with an exponent that can range from \(\pm 2^{10} = \pm 1024\) with a smallest value of tiny\(\approx 2^{-1024} \approx 5.56\times 10^{-309}\) and largest value of max\(=2^{1024}\approx 1.798\times 10^{208}\).

Try to compute these numbers! In particular, 2**1024.0 causes an overflow… how can you compute \(2^{1024}\approx 1.798\times 10^{208}\) using floating point numbers? Hint:

\[ 2^n = a\times 10^m, \qquad n\log_{10}(2) = \log_{10}(a) + m. \]

The actual values are a little different as we see below (e.g. tiny\(=2^{-1022}=2.225\times 10^{-308}\)) because of some subtleties in the actual implementation of the IEEE standard.

import numpy as np
print(np.finfo(float))

Machine parameters for float64
---------------------------------------------------------------
precision =  15   resolution = 1.0000000000000001e-15
machep =    -52   eps =        2.2204460492503131e-16
negep =     -53   epsneg =     1.1102230246251565e-16
minexp =  -1022   tiny =       2.2250738585072014e-308
maxexp =   1024   max =        1.7976931348623157e+308
nexp =       11   min =        -max
smallest_normal = 2.2250738585072014e-308   smallest_subnormal = 4.9406564584124654e-324
---------------------------------------------------------------

Density Matrices#

In quantum mechanics, density matrices with special properties play an important role. A density matrix \(\mat{\rho}\) has the following properties:

  • Hermtian: \(\mat{\rho} = \mat{\rho}^\dagger\). Real eigenvalues \(\lambda_n\) and orthonormal eigenvectors \(\braket{m|n} = \delta_{mn}\).

  • Positivity: \(\mat{\rho} \geq 0\). Non-negative eigenvalues \(\lambda_n \geq 0\).

  • Trace: \(\Tr\mat{\rho} = \sum_{n} \lambda_n = 1\).

We have also expressed these conditions in terms of properties of the eigenvalues \(\lambda_n\) and eigenvectors \(\ket{n}\).

  1. The most general 2×2 hermitian matrix \(\mat{H}\) can be expressed in terms of 4 arbitrary real numbers \(a\), \(b\), \(c\), \(d\):

    \[\begin{split} \newcommand{\pmat}[1]{\begin{pmatrix}#1\end{pmatrix}} \mat{H} = \pmat{a & b+ic\\b-ic & d}. \end{split}\]

    Hence the space of hermitian matrices is isomorphic to the 4-dimensional space\(\mathbb{R}^4\).

    What is the most general form of a 2×2 density matrix? What dimension is the space of such matrices?

  2. What is the dimension of the space of \(N\times N\) density matrices?