{
 "metadata": {
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np\n",
      "import numpy.linalg as la"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "n = 4 # crank up--8 is a good example\n",
      "\n",
      "nums = np.arange(n).astype(np.float64)\n",
      "nums"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 2,
       "text": [
        "array([ 0.,  1.,  2.,  3.])"
       ]
      }
     ],
     "prompt_number": 2
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "B = 1 + nums[:, np.newaxis] + nums\n",
      "B"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 3,
       "text": [
        "array([[ 1.,  2.,  3.,  4.],\n",
        "       [ 2.,  3.,  4.,  5.],\n",
        "       [ 3.,  4.,  5.,  6.],\n",
        "       [ 4.,  5.,  6.,  7.]])"
       ]
      }
     ],
     "prompt_number": 3
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "A = 1/B\n",
      "\n",
      "from scipy.linalg import hilbert\n",
      "A2 = hilbert(n)\n",
      "\n",
      "print \"difference: %g\" % la.norm(A-A2)\n",
      "print \"cond: %g\" % np.linalg.cond(A)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "difference: 0\n",
        "cond: 15513.7\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stderr",
       "text": [
        "/usr/lib/pymodules/python2.7/numpy/oldnumeric/__init__.py:11: ModuleDeprecationWarning: The oldnumeric module will be dropped in Numpy 1.9\n",
        "  warnings.warn(_msg, ModuleDeprecationWarning)\n"
       ]
      }
     ],
     "prompt_number": 4
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "b = np.ones((n,))\n",
      "\n",
      "b1 = b.copy()\n",
      "b1[-1] = 2.0/3.0\n",
      "x1 = la.solve(A, b1)\n",
      "print \"residual of x1: %g\"%(la.norm(b - np.dot(A, x1)))\n",
      "\n",
      "b2 = b.copy()\n",
      "b2[-1] = 0.667\n",
      "x2 = la.solve(A, b2)\n",
      "print \"residual of x2: %g\"%(la.norm(b - np.dot(A, x2)))\n",
      "\n",
      "print \"rel err in RHS: %g\"%(la.norm(b2-b1)/la.norm(b1))\n",
      "print \"rel err in solution: %g\"%(la.norm(x2-x1)/la.norm(x1))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "residual of x1: 0.333333\n",
        "residual of x2: 0.333\n",
        "rel err in RHS: 0.000179605\n",
        "rel err in solution: 0.0011524\n"
       ]
      }
     ],
     "prompt_number": 5
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x = np.ones((n,))\n",
      "b = np.dot(A, x)\n",
      "x3 = la.solve(A,b)\n",
      "print \"residual of x3: %g\"%(la.norm(b - np.dot(A, x3)))\n",
      "print \"rel err in solution: %g\"%(la.norm(x3-x)/la.norm(x))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "residual of x3: 1.11022e-16\n",
        "rel err in solution: 3.86179e-13\n"
       ]
      }
     ],
     "prompt_number": 6
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}