LinearSystems.ipynb

{
 "metadata": {
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "ESE - Numerical Methods I: Systems of Linear Equations"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%pylab inline\n",
      "import numpy as np\n",
      "\n",
      "# scipy linear algebra module\n",
      "import scipy.linalg as scal"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Populating the interactive namespace from numpy and matplotlib\n"
       ]
      }
     ],
     "prompt_number": 1
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Problems"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "1.)\n",
      "\n",
      "By evaluating the determinant, classify the following matrices as singular, ill conditioned,\n",
      "or well conditioned (using <a href=\"http://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.det.html#scipy.linalg.det\">`scal.det`</a>)\n",
      "\n",
      "\\begin{equation}\n",
      "A_{1} =\n",
      " \\begin{pmatrix}\n",
      "  1 & 2 & 3 \\\\\n",
      "  2 & 3 & 4 \\\\\n",
      "  3 & 4 & 5\n",
      " \\end{pmatrix}\n",
      "\\end{equation}\n",
      "\n",
      "\\begin{equation}\n",
      "A_{2} =\n",
      " \\begin{pmatrix}\n",
      "  2.11 & -0.80 & 1.77 \\\\\n",
      "  -1.84 & 3.03 & 1.21 \\\\\n",
      "  -1.57 & 5.25 & 4.30\n",
      " \\end{pmatrix}\n",
      "\\end{equation}\n",
      "\n",
      "\\begin{equation}\n",
      "A_{3} =\n",
      " \\begin{pmatrix}\n",
      "  2 & -1 & 0 \\\\\n",
      "  -1 & 2 & -1 \\\\\n",
      "  0 & -1 & 2\n",
      " \\end{pmatrix}\n",
      "\\end{equation}\n",
      "\n",
      "\\begin{equation}\n",
      "A_{4} =\n",
      " \\begin{pmatrix}\n",
      "  4 & 3 & -1 \\\\\n",
      "  7 & -2 & 3 \\\\\n",
      "  5 & -18 & 13\n",
      " \\end{pmatrix}\n",
      "\\end{equation}\n",
      "\n",
      "NB: You can use <a href=\"http://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.norm.html\">`scal.norm`</a> and the ratio of determinant to norm for this."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "A1 =  np.array([[1,2,3],[2,3,4],[3,4,5]])\n",
      "A2 =  np.array([[2.11,-0.80,1.77],[-1.84,3.03,1.21],[-1.57,5.25,4.3]])\n",
      "A3 =  np.array([[2,-1,0],[-1,2,-1,],[0,-1,2]])\n",
      "A4 =  np.array([[4,3,-1],[7,-2,3],[5,-18,13]])\n",
      "\n",
      "def conditioning(A):\n",
      "    d = scal.det(A)  # determinant\n",
      "    n = scal.norm(A) # frobenius matrix norm\n",
      "    cn = d/n\n",
      "    if cn == 0.:\n",
      "        cs = 'singular'\n",
      "    elif cn < 0.1:\n",
      "        cs = 'ill conditioned'\n",
      "    else:\n",
      "        cs = 'well conditioned'\n",
      "    print 'determinant is {0:.2e}, matrix norm is {1:.2e} and matrix is'.format(d,n), cs\n",
      "\n",
      "conditioning(A1)\n",
      "conditioning(A2)\n",
      "conditioning(A3)\n",
      "conditioning(A4)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "determinant is 2.22e-16, matrix norm is 9.64e+00 and matrix is ill conditioned\n",
        "determinant is 5.99e-01, matrix norm is 8.41e+00 and matrix is ill conditioned\n",
        "determinant is 4.00e+00, matrix norm is 4.00e+00 and matrix is well conditioned\n",
        "determinant is 0.00e+00, matrix norm is 2.46e+01 and matrix is singular\n"
       ]
      }
     ],
     "prompt_number": 2
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "2.)\n",
      "\n",
      "Solve the linear system $A x = b$ for x (using <a href=\"http://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.solve.html#scipy.linalg.solve\">`scal.solve`</a>). \n",
      "\n",
      "\\begin{equation}\n",
      "A =\n",
      " \\begin{pmatrix}\n",
      "  1 & 2 & 3 \\\\\n",
      "  4 & 5 & 6 \\\\\n",
      "  7 & 8 & 9\n",
      " \\end{pmatrix}\n",
      "\\end{equation}\n",
      "\n",
      "\\begin{equation}\n",
      "b =\n",
      " \\begin{pmatrix}\n",
      "  1  \\\\\n",
      "  2 \\\\\n",
      "  3\n",
      " \\end{pmatrix}\n",
      "\\end{equation}\n",
      "\n",
      "check the solution by assuring that `np.dot(A, x)-b` is numerically equal to `0.`, (using <a href=\"http://docs.scipy.org/doc/numpy/reference/generated/numpy.allclose.html\">`np.allclose`</a>)"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "A = np.array([[1,2,3], [4,5,6], [7,8,9]])\n",
      "b = np.array([1,2,3])\n",
      "\n",
      "x = scal.solve(A, b)\n",
      "print x\n",
      "np.allclose(np.dot(A, x)-b, 0.)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[-0.33333333  0.66666667  0.        ]\n"
       ]
      },
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 3,
       "text": [
        "True"
       ]
      }
     ],
     "prompt_number": 3
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "3.)\n",
      "\n",
      "Find $L$ and $U$ for\n",
      "\n",
      "\\begin{equation}\n",
      "A =\n",
      " \\begin{pmatrix}\n",
      "  4 & -1 & 0 \\\\\n",
      "  -1 & 4 & -1 \\\\\n",
      "  0 & -1 & 4\n",
      " \\end{pmatrix}\n",
      "\\end{equation}\n",
      "\n",
      "(using <a href=\"http://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.lu.html#scipy.linalg.lu\">`lu`</a>)"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "A = np.array([[1,-1,0], [-1,4,-1], [0,-1,4]])\n",
      "\n",
      "P,L,U = scal.lu(A)\n",
      "np.dot(L,U) # matrix multiplication"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 4,
       "text": [
        "array([[ 1., -1.,  0.],\n",
        "       [-1.,  4., -1.],\n",
        "       [ 0., -1.,  4.]])"
       ]
      }
     ],
     "prompt_number": 4
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "4.)\n",
      "\n",
      "Solve the linear system $A x = b$ for x (using the LU factorization <a href=\"http://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.lu_factor.html#scipy.linalg.lu_factor\">`scal.lu_factor`</a> and LU solver <a href=\"http://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.lu_solve.html#scipy.linalg.lu_solve\">`scal.lu_solve`</a>). \n",
      "\n",
      "\\begin{equation}\n",
      "A =\n",
      " \\begin{pmatrix}\n",
      "  1 & 2 & 3 \\\\\n",
      "  4 & 5 & 6 \\\\\n",
      "  7 & 8 & 9\n",
      " \\end{pmatrix}\n",
      "\\end{equation}\n",
      "\n",
      "\\begin{equation}\n",
      "b =\n",
      " \\begin{pmatrix}\n",
      "  1  \\\\\n",
      "  2 \\\\\n",
      "  3\n",
      " \\end{pmatrix}\n",
      "\\end{equation}\n",
      "\n",
      "check the solution by assuring that `np.dot(A, x)-b` is numerically equal to `0.`, (using <a href=\"http://docs.scipy.org/doc/numpy/reference/generated/numpy.allclose.html\">`np.allclose`</a>)"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "A = np.array([[1,2,3], [4,5,6], [7,8,9]])\n",
      "b = np.array([1,2,3])\n",
      "\n",
      "LU, P = scal.lu_factor(A)\n",
      "print LU\n",
      "x = scal.lu_solve((LU,P), b)\n",
      "print x\n",
      "np.allclose(np.dot(A, x)-b, 0.)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[[  7.00000000e+00   8.00000000e+00   9.00000000e+00]\n",
        " [  1.42857143e-01   8.57142857e-01   1.71428571e+00]\n",
        " [  5.71428571e-01   5.00000000e-01   1.11022302e-16]]\n",
        "[-0.33333333  0.66666667  0.        ]\n"
       ]
      },
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 5,
       "text": [
        "True"
       ]
      }
     ],
     "prompt_number": 5
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "5.) We now know the solution to \n",
      "\n",
      "\\begin{equation}\n",
      "A =\n",
      " \\begin{pmatrix}\n",
      "  2 & 3 & -1 \\\\\n",
      "  3 & 2 & -5 \\\\\n",
      "  1 & 4 & -1\n",
      " \\end{pmatrix}\n",
      "\\end{equation}\n",
      "\n",
      "\\begin{equation}\n",
      "b =\n",
      " \\begin{pmatrix}\n",
      "  1  \\\\\n",
      "  -1 \\\\\n",
      "  2\n",
      " \\end{pmatrix}\n",
      "\\end{equation}\n",
      "\n",
      "is \n",
      "\n",
      "\\begin{equation}\n",
      "x =\n",
      " \\begin{pmatrix}\n",
      "  1.84375  \\\\\n",
      "  0.40652   \\\\\n",
      "  1.46875\n",
      " \\end{pmatrix}\n",
      "\\end{equation}\n",
      "\n",
      "Without using functions from the `linalg` module (but e.g. `np.dot`), complete the code below so that the result is `True`."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "A = np.array([[2,-3,-1], [3,2,-5], [1,4,-1]],dtype=np.float) \n",
      "b = np.array([1,-1,2],dtype=np.float)\n",
      "\n",
      "def lu_dec(A):\n",
      "    \"\"\"\n",
      "    Returns the matrix A where [A] = [L\\U], i.e., [U] in the upper triangle and \n",
      "    the nondiagonal terms of [L] in the lower triangle.\n",
      "    \"\"\"\n",
      "    n = len(A)\n",
      "    for k in range(0,n-1):\n",
      "        for i in range(k+1,n):\n",
      "            if A[i,k] != 0.0:\n",
      "                lam = A[i,k]/A[k,k]\n",
      "                A[i,k+1:n] = A[i,k+1:n] - lam*A[k,k+1:n]\n",
      "                A[i,k] = lam\n",
      "    return A\n",
      "\n",
      "def lu_sol(LU,b):\n",
      "    \"\"\"\n",
      "    Returns the solution vector x as from [L][U]{x} = b, where [a] = [L\\U] is the matrix\n",
      "    returned from lu_dec.\n",
      "    \"\"\"\n",
      "    n = len(LU)\n",
      "    for k in range(1,n):\n",
      "        b[k] = b[k] - dot(LU[k,0:k],b[0:k])\n",
      "    for k in range(n-1,-1,-1):\n",
      "        b[k] = (b[k] - dot(LU[k,k+1:n],b[k+1:n]))/LU[k,k]\n",
      "    return b\n",
      "\n",
      "LU = lu_dec(A)\n",
      "x = lu_sol(LU, b)\n",
      "np.allclose(x, [1.84375,0.40625,1.46875])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 6,
       "text": [
        "True"
       ]
      }
     ],
     "prompt_number": 6
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 6
    }
   ],
   "metadata": {}
  }
 ]
}