{ "metadata": { "name": "" }, "nbformat": 3, "nbformat_minor": 0, "worksheets": [ { "cells": [ { "cell_type": "heading", "level": 1, "metadata": {}, "source": [ "Gaussian elimination II" ] }, { "cell_type": "heading", "level": 2, "metadata": {}, "source": [ "Importing functions we wrote in the last class : packages." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In python, one can just store a bunch of function definitions in a file. Then the file can be imported like any other package. We have stored our functions in `gauss_elim_old_fn.py`. We import it as" ] }, { "cell_type": "code", "collapsed": false, "input": [ "from gauss_elim_old_fn import *" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "[[1.0, 2.0, 3.0], [2.0, 1.0, 2.0], [3.0, 2.0, 1.0]]\n", " 1 2 3\n", " 2 1 2\n", " 3 2 1\n", " 1 2 3\n", " 2 1 2\n", " 3 2 1\n", "--------------------------\n", " 4 8 12\n", " 2 1 2\n", " 3 2 1\n", "-------------------------\n", " 4 8 12\n", " 2 1 2\n", " 3 2 1\n", " 4 8 12\n", " 2 1 2\n", " 3 2 1\n", "--------------------\n", " 0 6 8\n", " 2 1 2\n", " 3 2 1\n", "---------------------\n", " 0 6 8\n", " 2 1 2\n", " 3 2 1\n", "(1, 1)\n", "The matrix is the zero matrix.\n", "(-1, -1)\n", "The matrix is the zero matrix.\n", "(-1, -1)\n", "Please input a natural number.\n", "None\n", "1\n", "1\n", "2\n", "3\n", "5\n", "8\n", "13\n", "21\n", "34\n", "55\n", "Please input an positive integer.\n", "None\n", "Please input a natural number.\n", "None\n", "1\n", "1\n", "2\n", "3\n", "5\n", "8\n", "13\n", "21\n", "34\n", "55\n", "Please input a natural number.\n", "None\n", "9.33262154439e+157 Expected a non-negative integer.\n", "None Expected a non-negative integer.\n", "None 93326215443944152681699238856266700490715968264381621468592963895217599993229915608941463976156518286253697920827223758251185210916864000000000000000000000000\n", "(2, 1)\n", " 1 2 3 4\n", " 5 6 7 8\n", " 9 10 11 12\n", "\n", " 1 2 3 4\n", " 5 6 7 8\n", " 9 10 11 12\n", " 1 2 3 4\n", " 5 6 7 8\n", " 9 10 11 12\n", "---------------\n", " 1 2 3 4\n", " 5 6 7 -2\n", " 9 10 11 -3\n", " 0.00000 0.50000 1.00000 0.00000\n", " 2.00000 11.00000 0.00000 0.00000\n", " 1.00000 0.00000 0.00000 9.00000\n", "--------------------------------------------------\n", "The matrix is empty.\n", "Nothing to copy: S is empty\n", " 1.00000 5.50000 0.00000 0.00000\n", " 0.00000 1.00000 2.00000 0.00000\n", " 0.00000 0.00000 1.00000 0.81818\n" ] } ], "prompt_number": 1 }, { "cell_type": "code", "collapsed": false, "input": [ "def gauss_elim_prelim(matrix) :\n", " row = 0\n", " row_offset = 0\n", " col_offset = 0\n", " \n", " (r, c) = first_non_zero(matrix)\n", " if r == -1 :\n", " print \"Nothing to do: The matrix is zero.\"\n", " else :\n", " matrix = elem1(matrix, 0, r)\n", " matrix = sweep(matrix, (0, c))\n", " print_matrix(matrix, \"%6.2f\")" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 2 }, { "cell_type": "code", "collapsed": false, "input": [ "gauss_elim_prelim([[0,0,2,2],\\\n", " [1,4,4,4],\\\n", " [0,1,3,3]])" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ " 1.00 4.00 4.00 4.00\n", " 0.00 0.00 2.00 2.00\n", " 0.00 1.00 3.00 3.00\n" ] } ], "prompt_number": 3 }, { "cell_type": "code", "collapsed": false, "input": [ "def gauss_elim(matrix) :\n", " row = 0\n", " row_offset = 0\n", " col_offset = 0\n", " no_of_rows = len(matrix)\n", " if no_of_rows == 0 :\n", " print \"The matrix is empty.\"\n", " else :\n", " no_of_cols = len(matrix[0])\n", " if no_of_cols == 0 :\n", " print \"The matrix is empty.\"\n", " else :\n", " while row < no_of_rows :\n", " submatrix = []\n", " # Copy the submatrix to the right, bottom of (row_offset, col_offset)\n", " # to submatrix\n", " for i in range(no_of_rows - row_offset) :\n", " row_submatrix = []\n", " for j in range(no_of_cols - col_offset) :\n", " row_submatrix.append(matrix[i + row_offset][j + col_offset])\n", " submatrix.append(row_submatrix)\n", " #print_matrix(submatrix, \"%6.2f\")\n", " # print row, row_offset, col_offset, submatrix, matrix\n", " # print \"-\"*20\n", " \n", " (r, c) = first_non_zero(submatrix)\n", " if r != -1 :\n", " # print \"Nothing to do: The (sub)matrix is zero.\"\n", " # else :\n", " r += row_offset\n", " c += col_offset\n", " matrix = elem1(matrix, row_offset, r)\n", " matrix = sweep(matrix, (row_offset, c))\n", " row_offset += 1\n", " col_offset = c+1\n", " row += 1\n", " return matrix" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 4 }, { "cell_type": "code", "collapsed": false, "input": [ "B = read_matrix(\"files/B.txt\")\n", "print_matrix(B, \"%6.2f\")\n", "print \"-\"*40\n", "print_matrix(gauss_elim(B), \"%6.2f\")\n", "print \"=\"*40\n", "C = [[0,0,0,0], [0,1,1,0], [0,0,1,0]]\n", "print_matrix(C, \"%3d\")\n", "print \"-\"*40\n", "print_matrix(gauss_elim(C), \"%3d\")\n" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ " 0.00 2.00 3.00 2.00\n", " 1.00 4.00 5.00 5.00\n", " 3.00 0.00 2.00 6.00\n", "----------------------------------------\n", " 1.00 0.00 0.00 1.60\n", " 0.00 1.00 0.00 0.10\n", " 0.00 0.00 1.00 0.60\n", "========================================\n", " 0 0 0 0\n", " 0 1 1 0\n", " 0 0 1 0\n", "----------------------------------------\n", "The matrix is the zero matrix.\n", " 0 1 0 0\n", " 0 0 1 0\n", " 0 0 0 0\n" ] } ], "prompt_number": 5 }, { "cell_type": "heading", "level": 1, "metadata": {}, "source": [ "Solving a system of linear equations" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "First suppose we are given an equation of the type $Ax = b$ where $A$ is an $m \\times n$ matrix and $x$ is an $n$-vector and $b$ is an $m$-vector. We first form the matrix $A:b$ which is obtained by sticking $b$ after $A$." ] }, { "cell_type": "code", "collapsed": false, "input": [ "def join(A, B) :\n", " \"\"\"Returns A:B\"\"\"\n", " no_rows_A = len(A)\n", " no_rows_B = len(B)\n", " \n", " if no_rows_A != no_rows_B :\n", " print \"Matrices should have the same number of rows to be joined.\"\n", " retM = None\n", " else :\n", " retM = []\n", " if no_rows_A == 0 :\n", " retM = A\n", " else :\n", " no_cols_A = len(A[0])\n", " no_cols_B = len(B[0])\n", " for i in range(no_rows_A) :\n", " Mrow = []\n", " for j in range(no_cols_A) :\n", " Mrow.append(A[i][j])\n", " for k in range(no_cols_B) :\n", " Mrow.append(B[i][k])\n", " retM.append(Mrow)\n", " return retM" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 6 }, { "cell_type": "code", "collapsed": false, "input": [ "join([[4],[5]], [[4,5],[9,8]])" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 7, "text": [ "[[4, 4, 5], [5, 9, 8]]" ] } ], "prompt_number": 7 }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### To solve $Ax = b$\n", "Gaussian elimination of $A:b$, if $A$ is *invertible* is just finding a matrix $P$ such that $P(A:b)$ is of the form $I_3:c$. Then $P$ will be the inverse of $A$ and $c = Pb$ will be the solution." ] }, { "cell_type": "code", "collapsed": false, "input": [ "# Let A = \n", "A = [[0, 2, 3], [1, 0, 1], [0, 0, 1]]\n", "b = [[5], [2], [1]]\n", "# A:b = \n", "Ab = join(A, b)\n", "print_matrix(gauss_elim(Ab), \"%5.2g\")\n" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ " 1 0 0 1\n", " 0 1 0 1\n", " 0 0 1 1\n" ] } ], "prompt_number": 8 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Thus $x = \\begin{pmatrix}1\\\\1\\\\1\\end{pmatrix}$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To invert a matrix, one just appends identity on the right and does Gaussian elimination." ] }, { "cell_type": "code", "collapsed": false, "input": [ "def idmat(n) :\n", " \"\"\"returns the identity matrix of size n\"\"\"\n", " rmat = []\n", " for i in range(n) :\n", " row = []\n", " for j in range(n) :\n", " if i == j :\n", " row.append(1)\n", " else :\n", " row.append(0)\n", " rmat.append(row)\n", " return rmat" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 9 }, { "cell_type": "code", "collapsed": false, "input": [ "print_matrix(idmat(4), \"%2d\")" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ " 1 0 0 0\n", " 0 1 0 0\n", " 0 0 1 0\n", " 0 0 0 1\n" ] } ], "prompt_number": 10 }, { "cell_type": "code", "collapsed": false, "input": [ "def invert(A) :\n", " size = len(A)\n", " bigmat = join(A,idmat(size))\n", " biginv = gauss_elim(bigmat)\n", " inv = []\n", " for i in range(size) :\n", " row = []\n", " for j in range(size) :\n", " row.append(biginv[i][j + size])\n", " inv.append(row)\n", " return inv" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 11 }, { "cell_type": "code", "collapsed": false, "input": [ "UT1 = [[1,1,1],[0,1,1],[0,0,1]]\n", "print \"Inverse of\"\n", "print_matrix(UT1, \"%3d\")\n", "print \"is\"\n", "print_matrix(invert(UT1), \"%3d\")" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "Inverse of\n", " 1 1 1\n", " 0 1 1\n", " 0 0 1\n", "is\n", " 1 -1 0\n", " 0 1 -1\n", " 0 0 1\n" ] } ], "prompt_number": 12 } ], "metadata": {} } ] }