Computational_Physics_Problem_Solving_with_Python_Landau/CP06.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# *Chapter 6*<br> Matrix Computing\n",
    "\n",
    "| | | |\n",
    "|:---:|:---:|:---:|\n",
    "| ![image](Figs/Cover.png)|[From **COMPUTATIONAL PHYSICS**, 3rd Ed, 2015](http://physics.oregonstate.edu/~rubin/Books/CPbook/index.html) <br>RH Landau, MJ Paez, and CC Bordeianu (deceased) <br>Copyrights: <br> [Wiley-VCH, Berlin;](http://www.wiley-vch.de/publish/en/books/ISBN3-527-41315-4/) and [Wiley & Sons, New York](http://www.wiley.com/WileyCDA/WileyTitle/productCd-3527413154.html)<br>  R Landau, Oregon State Unv, <br>MJ Paez, Univ Antioquia,<br> C Bordeianu, Univ Bucharest, 2015.<br> Support by National Science Foundation.|![image](Figs/BackCover.png)|\n",
    "    \n",
    "**6 Matrix Computing**<br>\n",
    "[6.1 Problem 3: N-D Newton-Raphson; Two Masses on a String](#6.1)<br>\n",
    "&nbsp;&nbsp;&nbsp;&nbsp;[6.1.1 Theory: Statics](#6.1.1)<br>\n",
    "&nbsp;&nbsp;&nbsp;&nbsp;[6.1.2 Algorithm: Multidimensional Searching](#6.1.2)<br>\n",
    "[6.2 Why Matrix Computing?](#6.2)<br>\n",
    "[6.3 Classes of Matrix Problems (Maths)](#6.3)<br>\n",
    "&nbsp;&nbsp;&nbsp;&nbsp;6.3.1 Practical Matrix Computing](#6.3.1)<br>\n",
    "[6.4 Python Lists as Arrays](#6.4)<br>\n",
    "[6.5 Numerical Python (NumPy) Arrays](#6.5)<br>\n",
    "&nbsp;&nbsp;&nbsp;&nbsp;[6.5.1 NumPy’s linalg Package](#6.5.1)<br>\n",
    "[6.6 Exercise: Testing Matrix Programs](#6.6)<br>\n",
    "&nbsp;&nbsp;&nbsp;&nbsp;[6.6.1 Matrix Solution of the String Problem](#6.6.1)<br>\n",
    "&nbsp;&nbsp;&nbsp;&nbsp;[6.6.2 Explorations](#6.6.2)<br>\n",
    "\n",
    "*This chapter examines various aspects of computing with matrices, and\n",
    "in particular applications of the Python linear algebra packages.\n",
    "Because these packages are optimized and robust, we strongly recommend\n",
    "their use even if your programs are small (small programs often grow\n",
    "big). The two-mass-on-a-string problem is formulated as a matrix\n",
    "problem, and extends the Newton-Raphson search technique to be discussed\n",
    "in [Chapter 7](CP07.ipynb).*\n",
    "\n",
    "** This Chapter’s Lecture, Slide Web Links, Applets & Animations**\n",
    "\n",
    "| | |\n",
    "|---|---|\n",
    "|[All Lectures](http://physics.oregonstate.edu/~rubin/Books/CPbook/eBook/Lectures/index.html)|[![anything](Figs/RHLlectureMod4.png)](http://physics.oregonstate.edu/~rubin/Books/CPbook/eBook/Lectures/index.html)|\n",
    "\n",
    "| *Lecture (Flash)*| *Slides* | *Sections*|*Lecture (Flash)*| *Slides* | *Sections*|  \n",
    "|- - -|:- - -:|:- - -:|- - -|:- - -:|:- - -:|\n",
    "|[Matrix Compute](http://physics.oregonstate.edu/~rubin/Books/CPbook/eBook/Lectures/Modules/Matrices/Matrices.html)|[pdf](http://physics.oregonstate.edu/~rubin/Books/CPbook/eBook/Lectures/Slides/Slides_NoAnimate_pdf/Matrix.pdf)| 8.1 |[Matrix Compute II](http://physics.oregonstate.edu/~rubin/Books/CPbook/eBook/Lectures/Modules/Matrices/Matrices.html)|[pdf](http://physics.oregonstate.edu/~rubin/Books/CPbook/eBook/Lectures/Slides/Slides_NoAnimate_pdf/Matrix2.pdf)|8.4 |\n",
    "|[Trial and Error Searching](http://physics.oregonstate.edu/~rubin/Books/CPbook/eBook/Lectures/Modules/Searching/Searching.html)|[pdf](http://physics.oregonstate.edu/~rubin/Books/CPbook/eBook/Lectures/Slides/Slides_NoAnimate_pdf/Trial_Err.pdf)|7.7-7.10|[N-Dimensional Searching](http://physics.oregonstate.edu/~rubin/Books/CPbook/eBook/Lectures/Modules/NdimSearch/NdimSearch.html)|[pdf](http://physics.oregonstate.edu/~rubin/Books/CPbook/eBook/Lectures/Slides/Slides_NoAnimate_pdf/NDsearch.pdf)| 8.2.2|\n",
    "|[Interpolation](http://physics.oregonstate.edu/~rubin/Books/CPbook/eBook/Lectures/Modules/Interpolation/Interpolation.html)|[pdf](http://physics.oregonstate.edu/~rubin/Books/CPbook/eBook/Lectures/Slides/Slides_NoAnimate_pdf/Interp.pdf)|8.5 |[Least Square Fitting](http://physics.oregonstate.edu/~rubin/Books/CPbook/eBook/Lectures/Modules/Fitting/Fitting.html)| [pdf](http://physics.oregonstate.edu/~rubin/Books/CPbook/eBook/Lectures/Slides/Slides_NoAnimate_pdf/LeastSqFit.pdf)| 8.7 |\n",
    "|Applet: [Lagrange Interpolation](http://science.oregonstate.edu/~rubin/Books/CPbook/eBook/Applets/index.html)|-| 8.5 | Applet: [Spline Interpolation](http://science.oregonstate.edu/~rubin/Books/CPbook/eBook/Applets/index.html)|-| 8.5 |\n",
    "\n",
    "![image](Figs/Fig6_1.png) \n",
    "\n",
    "**Figure 6.1** Two weights connected by three\n",
    "pieces of string and suspended from a horizontal bar of length *L*. The\n",
    "lengths are all known, but the angles and the tensions in the strings\n",
    "are to be determined."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 6.1  Problem 3: N-D Newton-Raphson; Two Masses on a String<a id=\"6.1\"></a>\n",
    "\n",
    "\n",
    "**Problem:** Two weights (*W*<sub>1</sub>, *W*<sub>2</sub>)=(10, 20) are\n",
    "hung from three pieces of string with lengths\n",
    "(*L*<sub>1</sub>, *L*<sub>2</sub>, *L*<sub>3</sub>)=(3, 4, 4) and a\n",
    "horizontal bar of length *L* = 8 (Figure 6.1). Find the angles assumed\n",
    "by the strings and the tensions exerted by the strings.\n",
    "\n",
    "In spite of the fact that this is a simple problem requiring no more\n",
    "than first-year physics to formulate, the coupled transcendental\n",
    "equations that result are just about inhumanely painful to solve\n",
    "analytically\\[*Note:* Almost impossible anyway, as L. Molnar has\n",
    "supplied me with an analytic solution.\\]. However, we will show you how\n",
    "the computer can solve this problem, but even then only by a\n",
    "trial-and-error technique with no guarantee of success. Your **problem**\n",
    "is to test this solution for a variety of weights and lengths and then\n",
    "to extend it to the three-weight problem (not as easy as it may seem).\n",
    "In either case check the physical reasonableness of your solution; the\n",
    "deduced tensions should be positive and similar in magnitude to the\n",
    "weights of the spheres, and the deduced angles should correspond to a\n",
    "physically realizable geometry, as confirmed with a sketch such as\n",
    "Figure 2.2 B. Some of the exploration you should do is to see at what\n",
    "point your initial guess gets so bad that the computer is unable to find\n",
    "a physical solution.\n",
    "\n",
    "![image](Figs/Fig6_2.png)\n",
    "\n",
    "**Figure 6.2** A free body diagram for one\n",
    "weight in equilibrium. Balancing the forces in the *x* and *y*\n",
    "directions for all weights leads to the equations of static equilibrium.\n",
    "\n",
    "### 6.1.1  Theory: Statics<a id=\"6.1.1\"></a>\n",
    "\n",
    "We start with the geometric constraints that the horizontal length of the\n",
    "structure is *L* and that the strings begin and end at the same height\n",
    "(Figure 6.1):[[xml]](http://physics.oregonstate.edu/~rubin/Books/CPbook/eBook/xml/6.1.xml)\n",
    "\n",
    "$$\\begin{align}\n",
    "\\tag*{6.1}\n",
    "L_1\\cos\\theta_{1} + L_2\\cos\\theta_{2}+ L_3\\cos\\theta_{3} & = L,\\\\\n",
    "L_1\\sin\\theta_{1} + L_2\\sin\\theta_{2} -L_3\\sin\\theta_{3} & = 0,\\tag*{6.2}\\\\\n",
    "\\sin^2\\theta_1+\\cos^2\\theta_1 & =  1,\\tag*{6.3}\\\\\n",
    "\\sin^2\\theta_2+\\cos^2\\theta_2 & =  1,\\tag*{6.4}\\\\\n",
    "\\sin^2\\theta_3+\\cos^2\\theta_3 & =  1.\\tag*{6.5}\\end{align}$$\n",
    "\n",
    "Observe that the last three equations include trigonometric identities as\n",
    "independent equations because we are treating sin*θ* and cos*θ* as\n",
    "independent variables; this makes the search procedure easier to implement.\n",
    "The basics physics says that because there are no accelerations, the sum of the\n",
    "forces in the horizontal and vertical directions must equal zero (Figure 6.2):\n",
    "\n",
    "$$\\begin{align} T_{1}\\sin\\theta_{1} - T_{2}\\sin\\theta_{2} - W_{1} & =\n",
    "0,\\tag*{6.6}\\\\ T_{1}\\cos\\theta_{1} - T_{2}\\cos\\theta_{2} & = 0,\\tag*{6.7}\\\\\n",
    "T_{2}\\sin\\theta_{2} + T_{3}\\sin\\theta_{3} - W_{2} & = 0,\\tag*{6.8}\\\\\n",
    "T_{2}\\cos\\theta_{2} - T_{3}\\cos\\theta_{3} & = 0.\\tag*{6.9}\\end{align}$$\n",
    "\n",
    " Here *W*<sub>*i*</sub> is the weight of mass *i* and *T*<sub>*i*</sub>\n",
    "is the tension in string *i*. Note that because we do not have a rigid\n",
    "structure, we cannot assume the equilibrium of torques."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 6.1.2  Algorithm: Multidimensional Searching<a id=\"6.1.2\"></a>\n",
    "\n",
    "[![image](Figs/RHLlectureMod4.png)](http://physics.oregonstate.edu/~rubin/Books/CPbook/eBook/Lectures/Modules/NdimSearch/NdimSearch.html)\n",
    "\n",
    "Equations (6.1)-(6.9) are nine simultaneous nonlinear equations. While\n",
    "linear equations can be solved directly, nonlinear equations cannot\n",
    "\\[[Press et al.(94)](BiblioLinked.html#press)\\]. You can use the computer to *search* for a\n",
    "solution by guessing, but there is no guarantee of finding one.\n",
    "\n",
    "Unfortunately, not everything in life is logical, as we need to use a\n",
    "search technique now that will not be covered until [Chapter 7,\n",
    "*Trial-and-Error Searching & Data Fitting*](CP07.ipynb). While what we\n",
    "do next is self explanatory, you may want to look at\n",
    "[Chapter 7](CP07.ipynb) now if you are not at all familiar with\n",
    "searching.\n",
    "\n",
    "We apply to our set the same Newton-Raphson algorithm as used to solve a\n",
    "single equation by renaming the nine unknown angles and tensions as the\n",
    "subscripted variable *y*<sub>*i*</sub> and placing the variables\n",
    "together as a vector:\n",
    "\n",
    "$$\\tag*{6.10}\n",
    "\\vec{y} =\\begin{pmatrix}\n",
    "x_{1} \\\\ x_{2} \\\\ x_{3} \\\\ x_{4} \\\\ x_{5} \\\\ x_6 \\\\ x_{7}\\\\ x_{8}\n",
    "\\\\ x_9\n",
    "\\end{pmatrix} =\\begin{pmatrix}\n",
    "\\sin\\theta_{1} \\\\\n",
    "\\sin\\theta_{2}\\\\\n",
    "\\sin\\theta_{3} \\\\\n",
    "\\cos\\theta_{1} \\\\\n",
    "\\cos\\theta_2 \\\\\n",
    "\\cos\\theta_3 \\\\\n",
    "T_1 \\\\ T_{2}\\\\ T_3\n",
    "\\end{pmatrix}.$$\n",
    "\n",
    "The nine equations to be solved are written in a general form with zeros\n",
    "on the right-hand sides and placed in a vector:\n",
    "\n",
    "$$\\begin{align}\n",
    "\\tag*{6.11}\n",
    "f_i (x_1, x_2,\\ldots, x_N) & = 0,\\quad i=1, N,\\\\\n",
    "\\vec{f}(\\vec{y}) =\\begin{pmatrix}\n",
    "f_{1}(\\vec{y})\\\\ f_{2}(\\vec{y}) \\\\ f_3(\\vec{y})\\\\ f_4(\\vec{y})\\\\ f_5(\\vec{y})\\\\\n",
    "f_6(\\vec{y})\\\\ f_7(\\vec{y})\\\\ f_8(\\vec{y})\\\\ f_{9}(\\vec{y})\n",
    "\\end{pmatrix}  & =\\begin{pmatrix}\n",
    "3 x_4 + 4 x_5 + 4x_6 -8 \\\\\n",
    " 3 x_1 + 4 x_2 - 4 x_3  \\\\\n",
    " x_7 x_1 - x_8 x_2  -10 \\\\\n",
    " x_7 x_4 - x_8 x_5   \\\\\n",
    " x_8 x_2 + x_9 x_3 -20 \\\\\n",
    "x_8 x_5 - x_9 x_6\\\\ x_1^2 + x_4^2 -1 \\\\ x_{2}^2 + x_5^2 - 1 \\\\ x_3^2 + x_6^2 -1\n",
    "\\end{pmatrix}\n",
    " = \\vec{0}.\\tag*{6.12}\\end{align}$$\n",
    "\n",
    "The solution to these equations requires a set of nine *x*<sub>*i*</sub> values\n",
    "that make all nine *f*<sub>*i*</sub>’s vanish simultaneously. Although these\n",
    "equations are not very complicated (the physics after all is elementary), the\n",
    "terms quadratic in *x* make them nonlinear, and this makes it hard or\n",
    "impossible to find an analytic solution. The search algorithm guesses a solution,\n",
    "expands the nonlinear equations into linear form, solves the resulting linear\n",
    "equations, and continues to improve the guesses based on how close the\n",
    "previous one was to making $\\vec{f} =0$. (We discuss search algorithms using\n",
    "this procedure in [Chapter 7, *Trial-and-Error Searching & Data\n",
    "Fitting*](CP07.ipynb).)\n",
    "\n",
    "Explicitly, let the approximate solution at any one stage be the set\n",
    "*x*<sub>*i*</sub> and let us assume that there is an (unknown) set of\n",
    "corrections *Δ**x*<sub>*i*</sub> for which\n",
    "\n",
    "$$\\tag*{6.13} f_i (x_1+\\Delta x_1, x_2+\\Delta x_2, \\ldots, x_9+\\Delta x_9)= 0,\n",
    "\\quad i=1, 9.$$\n",
    "\n",
    "We solve for the approximate *Δ**x*<sub>*i*</sub>’s by assuming that our\n",
    "previous solution is close enough to the actual one for two terms in the\n",
    "Taylor series to be accurate:\n",
    "\n",
    "$$\\begin{align}\n",
    "\\tag*{6.14}\n",
    "f_i(x_1+\\Delta x_1,\\ldots x_9+\\Delta x_9) \\simeq f_i (x_1,\\ldots x_9) +\n",
    "\\sum_{j=1}^9 \\frac{\\partial f_i} {\\partial x_j} \\Delta x_j &= 0\\\\\n",
    "  i &= 1,\\ldots 9.\\end{align}$$\n",
    "\n",
    "We now have a solvable set of nine linear equations in the nine unknowns\n",
    "*Δ**x*<sub>*i*</sub>, which we express as a single matrix equation\n",
    "\n",
    "$$\\begin{align} f_1 + {\\partial f_1/\\partial x_1} \\Delta x_1 + {\\partial f_1/\\partial\n",
    "x_2}\n",
    "\\Delta x_2 + \\cdots + {\\partial f_1 / \\partial x_9}  \\Delta x_9 & =  0,   \\\\\n",
    "f_2 + {\\partial f_2 / \\partial x_1} \\Delta x_1 + {\\partial f_2 /\n",
    "\\partial x_2}  \\Delta x_2 + \\cdots + {\\partial f_2 / \\partial x_9}  \\Delta x_9 & =  0,   \\\\\n",
    "\\ddots &  \\\\\n",
    " f_9 + {\\partial f_9 /\\partial x_1}  \\Delta x_1 + {\\partial f_9 /\n",
    "\\partial x_2}  \\Delta x_2 + \\cdots + {\\partial f_9 / \\partial\n",
    "x_9} \\Delta x_9 & = 0, \\\\\n",
    "\\begin{pmatrix}\n",
    "f_1 \\\\ f_2\\\\\n",
    "\\ddots\\\\\n",
    " f_9\n",
    " \\end{pmatrix} +\\begin{pmatrix}\n",
    " {\\partial f_1 / \\partial x_1} & {\\partial f_1 / \\partial x_2}\n",
    "& \\cdots & {\\partial f_1 / \\partial x_9} \\\\ {\\partial f_2 /\\partial x_1} & {\\partial\n",
    "f_2 / \\partial x_2} & \\cdots & {\\partial f_2 / \\partial x_9} \\\\\n",
    "\\ddots \\\\\n",
    "{\\partial f_9 / \\partial x_1} & {\\partial f_9 /\n",
    "\\partial x_2} & \\cdots & {\\partial f_9 / \\partial x_9}\n",
    "\\end{pmatrix}\n",
    "\\begin{pmatrix}\n",
    "\\Delta x_1 \\\\\n",
    "\\Delta x_2 \\\\\n",
    "\\ddots \\\\        \\Delta x_9\n",
    "\\end{pmatrix} & =         0.             \\tag*{6.15}\\end{align}$$\n",
    "\n",
    "Note now that the derivatives and the *f*’s are all evaluated at known values of\n",
    "the *x*<sub>i</sub>’s, so that only the vector of the *Δx*<sub>i</sub>\n",
    "values is unknown. We write this equation in matrix notation as[[xml]](http://physics.oregonstate.edu/~rubin/Books/CPbook/eBook/xml/6.16.xml) \n",
    "\n",
    "$$\\begin{align}\n",
    "\\tag*{6.16}\n",
    " \\vec{f}+ {F'}  \\vec{\\Delta x} = 0,\\quad &\\Rightarrow\\quad\n",
    "{F'}\\vec{\\Delta x} = - \\vec{f}, \\\\\n",
    " &\\vec{\\Delta x}  =\\begin{pmatrix}\n",
    " \\Delta x_1 \\\\\n",
    " \\Delta x_2 \\\\\n",
    " \\ddots\\\\\n",
    " \\Delta x_9\n",
    " \\end{pmatrix},\\quad \\vec{f} =\\begin{pmatrix}\n",
    " f_1 \\\\\n",
    " f_2\\\\\n",
    "\\ddots \\\\\n",
    "f_9\n",
    "\\end{pmatrix}, \\quad  &{F'} =\\begin{pmatrix}\n",
    "{\\partial f_1 / \\partial x_1} & \\cdots & {\\partial f_1 / \\partial x_9} \\\\ {\\partial f_2\n",
    "/ \\partial x_1} & \\cdots & {\\partial f_2 / \\partial x_9} \\\\\n",
    "\\ddots \\\\\n",
    "{\\partial f_9/\\partial x_1} & \\cdots & {\\partial f_9 /\\partial x_9}\n",
    "\\end{pmatrix}.\\end{align}$$\n",
    "\n",
    "Here we use arrows to emphasize the vector nature of the columns of\n",
    "*f*<sub>i</sub> and *Δx*<sub>i</sub> values, and call the matrix\n",
    "of the derivatives *F*′ (it is also sometimes called *J* because it is\n",
    "the *Jacobian* matrix).\n",
    "\n",
    "The equation ${F'}\\vec{\\Delta x} = - \\vec{f}$ is in the standard form for the\n",
    "solution of a linear equation (often written ${A}\\vec{x}=\\vec{b}$), where\n",
    "$\\vec{\\Delta x}$ is the vector of unknowns and . Matrix equations are solved\n",
    "using the techniques of linear algebra, and in the sections to follow we shall\n",
    "show how to do that. In a formal sense, the solution of (6.16) is obtained by\n",
    "multiplying both sides of the equation by the inverse of the *F*′ matrix:\n",
    "\n",
    "$$\\tag*{6.17}\n",
    "\\vec{\\Delta x} = - {F'}^{-1}\\vec{f},$$\n",
    "\n",
    "where the inverse must exist if there is to be a unique solution.\n",
    "Although we are dealing with matrices now, this solution is identical in\n",
    "form to that of the 1-D problem, *Δx* = −(1/*f*′)*f*. In fact, one of\n",
    "the reasons we use formal or abstract notation for matrices is to reveal\n",
    "the simplicity that lies within.\n",
    "\n",
    "As we indicate for the single-equation Newton-Raphson method in §7.3,\n",
    "even in a case such as this where we can deduce analytic expressions for\n",
    "the derivatives ∂*f*<sub>i</sub>/∂*x*<sub>j</sub>, there are\n",
    "9 × 9 = 81 such derivatives for this (small) problem, and entering them\n",
    "all would be both time-consuming and error-prone. In contrast,\n",
    "especially for more complicated problems, it is straightforward to\n",
    "program a forward-difference approximation for the derivatives,\n",
    "\n",
    "$$\\tag*{6.18}\n",
    "\\frac{\\partial f_i}{\\partial x_j} \\simeq \\frac{f_i(x_j + \\delta\n",
    "x_j) - f_i(x_j)}{\\delta x_j},$$\n",
    "\n",
    "where each individual *x*<sub>*j*</sub> is varied independently because\n",
    "these are partial derivatives and *δx*<sub>j</sub> are some\n",
    "arbitrary changes you input. While a central-difference approximation\n",
    "for the derivative would be more accurate, it would also require more\n",
    "evaluations of the *f*’s, and once we find a solution it does not matter\n",
    "how accurate our algorithm for the derivative was.\n",
    "\n",
    "As also discussed for the 1-D Newton-Raphson method (§7.3.1), the method\n",
    "can fail if the initial guess is not close enough to the zero of *f*\n",
    "(here all *N* of them) for the *f*’s to be approximated as linear. The\n",
    "*backtracking* technique may be applied here as well, in the present\n",
    "case, progressively decreasing the corrections *Δx*<sub>i</sub>\n",
    "until\n",
    "|*f*|<sup>2</sup> = |*f*<sub>1</sub>|<sup>2</sup> + |*f*<sub>2</sub>|<sup>2</sup> + ⋯ + |*f*<sub>*N*</sub>|<sup>2</sup>\n",
    "decreases."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 6.2  Why Matrix Computing?<a id=\"6.2\"></a>\n",
    "\n",
    "Physical systems are often modeled by systems of simultaneous equations\n",
    "written in matrix form. As the models are made more realistic, the\n",
    "matrices correspondingly become larger, and it becomes more important to\n",
    "use a good linear algebra library. Computers are unusually good with\n",
    "matrix manipulations because those manipulations typically involve the\n",
    "continued repetition of a small number of simple instructions, and\n",
    "algorithms exist to do this quite efficiently. Further speedup may be\n",
    "achieved by *tuning* the codes to the computer’s architecture, as\n",
    "discussed in [Chapter 11, *Applied HPC: Optimization, Tuning & GPU\n",
    "Programming*.](CP11.ipynb)\n",
    "\n",
    "Industrial-strength subroutines for matrix computing are found in\n",
    "well-established scientific libraries. These subroutines are usually an\n",
    "order of magnitude or more faster than the elementary methods found in\n",
    "linear algebra texts,\\[*Note:* Although we prize the book \\[[Press et\n",
    "al.(94)](BiblioLinked.html#press)\\] and what it has accomplished, we cannot recommend taking\n",
    "subroutines from it. They are neither optimized nor documented for easy,\n",
    "stand-alone use, whereas the subroutine libraries recommended in this\n",
    "chapter are.\\] are usually designed to minimize round-off error, and are\n",
    "often “robust,” that is, have a high chance of being successful for a\n",
    "broad class of problems. For these reasons we recommend that you *do not\n",
    "write your own matrix methods* but instead get them from a library. An\n",
    "additional value of library routines is that you can often run the same\n",
    "program either on a desktop machine or on a parallel supercomputer, with\n",
    "matrix routines automatically adapting to the local architecture. The\n",
    "thoughtful reader may be wondering when a matrix is “large” enough to\n",
    "require the use of a library routine. One rule of thumb is “if you have\n",
    "to wait for the answer”, and another rule is if the matrices you are\n",
    "using take up a good fraction of your computer’s random-access memory\n",
    "(RAM)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 6.3  Classes of Matrix Problems (Math) <a id=\"6.3\"></a>\n",
    "\n",
    "It helps to remember that the rules of mathematics apply even to the\n",
    "world’s most powerful computers. For example, you *should* encounter\n",
    "problems solving equations if you have more unknowns than equations, or\n",
    "if your equations are not linearly independent. But do not fret. While\n",
    "you cannot obtain a unique solution when there are not enough equations,\n",
    "you may still be able to map out a space of allowable solutions. At the\n",
    "other extreme, if you have more equations than unknowns, you have an\n",
    "*overdetermined* problem, which may not have a unique solution. An\n",
    "overdetermined problem is sometimes treated using data fitting\n",
    "techniques in which a solution to a sufficient set of equations is\n",
    "found, tested on the unused equations, and then improved if needed. Not\n",
    "surprisingly, this latter technique is known as the *linear\n",
    "least-squares method* (as in [Chapter 7, *Trial-and-Error Searching and\n",
    "Data Fitting*](CP07.ipynb)) because the technique minimizes the\n",
    "disagreement with the equations.\n",
    "\n",
    "The most basic matrix problem is a system of linear equations:\n",
    "\n",
    "$$\\tag*{6.19} {A} \\vec{x} = \\vec{b},$$\n",
    "\n",
    "where *A* is a known *N* × *N* matrix, $\\vec{x}$ is an unknown vector of\n",
    "length *N*, and $\\vec{b}$ is a known vector of length *N*. The obvious way to\n",
    "solve this equation is to determine the inverse of *A* and then form the solution\n",
    "by multiplying both sides of (6.19) by *A*<sup>−1</sup>:\n",
    "\n",
    "$$\\begin{align}\n",
    "\\tag*{6.20}\n",
    "\\vec{x} = {A}^{-1}   \\vec{b}.\\end{align}$$\n",
    "\n",
    "Both the direct solution of (6.19) and the determination of a matrix’s inverse\n",
    "are standards in a matrix subroutine library. A more efficient way to solve (6.19)\n",
    "is by Gaussian elimination or lower-upper (LU) decomposition. This yields the\n",
    "vector $\\vec{x}$ without explicitly calculating *A*<sup>−1</sup>. However,\n",
    "sometime you may want the inverse for other purposes, in which case (6.20) is\n",
    "preferred.\n",
    "\n",
    "If you have to solve the matrix equation\n",
    "\n",
    "$$\\tag*{6.21} {A} \\vec{x} = \\lambda\n",
    "\\vec{x},$$\n",
    "\n",
    "with $\\vec{x}$ an unknown vector and *λ* an unknown parameter, then the\n",
    "direct solution (6.20) will not be of much help because the matrix $\\vec{b}=\n",
    "\\lambda\n",
    "\\vec{x}$ contains the unknowns *λ* and $\\vec{x}$. Equation (6.21) is\n",
    "the *eigenvalue problem*. It is harder to solve than (6.19) because\n",
    "solutions exist for only certain, if any, values of *λ*. To find a\n",
    "solution, we use the identity matrix to rewrite (6.21) as\n",
    "\n",
    "$$\\tag*{6.22} [{A}-\\lambda {I}]\\vec{x} = 0.$$\n",
    "\n",
    "We see that multiplication of (6.22) by \\[*A* − *λI*\\]<sup>−1</sup>\n",
    "yields the *trivial solution*\n",
    "\n",
    "$$\\tag*{6.23}\n",
    "\\vec{x} = 0 \\quad (\\mbox{trivial solution}).$$\n",
    "\n",
    "While the trivial solution is a bona fide solution, it is nonetheless\n",
    "trivial. A more interesting solution requires the existence of a\n",
    "condition that forbids us from multiplying both sides of (6.22) by\n",
    "\\[*A* − *λI*\\]<sup>−1</sup>.That condition is the nonexistence of the\n",
    "inverse, and if you recall that Cramer’s rule for the inverse requires\n",
    "division by det\\[*A* − *λI*\\], it is clear that the inverse fails to\n",
    "exist (and in this way eigenvalues *do* exist) when\n",
    "\n",
    "$$\\tag*{6.24}\n",
    "\\det[A − λI] = 0.$$\n",
    "\n",
    "\n",
    "\n",
    "The traditional way to solve the eigenvalue problem (6.21) for both eigenvalues\n",
    "and eigenvectors is by *diagonalization*. This is equivalent to successive\n",
    "changes of basis vectors, each change leaving the eigenvalues unchanged while\n",
    "continually decreasing the values of the off-diagonal elements of *A*. The\n",
    "sequence of transformations is equivalent to continually operating on the\n",
    "original equation with the transformation matrix *U* until one is found for\n",
    "which *UAU*<sup>−1</sup> is diagonal:\n",
    "\n",
    " $$\\begin{align}\n",
    "\\tag*{6.25}\n",
    "{U A} ({U}^{-1}{U}) \\vec{x} & = \\lambda {U}\\vec{x} ,\\\\ ({U} {A} {U}^{-1})\n",
    "({U}\\vec{x}) & =\n",
    "\\lambda {U}\\vec{x},\\tag*{6.26}\\\\\n",
    "{U} {A} {U}^{-1} &=\\begin{pmatrix}\n",
    "\\lambda_1^{'} &  & \\cdots& 0\\\\\n",
    "0& \\lambda_2^{'} &\\cdots& 0\\\\ 0 & 0& \\lambda_3^{'} & \\cdots\\\\ 0 &\n",
    "\\cdots & & \\lambda_N^{'}\n",
    "\\end{pmatrix}.\\tag*{6.27}\\end{align}$$\n",
    "\n",
    " The diagonal values of *UAU*<sup>−1</sup> are the eigenvalues with\n",
    "eigenvectors\n",
    "\n",
    "$$\\tag*{6.28}\n",
    "\\vec{x}_{i} =  {U}^{-1}\n",
    "\\hat{e}_i,$$\n",
    "\n",
    "that is, the eigenvectors are the columns of the matrix\n",
    "*U*<sup>−1</sup>. A number of routines of this type are found in\n",
    "subroutine libraries.\n",
    "\n",
    "### 6.3.1  Practical Matrix Computing<a id=\"6.3.1\"></a>\n",
    "\n",
    "Many scientific programming bugs arise from the improper use of\n",
    "arrays.\\[*Note:* Even a vector *V*(*N*) is called an array, albeit a 1-D\n",
    "one.\\] This may be as a result of the extensive use of matrices in\n",
    "scientific computing or to the complexity of keeping track of indices\n",
    "and dimensions.In any case, here are some rules of thumb to observe.\n",
    "\n",
    "**Computers are finite:** Unless you are careful, your matrices may so\n",
    "much memory that your computation will slow down significantly,\n",
    "especially if it starts to use virtual memory. As a case in point, let’s\n",
    "say that you store data in a 4-D array with each index having a\n",
    "<span>*physical dimension*</span> of 100: `A[100] [100] [100] [100]`.\n",
    "This array of (100)<sup>4</sup> 64-byte words occupies ≃1*GB* of\n",
    "memory.\n",
    "\n",
    "**Processing time:** Matrix operations such as inversion require on the\n",
    "order of *N*<sup>3</sup> steps for a square matrix of dimension *N*.\n",
    "Therefore, doubling the dimensions of a 2-D square matrix (as happens\n",
    "when the number of integration steps is doubled) leads to an *eightfold*\n",
    "increase in processing time.\n",
    "\n",
    "**Paging:** Many operating systems have [*virtual\n",
    "memory*](http://www.science.oregonstate.edu/~rubin/Books/CPbook/eBook/GlossarySound/virtualmemory.wav)\n",
    "in which disk space is used when a program runs out of RAM (see\n",
    "[Chapter 10, *High-Performance Hardware & Parallel\n",
    "Computers*](CP10.ipynb), for a discussion of how computers arrange\n",
    "memory). This is a slow process that requires writing a full *page* of\n",
    "words to the disk. If your program is near the memory limit at which\n",
    "paging occurs, even a slight increase in a matrix’s dimension may lead\n",
    "to an order-of-magnitude increase in execution time.\n",
    "\n",
    "**Matrix storage:** While we think of matrices as multidimensional\n",
    "blocks of stored numbers, the computer stores them as linear strings.\n",
    "For instance, a matrix `a[3,3]` in Python is stored in [*row-major\n",
    "order*](http://www.science.oregonstate.edu/~rubin/Books/CPbook/eBook/GlossarySound/rowmajor.wav)\n",
    "(Figure 6.3 A):\n",
    "\n",
    "$$\\tag*{6.29} a_{0,0}\\ \\ a_{0,1} \\ \\ a_{0,2} \\ \\ a_{1,0}\\ \\ a_{1,1}\\ \\ a_{1,2}\\\n",
    "\\ a_{2,0} \\ \\ a_{2,1}\\ \\ a_{2,2}\\ \\ \\ldots ,$$\n",
    "\n",
    "while in Fortran, starting subscripts at 0, it is stored in\n",
    "[*column-major\n",
    "order*](http://www.science.oregonstate.edu/~rubin/Books/CPbook/eBook/GlossarySound/columnmajor.wav)\n",
    "(Figure 6.3 B):\n",
    "\n",
    "$$\\tag*{6.30}\n",
    " a_{0,1}\\ \\ a_{1,0}\\ \\  a_{2,0} \\  \\ a_{0,1} \\  \\ a_{1,1}\\  \\ a_{2,1}\\ \\ a_{0,2}\\  \\ a_{1,2}\\\n",
    "\\ a_{2,2}\\ \\  \\ldots.$$\n",
    "\n",
    "It is important to keep this linear storage scheme in mind in order to\n",
    "write proper code and to permit the mixing of Python and Fortran\n",
    "programs.\n",
    "\n",
    "![image](Figs/Fig6_3.png) \n",
    "\n",
    "**Figure 6.3** *A, Left:* Row-major order used\n",
    "for matrix storage in Python, C and Java. *B, Right:* Column-major order\n",
    "used for matrix storage in Fortran. The table on the bottom shows how\n",
    "successive matrix elements are actually stored in a linear fashion in\n",
    "memory.\n",
    "\n",
    "When dealing with matrices, you have to balance the clarity of the\n",
    "operations being performed against the efficiency with which the\n",
    "computer performs them. For example, having one matrix with many indices\n",
    "such as `V[L,Nre,Nspin,k,kp,Z,A]` may be neat packaging, but it may\n",
    "require the computer to jump through large blocks of memory to get to\n",
    "the particular values needed (large *strides*) as you vary `k`, `kp`,\n",
    "and `Nre`. The solution would be to have several matrices such as `V1`,\n",
    "`V2[Nre,Nspin,k,kp,Z,A]`, and `V3[Nre,Nspin,k,kp,Z,A]`.\n",
    "\n",
    "**Subscript 0:** It is standard in Python, C and Java to have array\n",
    "indices begin with the value 0. While this is now permitted in Fortran,\n",
    "the standard in Fortran and in most mathematical equations has been to\n",
    "start indices at 1. On that account, in addition to the different\n",
    "locations in memory as a result of row-major and column-major ordering,\n",
    "the same matrix element may be referenced differently in the different\n",
    "languages:\n",
    "\n",
    "|Location | Python/C Element | Fortran Element |\n",
    "|- - -|:- - -:|:- - -:|\n",
    "|Lowest| `a[0,0]` | `a(1,1)` |\n",
    "||`a[0,1]` | `a(2,1)`|\n",
    "||`a[1,0]` | `a(3,1)`|\n",
    "||`a[1,1]` | `a(1,2)`|\n",
    "||`a[2,0]` | `a(2,2)`|\n",
    "|Highest| `a[2,1]` |`a(3,2)`|\n",
    "**Tests:** *Always* test a library routine on a small problem whose\n",
    "answer you know (such as the exercises in §6.6). Then you’ll know if you\n",
    "are supplying it with the right arguments and if you have all the links\n",
    "working."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 6.4  Python Lists as Arrays<a id=\"6.4\"></a>\n",
    "\n",
    "A *list* is Python’s built-in sequence of numbers or objects. Although\n",
    "called a “list” it is similar to what other computer languages call an\n",
    "“array”. It may be easier for you to think of a Python list as a\n",
    "container that holds a bunch of items in a definite order. (Soon we will\n",
    "describe the higher-level `array` data types available with the *NumPy*\n",
    "package.) In this section we review some of Python’s native *list*\n",
    "features.\n",
    "\n",
    "Python interprets a sequence of ordered items, $L = l_0,l_1, \\ldots,\n",
    "l_{N-1}$, as a *list* and represents it with a single symbol `L`:\n",
    "\n",
    "\n",
    "\n",
    "     >>> L = [1, 2, 3]                   # Create list\n",
    "     >>> L[0]                                # Print element 0 (first)\n",
    "    1                                             # Python output\n",
    "    >>> L                                      # Print entire list\n",
    "     [1, 2, 3]                                # Output\n",
    "    >>> L[0]= 5                           # Change element 0\n",
    "    >>> L\n",
    "    [5, 2, 3]\n",
    "    >>> len(L)                           # Length of list\n",
    "    3\n",
    "    >>> for items in L: print items      # for loop over items\n",
    "    5\n",
    "    2\n",
    "    3\n",
    "\n",
    "We observe that square brackets with comma separators \\[1, 2, 3\\] are\n",
    "used for lists, and that a square bracket is also used to indicate the\n",
    "index for a list item, as in line 2 (`L[0]`). Lists contain sequences of\n",
    "arbitrary objects that are *mutable* or changeable. As we see in line 7\n",
    "in the `L` command, an entire list can be referenced as a single object,\n",
    "in this case to obtain its printout.\n",
    "\n",
    "Python also has a built-in type of list know as a *tuple* whose elements\n",
    "are not mutable. Tuples are indicated by round parenthesis (.., .., .),\n",
    "with individual elements still referenced by square brackets:\n",
    "\n",
    "    >>> T = (1, 2, 3, 4)                 # Create tuple\n",
    "    >>> T[3]                                  # Print element 3\n",
    "    4\n",
    "    >>> T\n",
    "    (1, 2, 3, 4)                             # Print entire tuple\n",
    "    >>> T[0] = 5                          # Attemp to change element 0\n",
    "    Traceback (most recent call last):\n",
    "        T[0] = 5\n",
    "    TypeError: 'tuple' object does not support item assignment\n",
    "\n",
    "Note that the error message that arises when we try to change an element\n",
    "of a tuple.\n",
    "\n",
    "Most languages require you to specify the size of an\n",
    "[*array*](http://www.science.oregonstate.edu/~rubin/Books/CPbook/eBook/GlossarySound/array.wav)\n",
    "before you can start storing objects in it. In contrast, Python lists\n",
    "are *dynamic*, which means that their sizes adjust as needed. In\n",
    "addition, while a list is essentially one dimensional because it is a\n",
    "sequence, a compound list can be created in Python by having the\n",
    "individual elements themselves as lists:\n",
    "\n",
    "    >>> L = [[1,2], [3,4], [5,6]]        # A list of lists\n",
    "    >>> L\n",
    "    [[1, 2], [3, 4], [5, 6]]\n",
    "    >>> L[0]                                        # The first element\n",
    "    [1, 2]\n",
    "\n",
    "Python can perform a large number of operations on lists, for example,\n",
    "\n",
    "|**Operation** | **Effect**|**Operation** | **Effect**|\n",
    "|- - -|- - -|- - -|- - -|\n",
    "|L = \\[1, 2, 3, 4\\] | Form list |L1 + L2 | Concatenate lists|\n",
    "|L\\[i\\] | i<sup>th</sup> element | len(L) | Length of list L|\n",
    "|i in L | True if i in L | L\\[i:j\\] | Slice from i to j|\n",
    "|for i in L | Iteration index | L.append(x) | Append x to end of L|\n",
    "|L.count(x) | Number of x’s in L|L.index(x) | Location of 1st x in L|\n",
    "|L.remove(x) | Remove 1st x in L|L.reverse()| Reverse elements in L|\n",
    "|L.sort()| Order elements in L|||"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 6.5  Numerical Python (NumPy) Arrays<a id=\"6.5\"></a>\n",
    "\n",
    "Although basic Python does have an `array` data type, it is rather\n",
    "limited and we suggest using NumPy arrays, which will converts Python\n",
    "lists into arrays. Because *we often use NumPy’s `array` command* to\n",
    "form computer representations of vectors and matrices, you will have to\n",
    "import *NumPy* in your programs (or <span>Visual</span>, which includes\n",
    "<span>NumPy</span>) to use those examples. For instance, here we show\n",
    "the results of running our program `Matrix.py` from a shell:\n",
    "\n",
    "    >>> from visual import *                      # Load Visual package\n",
    "    >>> vector1 = array([1, 2, 3, 4, 5])     # Fill 1-D array\n",
    "    >>> print('vector1 =',vector1)            # Print array (parens if Python 3)\n",
    "    vector1 =  [1 2 3 4 5]                            # Output\n",
    "    >>> vector2 = vector1 + vector1        # Add 2 vectors\n",
    "    >>> print('vector2=',vector2)             # Print vector2\n",
    "    vector2= [ 2  4  6  8 10]                       # Output\n",
    "    >>> vector2 = 3 * vector1                    # Mult array by scalar\n",
    "    >>> print ('3 * vector1 = ', vector2)    # Print vector\n",
    "    3 * vector1 =  [ 3  6  9 12 15]              # Output\n",
    "    >>> matrix1 = array(([0,1],[1,3]))       # An array of arrays\n",
    "    >>> print(matrix1)                                # Print matrix1\n",
    "    [[0 1]\n",
    "     [1 3]]\n",
    "    >>> print ('vector1.shape= ',vector1.shape)\n",
    "    vector1.shape =  (5)\n",
    "    >>> print  (matrix1 * matrix1)           # Matrix multiply\n",
    "     [[0 1]\n",
    "      [1 9]]\n",
    "\n",
    "We see here that we have initialized an array object, have added two 1-D array\n",
    "objects together and have printed out the result. Likewise we see that\n",
    "multiplying an array by a constant does in fact multiply each element by that\n",
    "constant (line 8). We then construct a “matrix” as a 1-D array of two 1-D arrays,\n",
    "and when we print it out, we note that it does indeed look like a matrix.\n",
    "However, when we multiply this matrix by itself, the result is not the\n",
    "$\\left[\\begin{array}{cc} 1&3\\\\3&10\\end{array}\\right]$ that one normally\n",
    "expects from matrix multiplication. So if you need actual mathematical\n",
    "matrices, then you need to use NumPy!\n",
    "\n",
    "Now we give some examples of the use of NumPy, but do refer the reader\n",
    "to the NumPy Tutorial \\[[NumPyTut(12)](BiblioLinked.html#numpytut)\\] and to the articles in *Computing\n",
    "in Science & Engineering* \\[[CiSE(15)](BiblioLinked.html#CiSE)\\] for more information. To start,\n",
    "we note that a NumPy array can hold up to 32 dimensions (32 indices),\n",
    "but each element must be of the same type (a *uniform* array). The\n",
    "elements are not restricted to just floating-point numbers or integers,\n",
    "but can be any object, as long as all elements are of this same type.\n",
    "(Compound objects may be useful, for example, for storing parts of data\n",
    "sets.) There are various ways to create arrays, with square brackets\n",
    "\\[…\\] used for indexing in all cases. First we start with a Python list\n",
    "(tupples work as well) and create an array from it:\n",
    "\n",
    "    >>> from numpy  import *\n",
    "    >>> a = array( [1, 2, 3, 4] )          # Array from a list\n",
    "    >>> a                                               # Check with print\n",
    "    array([1, 2, 3, 4])\n",
    "\n",
    "Notice that it is essential to have the square brackets within the round\n",
    "brackets because the square brackets produce the list object while the\n",
    "round brackets indicate a function argument. Note too that because the\n",
    "data in our original list were all integers, the created array is a 32\n",
    "bit integer data type, which we can check by affixing the `dtype`\n",
    "method:\n",
    "\n",
    "    >>> a.dtype\n",
    "    dtype('int32')\n",
    "\n",
    "If we had started with floating-point numbers, or a mix of floats and\n",
    "ints, we would have ended up with floating-point arrays:\n",
    "\n",
    "    >>> b = array([1.2, 2.3, 3.4])\n",
    "    >>> b\n",
    "    array([ 1.2,  2.3,  3.4])\n",
    "    >>> b.dtype\n",
    "    dtype('float64')\n",
    "\n",
    "When describing NumPy arrays, the number of “dimensions”, `ndim`, means\n",
    "the number of indices, which as we said can be as high as 32. What might\n",
    "be called the “size” or “dimensions” of a matrix in mathematics is\n",
    "called the *shape* of a NumPy array. Furthermore, NumPy does have a\n",
    "`size` method that returns the total number of elements. Because\n",
    "Python’s lists and tupples are all one dimensional, if we want an array\n",
    "of a particular shape, we can attain that by affixing the `reshape`\n",
    "method when we create the array. Where Python has a `range` function to\n",
    "generate a sequence of numbers, NumPy has an `arange` function that\n",
    "creates an array, rather than a list. Here we use it and then reshape\n",
    "the 1D array into a 3 × 4 array:\n",
    "\n",
    "    >>> import numpy as np\n",
    "    >>> np.arange(12)                    # List of 12 ints in 1D array\n",
    "    array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11])\n",
    "    >>> np.arange(12).reshape((3,4))     # Create, shape to 3x4 array\n",
    "    array([[ 0,  1,  2,  3],\n",
    "           [ 4,  5,  6,  7],\n",
    "           [ 8,  9, 10, 11]])\n",
    "    >>> a = np.arange(12).reshape((3,4))         # Give array a name\n",
    "    >>> a\n",
    "    array([[ 0,  1,  2,  3],\n",
    "           [ 4,  5,  6,  7],\n",
    "           [ 8,  9, 10, 11]])\n",
    "    >>> a.shape                                             # Shape = ?\n",
    "    (3L, 4L)\n",
    "    >>> a.ndim                                             # Dimension?\n",
    "    2\n",
    "    >>> a.size                        # Size of a (number of elements)?\n",
    "    12\n",
    "\n",
    "Note that here we imported NumPy as the object `np`, and then affixed\n",
    "the `arange` and `reshape` methods to this object. We then checked the\n",
    "shape of `a`, and found it to have three rows and four columns of long\n",
    "integers (Python 3 may just say `ints`). Note too, as we see on line 9,\n",
    "NumPy uses parentheses () to indicate the shape of an array, and so\n",
    "`(3L,4L)` indicates an array with 3 rows and 4 columns of long ints.\n",
    "\n",
    "Now that we have shapes on our minds, we should note that NumPy offers a\n",
    "number of ways to change shapes. For example we can transpose an array\n",
    "with the `.T` method, or `reshape` into a vector:\n",
    "\n",
    "    >>> from numpy import *\n",
    "    >>> a = arange(12).reshape((3,4))                   # Give array a name\n",
    "    >>> a\n",
    "    array([[ 0,  1,  2,  3],\n",
    "           [ 4,  5,  6,  7],\n",
    "           [ 8,  9, 10, 11]])\n",
    "    >>> a.T                                                                    # Transpose\n",
    "    array([[ 0,  4,  8],\n",
    "           [ 1,  5,  9],\n",
    "           [ 2,  6, 10],\n",
    "           [ 3,  7, 11]])\n",
    "    >>> b = a.reshape( (1,12) )                              # Form vector length 12\n",
    "    >>> b\n",
    "    array([[ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11]])\n",
    "\n",
    "And again, `(1,12)` indicates an array with one row and 12 columns. Yet\n",
    "another handy way to take a matrix and extract just what you want from\n",
    "it is to use Python’s *slice* operator `start:stop:step:` to take a\n",
    "slice out of an array:\n",
    "\n",
    "    >>> a\n",
    "    array([[ 0,  1,  2,  3],\n",
    "           [ 4,  5,  6,  7],\n",
    "           [ 8,  9, 10, 11]])\n",
    "    >>> a[:2, :]                                         # First 2 rows\n",
    "    array([[0, 1, 2, 3],\n",
    "           [4, 5, 6, 7]])\n",
    "    >>> a[:,1:3]                                        # Columns 1-3\n",
    "    array([[ 1,  2],\n",
    "           [ 5,  6],\n",
    "           [ 9, 10]])\n",
    "\n",
    "We note here that Python indices start counting from 0, and so 1:3 means\n",
    "indices 0, 1, 2 (without the 3). As we discuss in [Chapter 11, *Applied\n",
    "HPC: Optimization, Tuning & GPU Programming*](CP11.ipynb), slicing can\n",
    "be very useful in speeding up programs by picking out and placing in\n",
    "memory just the specific data elements from a large data set that need\n",
    "to be processed. This avoids the time-consuming jumping through large\n",
    "segments of memory as well as excessive reading from disk.\n",
    "\n",
    "Finally, we remind you that while all elements in a NumPy array must be\n",
    "of the same data type, that data type can be compound. For example, an\n",
    "array of arrays:\n",
    "\n",
    "    >>> from numpy import *\n",
    "    >>> M = array( [ (10, 20), (30,40), (50, 60) ] )       # Array of 3 arrays\n",
    "    >>> M\n",
    "    array([[10, 20],\n",
    "           [30, 40],\n",
    "           [50, 60]])\n",
    "    >>> M.shape\n",
    "    (3L, 2L)\n",
    "    >>> M.size\n",
    "    6\n",
    "    >>> M.dtype\n",
    "    dtype('int32')\n",
    "\n",
    "Furthermore, an array can be composed of complex numbers by specifying\n",
    "the `complex` data type as an option on the `array` command. NumPy then\n",
    "uses the *j* symbol for the imaginary number *i*:\n",
    "\n",
    "    >>> c = array( [ [1,complex(2,2)], [complex(3,2),4] ], dtype=complex )\n",
    "    >>> c\n",
    "    array([[ 1.+0.j,  2.+2.j],\n",
    "           [ 3.+2.j,  4.+0.j]])\n",
    "\n",
    "In the next section we discuss using true mathematical matrices with\n",
    "NumPy, which is one use of an array object. Here we note that if you\n",
    "wanted the familiar matrix product from two arrays, you would use the\n",
    "`dot` function, whereas `*` is used for an element-by-element product:\n",
    "\n",
    "    >>> matrix1= array( [[0,1], [1,3]])\n",
    "    >>> matrix1\n",
    "    array([[0, 1],\n",
    "           [1, 3]])\n",
    "    >>> print ( dot(matrix1,matrix1) )          # Matrix or dot product\n",
    "    [[ 1  3]\n",
    "     [ 3 10]]\n",
    "    >>> print (matrix1 * matrix1)               # Element-by-element product\n",
    "    [[0 1]\n",
    "     [1 9]]\n",
    "\n",
    "NumPy is actually optimized to work well with arrays, and in part this\n",
    "is because arrays are handled and processed much as if they were simple,\n",
    "scalar variables\\[*Note:* We thank Bruce Sherwood for helpful comments\n",
    "on these points.\\]. For example, here is another example of *slicing*, a\n",
    "technique that is also used in ordinary Python with lists and tuples, in\n",
    "which two indices separated by a colon indicate a range:\n",
    "\n",
    "    from visual import *\n",
    "    stuff = zeros(10, float)\n",
    "    t = arange(4)\n",
    "    stuff[3:7] = sqrt(t+1)\n",
    "\n",
    "Here we start by creating the NumPy array `stuff` of floats, all of\n",
    "whose 10 elements are initialized to zero. Then we create the array `t`\n",
    "containing the four elements \\[0, 1, 2, 3\\] by assigning 4 variables\n",
    "uniformly in the range 0-4 (the “a” in `arange` creates floating-point\n",
    "variables, `range` creates integers). Next we use a slice to assign\n",
    "\\[sqrt(0+1), sqrt(1+1), sqrt(2+1), sqrt(3+1)\\] = \\[1, 1.414, 1.732, 2\\]\n",
    "to the middle elements of the `stuff` array. Note that the NumPy version\n",
    "of the `sqrt` function, one of many universal function (*ufunctions*)\n",
    "supported by NumPy\\[*Note:* A ufunction is a function that operates on\n",
    "N-D arrays in an element-by-element fashion, supporting array\n",
    "broadcasting, type casting, and several other standard features. In\n",
    "other words, a ufunc is a vectorized wrapper for a function that takes a\n",
    "fixed number of scalar inputs and produces a fixed number of scalar\n",
    "outputs.\\], has the amazing property of automatically outputting an\n",
    "array whose length is that of its argument, in this case, the array `t`.\n",
    "In general, major power in NumPy comes from its *broadcasting*\n",
    "operation, an operation in which values are assigned to multiple\n",
    "elements via a single assignment statement. Broadcasting permits Python\n",
    "to *vectorize* array operations, which means that the same operation can\n",
    "be performed on different array elements in parallel (or nearly so).\n",
    "Broadcasting also speeds up processing because array operations occur in\n",
    "C instead of Python, and with a minimum of array copies being made. Here\n",
    "is a sample of broadcasting:\n",
    "\n",
    "    w = zeros(100, float)\n",
    "    w = 23.7\n",
    "\n",
    "The first line creates the NumPy array `w`, and the second line\n",
    "“broadcasts” the value 23.7 to all elements in the array. There are many\n",
    "possible array operations in NumPy and various rules pertaining to them;\n",
    "we recommend that the serious user explore the extensive NumPy\n",
    "documentation for additional information."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 6.5.1  NumPy’s linalg Package<a id=\"6.5.1\"></a>\n",
    "\n",
    "The array objects of NumPy and Visual are not the same as mathematical\n",
    "matrices, although an array can be used to represent a matrix.\n",
    "Fortunately, there is the `LinearAlgebra` package that treats 2-D arrays\n",
    "(a 1-D array of 1-D arrays) as mathematical matrices, and also provides\n",
    "a simple interface to the powerful LAPACK linear algebra library\n",
    "\\[[Anderson et al.(13)](BiblioLinked.html#folding)\\]. As we keep saying, there is much to be gained\n",
    "in speed and reliability from using these libraries rather than writing\n",
    "your own matrix routines.\n",
    "\n",
    "Our first example from linear algebra is the standard matrix equation\n",
    "\n",
    "$$\\tag*{6.31}\n",
    " {A}\\vec{x} = \\vec{b},$$\n",
    "\n",
    "where we have used a bold character to represent a matrix and the vector sign\n",
    "to represent a 1D matrix (a vector). Equation (6.31) describes a set of linear\n",
    "equations with $\\vec{x}$ an unknown vector and *A* a known matrix. Now we\n",
    "take *A* to be 3 × 3, $\\vec{b}$ to be 3 × 1, and let the program figure out that\n",
    "$\\vec{x}$ must be a 3 × 1 vector\\[*Note:* Don’t be bothered by the fact that\n",
    "although we think as these vectors as 3 × 1, they sometimes get printed out as\n",
    "1 × 3; think of all the trees that get saved!\\]. We start by importing all the\n",
    "packages, by inputting a matrix and a vector, and then by printing out *A* and\n",
    "$\\vec{x}$:\n",
    "\n",
    "    >>> from numpy import *\n",
    "    >>> from numpy.linalg import*\n",
    "    >>> A = array( [ [1,2,3], [22,32,42], [55,66,100] ] )  # Array of arrays\n",
    "    >>> print ('A =', A)\n",
    "    A = [[  1   2   3]\n",
    "         [ 22  32  42]\n",
    "         [ 55  66 100]]\n",
    "    >>> b = array([1,2,3])\n",
    "    >>> print ('b =', b)\n",
    "    b = [1 2 3]\n",
    "\n",
    "Because we have the matrices **A** and $\\vec{b}$, we can go ahead and solve\n",
    "$ {A}\\vec{x} = \\vec{b}$ using NumPy’s `solve` command, and then test how close\n",
    "$ {A}\\vec{x} - \\vec{b}$ is to a zero vector:\n",
    "\n",
    "    >>> from numpy.linalg import solve\n",
    "    >>> x = solve(A, b)                                     # Finds solution\n",
    "    >>> print ('x =', x)\n",
    "    x = [ -1.4057971  -0.1884058   0.92753623]              # The solution\n",
    "    >>> print ('Residual =',  dot(A, x) - b)                # LHS-RHS\n",
    "\n",
    "    Residual = [4.44089210e-16   0.00000000e+00  -3.55271368e-15]\n",
    "\n",
    "This is really quite impressive. We have solved the entire set of linear\n",
    "equations (by elimination) with just the single command `solve`,\n",
    "performed a matrix multiplication with the single command `dot`, did a\n",
    "matrix subtraction with the usual - operator, and are left with a\n",
    "residual essentially equal to machine precision.\n",
    "\n",
    "Although there are more efficient numerical approaches, a direct way to\n",
    "solve\n",
    "\n",
    "$$\\tag*{6.32}\n",
    " {A}\\vec{x} = \\vec{b}$$\n",
    "\n",
    "is to calculate the inverse *A*<sup>−1</sup>, and then multiply both\n",
    "sides of the equation by the inverse, yielding\n",
    "\n",
    "$$\\tag*{6.33}\n",
    "  \\vec{x} =  y{A}^{-1}\\vec{b}$$\n",
    "\n",
    "    >>> from numpy.linalg import inv\n",
    "    >>>  dot(inv(A), A)                               # Test inverse\n",
    "\n",
    "    array([[  1.00000000e+00,  -1.33226763e-15,  -1.77635684e-15],\n",
    "           [  8.88178420e-16,   1.00000000e+00,   0.00000000e+00],\n",
    "           [ -4.44089210e-16,   4.44089210e-16,   1.00000000e+00]])\n",
    "    >>> print ('x =', multiply(inv(A), b))\n",
    "    x = [-1.4057971  -0.1884058   0.92753623]             # Solution\n",
    "    >>> print ('Residual =',  dot(A, x) - b)\n",
    "\n",
    "    Residual = [  4.44089210e-16   0.00000000e+00  -3.55271368e-15]\n",
    "\n",
    "Here we first tested that `inv(A)` is in fact the inverse of `A` by\n",
    "seeing if `A` times `inv(A)` equals the identity matrix. Then we used\n",
    "the inverse to solve the matrix equation directly, and got the same\n",
    "answer as before (some error at the level of machine precision is just\n",
    "fine).\n",
    "\n",
    "Our second example occurs in the solution for the principal-axes system\n",
    "of a cube, and requires us to find a coordinate system in which the\n",
    "inertia tensor is diagonal. This entails solving the eigenvalue problem,\n",
    "\n",
    "$$\\tag*{6.34}\n",
    "  \\vec{\\omega}   = \\lambda \\vec{\\omega},$$\n",
    "\n",
    "where *I* is the inertia matrix (tensor), $\\vec{\\omega}$ is an unknown\n",
    "eigenvector, and *λ* is an unknown eigenvalue. The program `Eigen.py` solves\n",
    "for the eigenvalues and vectors, and shows how easy it is to deal with matrices.\n",
    "Here it is in an abbreviated interpretive mode:\n",
    "\n",
    "    >>> from numpy import*\n",
    "    >>> from numpy.linalg import eig\n",
    "    >>> I = array( [[2./3,-1./4], [-1./4,2./3]] )\n",
    "    >>> print('I =\\n', I)\n",
    "    I =\n",
    "     [[ 0.66666667 -0.25      ]\n",
    "     [-0.25        0.66666667]]\n",
    "    >>> Es, evectors = eig(A)                   # Solves eigenvalue problem\n",
    "    >>> print('Eigenvalues =', Es, '\\n Eigenvector Matrix =\\n', evectors)\n",
    "    Eigenvalues =   [ 0.91666667  0.41666667]\n",
    "     Eigenvector Matrix =\n",
    "     [[ 0.70710678  0.70710678]\n",
    "     [-0.70710678  0.70710678]]\n",
    "    >>> Vec = array([ evectors[0, 0], evectors[1, 0] ] )\n",
    "    >>> LHS = dot(I, Vec)                        # Matrix x vector\n",
    "    >>> RHS = Es[0]*Vec                          # Scalar mult\n",
    "    >>> print('LHS - RHS =', LHS-RHS)            # Test for zero\n",
    "    LHS - RHS = [  1.11022302e-16  -1.11022302e-16]\n",
    "\n",
    "We see here how, after we set up the array `I` on line 3, we then solve\n",
    "for its eigenvalues and eigenvectors with the single statement\n",
    "`Es, evectors = eig(I)` on line 8. We then extract the first eigenvector\n",
    "on line 14 and use it, along with the first eigenvalue, to check that\n",
    "(6.34) is in fact satisfied to machine precision.\n",
    "\n",
    "Well, we think by now you have some idea of the use of NumPy. In\n",
    "Table 6.1 we indicate some more of what’s available.\n",
    "\n",
    "**Table 6.1** The operators of NumPy and their effects.\n",
    "\n",
    "|*Operator* | *Effect*|*Operator*|*Effect*| \n",
    "|- - -|- - - |- - -|- - - | \n",
    "|dot(a,b\\[,out\\]) | Dot product arrays | vdot(a, b) | Dot product | \n",
    "|inner(a, b) | Inner product arrays |outer(a, b) | Outer product | \n",
    "|tensordot(a, b) | Tensor dot product | einsum( ) |Einstein sum | \n",
    "|linalg.matrix_power(M, n) | Matrix to power n | kron(a, b) | Kronecker product| \n",
    "|linalg.cholesky(a) | Cholesky decomp |linalg.qr(a)| QR factorization | \n",
    "|linalg.svd(a ) | Singular val decomp |linalg.eig(a) | Eigenproblem | \n",
    "|linalg.eigh(a) | Hermitian eigen |linalg.eigvals(a)| General eigen |\n",
    "|linalg.eigvalsh(a) | Hermitian eigenvals |linalg.norm(x) |Matrix norm | \n",
    "|linalg.cond(x) | Condition number |linalg.det(a) | Determinant | \n",
    "|linalg.slogdet(a) | Sign & log(det) |trace(a) | Diagnol sum | \n",
    "|linalg.solve(a, b) | Solve equation |linalg.tensorsolve(a, b) | Solve a x = b | \n",
    "|linalg.lstsq(a, b) |Least-squares solve |linalg.inv(a) | Inverse | \n",
    "|linalg.pinv(a) | Penrose inverse|linalg.tensorinv(a)| Inverse N-D array |"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 6.6  Exercise: Testing Matrix Programs <a id=\"6.6\"></a>\n",
    "\n",
    "Before you direct the computer to go off crunching numbers on a million\n",
    "elements of some matrix, it’s a good idea to try out your procedures on\n",
    "a small matrix, especially one for which you know the right answer. In\n",
    "this way it will take you only a short time to realize how hard it is to\n",
    "get the calling procedure perfect! Here are some exercises.\n",
    "\n",
    "1.  Find the numerical inverse of \n",
    "\n",
    "       ${A} =\\begin{pmatrix}\n",
    "            +4&-2&+1\\\\\n",
    "            +3&+6&-4\\\\\n",
    "            +2&+1&+8\\end{pmatrix}.$\n",
    "\n",
    "    -   As a general check, applicable even if you do not know the\n",
    "        analytic answer, check your inverse in both directions; that is,\n",
    "        check that *AA*<sup>−1</sup> = *A*<sup>−1</sup>*A* = *I*, and\n",
    "        note the number of decimal places to which this is true. This\n",
    "        also gives you some idea of the precision of your calculation.\n",
    "\n",
    "    -   Determine the number of decimal places of agreement there is\n",
    "        between your numerical inverse and the analytic result:\n",
    "        \n",
    "        $ {A}^{-1} = \\dfrac{1}{263}\\begin{pmatrix}\n",
    "                +52&+17&+2\\\\\n",
    "                -32&+30&+19\\\\\n",
    "                -9&-8&+30 \\end{pmatrix}.$ \n",
    "                \n",
    "        Is this similar to the\n",
    "        error in *AA*<sup>−1</sup>?\n",
    "\n",
    "2.  Consider the same matrix *A* as before, here being used to describe\n",
    "    three simultaneous linear equations, ${A}\\vec{x} =\n",
    "    \\vec{b}$, or explicitly,\n",
    "\n",
    "    $$\\tag*{6.35}\n",
    "    \\begin{pmatrix}\n",
    "    a_{00} & a_{01}& a_{02}\\\\\n",
    "    a_{10} & a_{11}& a_{12}\\\\\n",
    "    a_{20} & a_{21}& a_{22}\n",
    "    \\end{pmatrix}\n",
    "        \\begin{pmatrix}\n",
    "        x_0\\\\\n",
    "        x_1\\\\\n",
    "        x_2\\end{pmatrix} =\n",
    "        \\begin{pmatrix}\n",
    "         b_0\\\\\n",
    "         b_1\\\\\n",
    "         b_2\n",
    "         \\end{pmatrix} .$$\n",
    "\n",
    "    Now the vector $\\vec{b}$ on the\n",
    "    [*RHS*](http://www.science.oregonstate.edu/~rubin/Books/CPbook/eBook/GlossarySound/rhs.wav)\n",
    "    is assumed known, and the problem is to solve for the vector\n",
    "    $\\vec{x}$. Use an appropriate subroutine to solve these equations\n",
    "    for the three different $\\vec{x}$ vectors appropriate to these\n",
    "    three different $\\vec{b}$ values on the RHS:\n",
    "\n",
    "    $$\\begin{align}\n",
    "    \\vec{b}_{1} =\n",
    "    \\begin{pmatrix}\n",
    "    +12\\\\\n",
    "    -25\\\\\n",
    "     +32 \\end{pmatrix}\\!, \\quad  \\vec{b}_{2} =\\begin{pmatrix}\n",
    "     +4\\\\\n",
    "     -10\\\\\n",
    "     +22\\end{pmatrix}\\!, \\quad  \\vec{b}_{3} =\\begin{pmatrix}\n",
    "     +20\\\\\n",
    "     -30\\\\\n",
    "     +40\n",
    "     \\end{pmatrix}\\!.\\end{align}$$\n",
    "\n",
    "    The solutions should be\n",
    "\n",
    "    $$\\tag*{6.36}\n",
    "    \\vec{x}_{1} =\\begin{pmatrix}\n",
    "    +1\\\\\n",
    "    -2\\\\\n",
    "    +4\\end{pmatrix}\\!, \\quad \\vec{x}_{2} =\\begin{pmatrix}\n",
    "    +0.312\\\\\n",
    "    -0.038\\\\\n",
    "    +2.677\\end{pmatrix}\\!, \\quad \\vec{x}_{3} =\\begin{pmatrix}\n",
    "    +2.319 \\\\\n",
    "    -2.965 \\\\\n",
    "    +4.790\\end{pmatrix}\\!.$$\n",
    "\n",
    "3. Consider the matrix \n",
    "    \n",
    "      $A =\\begin{pmatrix}\n",
    "        \\alpha & \\beta\\\\\n",
    "       -\\beta & \\alpha\\\\\n",
    "        \\end{pmatrix}$\n",
    "    \n",
    "    where you are free to use\n",
    "    any values you want for *α* and *β*. Use a numerical eigenproblem\n",
    "    solver to show that the eigenvalues and eigenvectors are the complex\n",
    "    conjugates\n",
    "\n",
    "    $$\\tag*{6.37}\n",
    "    {\\vec{x}_{1,2}} =\\begin{pmatrix}\n",
    "    +1 \\\\\n",
    "    \\mp i\n",
    "    \\end{pmatrix} , \\quad \\lambda_{1,2}=\\alpha \\mp i\\beta .$$\n",
    "\n",
    "4.  Use your eigenproblem solver to find the eigenvalues of the matrix\n",
    "\n",
    "    $$\\tag*{6.38}\n",
    "     {A}=\\begin{pmatrix}\n",
    "        -2&+2&-3\\\\\n",
    "        +2&+1&-6\\\\\n",
    "        -1&-2&+0  \\end{pmatrix}.$$\n",
    "\n",
    "    -   Verify that you obtain the eigenvalues *λ*<sub>1</sub> = 5,\n",
    "        *λ*<sub>2</sub> = *λ*<sub>3</sub> = −3. Notice that double roots\n",
    "        can cause problems. In particular, there is a uniqueness issue\n",
    "        with their eigenvectors because any combination of these\n",
    "        eigenvectors is also an eigenvector.\n",
    "\n",
    "    -   Verify that the eigenvector for *λ*<sub>1</sub> = 5 is\n",
    "        proportional to\n",
    "\n",
    "        $$\\tag*{6.39}\n",
    "        \\vec{x}_{1} = \\frac{1}{\\sqrt{6}}\\begin{pmatrix}\n",
    "            -1\\\\\n",
    "            -2\\\\\n",
    "            +1\n",
    "            \\end{pmatrix}.$$\n",
    "\n",
    "    -   The eigenvalue −3 corresponds to a double root. This means that\n",
    "        the corresponding eigenvectors are degenerate, which in turn\n",
    "        means that they are not unique. Two linearly independent ones\n",
    "        are\n",
    "\n",
    "        $$\\tag*{6.40}\n",
    "        \\vec{x}_{2} =        \\frac{1} {\\sqrt{5}}\n",
    "        \\begin{pmatrix}\n",
    "            -2\\\\\n",
    "            +1\\\\\n",
    "            +0 \\end{pmatrix}, \\quad \\vec{x}_{3} = \\frac{1} {\\sqrt{10}}\\begin{pmatrix}\n",
    "            3\\\\\n",
    "            0\\\\\n",
    "         1 \\end{pmatrix}.$$\n",
    "\n",
    "        In this case it’s not clear what your eigenproblem solver will\n",
    "        give for the eigenvectors. Try to find a relationship between\n",
    "        your computed eigenvectors with the eigenvalue −3 and these two\n",
    "        linearly independent ones.\n",
    "\n",
    "5.  Imagine that your model of some physical system results in *N* = 100\n",
    "    coupled linear equations in *N* unknowns:\n",
    "\n",
    "    $$\\begin{align}\n",
    "    a_{00}y_0 + a_{01}y_1 + \\cdots + a_{0(N-1)}y_{N-1}  & =  b_0 , \\\\\n",
    "    a_{10}y_0 + a_{11}y_1 + \\cdots + a_{1(N-1)}y_{N-1} & =  b_1,\\\\\n",
    "    \\cdots   &  \\\\\n",
    "    a_{(N-1)0}y_0 + a_{(N-1)1}y_1 + \\cdots + a_{(N-1)(N-1)}y_{N-1} & =  b_{N-1}.\n",
    "         \\end{align}$$\n",
    "\n",
    "    In many cases, the *a* and *b* values are known, so your exercise is\n",
    "    to solve for all the *x* values, taking *a* as the *Hilbert* matrix\n",
    "    and $\\vec{b}$ as its first column:\n",
    "\n",
    "    $$\\begin{align}\n",
    "    &\\displaystyle [a_{ij}]        = {a}= \\left[ \\frac{1} {\n",
    "    i+j-1}\\right]\n",
    "     =\\begin{pmatrix}\n",
    "    1 &  \\frac{1}{2} & \\frac{1} { 3} &\\frac{1}{4} & \\cdots &\n",
    "    \\frac{1}{100} \\\\\n",
    "    \\frac{1}{2} & \\frac{1}{3}&  \\frac{1}{4} & \\frac{1}{5}& \\cdots &\n",
    "    \\frac{1}{101}\\\\\n",
    "    \\ddots\\\\\n",
    "    \\frac{1}{100} &\\frac{1}{101} &\\cdots & & \\cdots&\n",
    "    \\frac{1}{199}\\end{pmatrix},\\\\\n",
    "    &\\displaystyle[b_{i}] = \\vec{b} = \\left[\\frac{1}{i}\\right] =\n",
    "    \\begin{pmatrix}\n",
    "    \\\\\n",
    "    1 \\\\\n",
    "     \\frac{1}{2} \\\\\n",
    "     \\frac{1}{3}\\\\\n",
    "    \\ddots \\\\\n",
    "    \\frac{1}{100}\\\\\n",
    "    \\end{pmatrix}.\\end{align}$$\n",
    "\n",
    "    Compare to the analytic solution\n",
    "\n",
    "    $$\\tag*{6.41}\n",
    "    \\begin{pmatrix}\n",
    "    y_1\\\\\n",
    "    y_2 \\\\\n",
    "    \\ddots\\\\\n",
    "    y_N\\end{pmatrix} =\\begin{pmatrix}\n",
    "    1 \\\\\n",
    "    0 \\\\\n",
    "    \\ddots \\\\\n",
    "    0\\end{pmatrix}.$$\n",
    "\n",
    "6.  **Dirac Gamma Matrices** The Dirac equation extends quantum\n",
    "    mechanics to include relativity and spin 1/2. The extension of the\n",
    "    Hamiltonian operator for an electron requires it to contain\n",
    "    matrices, and those matrices are expressed in terms of 4 × 4 *γ*\n",
    "    matrices that can be represented in terms of the familiar 2 × 2\n",
    "    Pauli matrices *σ*<sub>i</sub>:\n",
    "    \n",
    "    $$\\begin{align}\n",
    "    \\tag*{6.42}\n",
    "    \\gamma_{i} & =\\begin{pmatrix}\n",
    "    0 &\\sigma_i\\\\\n",
    "    -\\sigma_i&nn0\\\\\n",
    "    \\end{pmatrix}\\!,\\ \\ \\ i=1,2,3, \\\\\n",
    "    \\sigma_{1} & =\\begin{pmatrix}\n",
    "    0 &1\\\\\n",
    "    1&0\\\\ \\end{pmatrix}\\!, \\quad  \\sigma_{2} =\\begin{pmatrix}\n",
    "     0&-i\\\\\n",
    "     i&0\\\\ \\end{pmatrix}\\!, \\quad  \\sigma_{3} =\\begin{pmatrix}\n",
    "     1&0\\\\\n",
    "     0&-1\\\\ \\end{pmatrix}\\!.\\tag*{6.43}\\end{align}$$\n",
    "     \n",
    "     Confirm the following properties of the *γ* matrices:\n",
    "\n",
    "$$\\begin{align}\n",
    "\\tag*{6.44}\n",
    "\\gamma_2^\\dagger=\\gamma_2^{-1} & = -\\gamma_2, \\\\\n",
    "\\gamma_1 \\gamma_2  & =-i\\begin{pmatrix}\n",
    "\\sigma_3&0\\\\\n",
    "0&\\sigma_3\\\\\n",
    "\\end{pmatrix}\\! .\\tag*{6.45}\\end{align}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### NewtonNDanimate.py, Notebook Version"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div id=\"scene0\"><div id=\"glowscript\" class=\"glowscript\"></div></div>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/javascript": [
       "require.undef(\"nbextensions/glow.1.0.min\");"
      ],
      "text/plain": [
       "<IPython.core.display.Javascript object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/javascript": [
       "require.undef(\"nbextensions/jquery-ui.custom.min\");"
      ],
      "text/plain": [
       "<IPython.core.display.Javascript object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/javascript": [
       "require.undef(\"nbextensions/glow.1.1.min\");"
      ],
      "text/plain": [
       "<IPython.core.display.Javascript object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/javascript": [
       "require.undef(\"nbextensions/glowcomm\");"
      ],
      "text/plain": [
       "<IPython.core.display.Javascript object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/javascript": [
       "require([\"nbextensions/glowcomm\"], function(){console.log(\"glowcomm loaded\");})"
      ],
      "text/plain": [
       "<IPython.core.display.Javascript object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# NewtonNDanimate.py, Notebook Version\n",
    "\n",
    "#from __future__ import division,print_function\n",
    "from ivisual import *\n",
    "from IPython.display import IFrame\n",
    "from numpy.linalg import solve\n",
    "from numpy import array, take\n",
    "\n",
    "scene=canvas(title=\"Strings and masses configuration\")\n",
    "scene.width=500\n",
    "scene.height=500\n",
    "tempe = curve(x=range(0,500),color=color.black)\n",
    "n = 9\n",
    "eps = 1e-6\n",
    "deriv = zeros( (n, n), float)\n",
    "f = zeros( (n), float)\n",
    "x = array([0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 1., 1., 1.])\n",
    "def plotconfig():\n",
    "    for obj in scene.objects:\n",
    "        obj.visible=0               #to erase the previous configuration\n",
    "    L1=3.0\n",
    "    L2=4.0\n",
    "    L3=4.0\n",
    "    xa=L1*x[3]                   #L1*cos(th1)\n",
    "    ya=L1*x[0]                  #L1 sin(th1)\n",
    "    xb=xa+L2*x[4]               #L1*cos(th1)+L2*cos(th2)\n",
    "    yb=ya+L2*x[1]               #L1*sin(th1)+L2*sen(th2)\n",
    "    xc=xb+L3*x[5]                #L1*cos(th1)+L2*cos(th2)+L3*cos(th3)\n",
    "    yc=yb-L3*x[2]          #L1*sin(th1)+L2*sen(th2)-L3*sin(th3)\n",
    "    mx=100.0                #for linear coordinate transformation\n",
    "    bx=-500.0               # from 0=< x =<10\n",
    "    my=-100.0               #to    -500 =<x_window=>500\n",
    "    by=400.0                #same transformation for y\n",
    "    xap=mx*xa+bx            #to keep aspect ratio\n",
    "    yap=my*ya+by\n",
    "    ball1=sphere(pos=(xap,yap,0), color=color.cyan,radius=15) \n",
    "    xbp=mx*xb+bx\n",
    "    ybp=my*yb+by\n",
    "    ball2=sphere(pos=(xbp,ybp,0), color=color.cyan,radius=25) \n",
    "    xcp=mx*xc+bx\n",
    "    ycp=my*yc+by\n",
    "    x0=mx*0+bx\n",
    "    y0=my*0+by\n",
    "    line1=curve(pos=[(x0,y0),(xap,yap)], color=color.yellow,radius=4)\n",
    "    line2=curve(pos=[(xap,yap),(xbp,ybp)], color=color.yellow,radius=4)\n",
    "    line3=curve(pos=[(xbp,ybp),(xcp,ycp)], color=color.yellow,radius=4)\n",
    "    topline=curve(pos=[(x0,y0),(xcp,ycp)], color=color.red,radius=4)\n",
    "    rate(1)\n",
    "\n",
    "def F(x, f):                 # Define F function\n",
    "    f[0] = 3*x[3]  +  4*x[4]  +  4*x[5]  -  8.0\n",
    "    f[1] = 3*x[0]  +  4*x[1]  -  4*x[2]\n",
    "    f[2] = x[6]*x[0]  -  x[7]*x[1]  -  10.0\n",
    "    f[3] = x[6]*x[3]  -  x[7]*x[4]\n",
    "    f[4] = x[7]*x[1]  +  x[8]*x[2]  -  20.0\n",
    "    f[5] = x[7]*x[4]  -  x[8]*x[5]\n",
    "    f[6] = pow(x[0], 2)  +  pow(x[3], 2)  -  1.0\n",
    "    f[7] = pow(x[1], 2)  +  pow(x[4], 2)  -  1.0\n",
    "    f[8] = pow(x[2], 2)  +  pow(x[5], 2)  -  1.0\n",
    "    \n",
    "def dFi_dXj(x, deriv, n): # Define derivative function\n",
    "    h = 1e-4                                             \n",
    "    for j in range(0, n):\n",
    "        temp = x[j]\n",
    "        x[j] = x[j] +  h/2.\n",
    "        F(x, f)                                                 \n",
    "        for i in range(0, n):  deriv[i, j] = f[i] \n",
    "        x[j] = temp\n",
    "         \n",
    "    for j in range(0, n):\n",
    "        temp = x[j]\n",
    "        x[j] = x[j] - h/2.\n",
    "        F(x, f)\n",
    "        for i in range(0, n): deriv[i, j] = (deriv[i, j] - f[i])/h\n",
    "        x[j] = temp\n",
    "\n",
    "for it in range(1, 100):\n",
    "    rate(1)         # 1 second between graphs\n",
    "    F(x, f)                              \n",
    "    dFi_dXj(x, deriv, n)   \n",
    "    B = array([[ - f[0]], [ - f[1]], [ - f[2]], [ - f[3]], [ - f[4]], [ - f[5]], [ - f[6]], [ - f[7]], [ - f[8]]])      \n",
    "    sol = solve(deriv, B)\n",
    "    dx = take(sol, (0, ), 1)             # take the first column of matrix sol\n",
    "    for i in range(0, n):\n",
    "        x[i]  = x[i]  +  dx[i]\n",
    "    plotconfig()\n",
    "    errX = errF = errXi = 0.0\n",
    "      \n",
    "    for i in range(0, n):\n",
    "        if ( x[i] !=  0.): errXi = abs(dx[i]/x[i])\n",
    "        else:  errXi = abs(dx[i])\n",
    "        if ( errXi > errX): errX = errXi                            \n",
    "        if ( abs(f[i]) > errF ):  errF = abs(f[i])        \n",
    "        if ( (errX <=  eps) and (errF <=  eps) ): break\n",
    "        \n",
    "print('Number of iterations = ', it)\n",
    "print('Solution:')\n",
    "for i in range(0, n):\n",
    "        print('x[', i, '] =  ', x[i])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 6.6.1  Matrix Solution of the String Problem<a id=\"6.6.1\"></a>\n",
    "\n",
    "In § 6.1 we set up the solution to our problem of two masses on a\n",
    "string. Now we have the matrix tools needed to solve it. Your\n",
    "**problem** is to check out the physical reasonableness of the solution\n",
    "for a variety of weights and lengths. You should check that the deduced\n",
    "tensions are positive and that the deduced angles correspond to a\n",
    "physical geometry (e.g., with a sketch). Because this is a physics-based\n",
    "problem, we know that the sine and cosine functions must be less than 1\n",
    "in magnitude and that the tensions should be similar in magnitude to the\n",
    "weights of the spheres. Our solution `NewtonNDanimate.py` (Listing 6.1)\n",
    "shows graphically the step-by-step search for the solution.\n",
    "\n",
    "[**Listing 6.1  NewtonNDanimate.py**](http://www.science.oregonstate.edu/~rubin/Books/CPbook/Codes/PythonCodes/NewtonNDanimate.py) shows the step-by-step search for\n",
    "the solution of the two-mass-on-a-string problem via Newton-Raphson\n",
    "search.\n",
    "\n",
    "### 6.6.2  Explorations<a id=\"6.6.2\"></a>\n",
    "\n",
    "1.  See at what point your initial guess for the angles of the strings\n",
    "    gets so bad that the computer is unable to find a physical solution.\n",
    "\n",
    "2.  A possible problem with the formalism we have just laid out is that\n",
    "    by incorporating the identity\n",
    "    sin<sup>2</sup>*θ*<sub>*i*</sub> + cos<sup>2</sup>*θ*<sub>*i*</sub> = 1\n",
    "    into the equations we may be discarding some information about the\n",
    "    sign of sin*θ* or cos*θ*. If you look at Figure 6.1, you can observe\n",
    "    that for some values of the weights and lengths, *θ*<sub>2</sub> may\n",
    "    turn out to be negative, yet cos*θ* should remain positive. We can\n",
    "    build this condition into our equations by replacing\n",
    "    *f*<sub>7</sub> − *f*<sub>9</sub> with *f*’s based on the form\n",
    "\n",
    "    $$f_7 =  x_4 -\\sqrt{1-x_1^2},\\quad f_8 = x_5 - \\sqrt{1-x_2^2},\\quad\n",
    "    f_9 =  x_6 - \\sqrt{1-x_3^2}.$$\n",
    "\n",
    "    See if this makes any difference in the solutions obtained.\n",
    "\n",
    "3.  ⊙ Solve the similar three-mass problem. The approach is the same,\n",
    "    but the number of equations is larger."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}