{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Chapter 7 - Direct solutions of linear equations" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Example No. 7_01 Pg No. 211" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "A =\n", "[[ 3 2 1 10]\n", " [ 0 2 2 8]\n", " [ 0 2 3 11]]\n", "After transformation\n", "A=\n", "[[ 3 2 1 10]\n", " [ 0 2 2 8]\n", " [ 0 0 1 3]]\n", "x = 1\n", "y = 1\n", "z = 3\n" ] } ], "source": [ "from numpy import mat\n", "#Elimination Process\n", "\n", "A = mat([[3, 2, 1, 10],[2, 3, 2, 14],[1, 2, 3, 14]])\n", "A[1,:] = A[1,:] - A[0,:]*A[1,0]/A[0,0]\n", "A[2,:] = A[2,:] - A[0,:]*A[2,0]/A[0,0]\n", "print 'A =\\n',A\n", "\n", "print 'After transformation'\n", "A[2,:] = A[2,:] - A[1,:]*A[2,1]/A[1,1]\n", "print 'A=\\n',A\n", "\n", "z = A[2,3]/A[2,2]\n", "y = (A[1,3] - A[1,2]*z)/A[1,1]\n", "x = (A[0,3] - A[0,1]*y - A[0,2]*z)/A[0,0]\n", "print 'x =',x \n", "print 'y =',y \n", "print 'z =',z " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Example No. 7_02 Pg No. 214" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[ 3 6 1 16]\n", " [ 2 4 3 13]\n", " [ 1 3 2 9]]\n", "since Aug(2,2) = 0 elimination is not possible,so reordering the matrix\n", "[[ 3 6 1 16]\n", " [ 2 4 3 13]\n", " [ 1 3 2 9]]\n", "Elimination is complete and by back substitution the solution is\n", "\n", "x3 = 1, x2 = 2 , x1 = 1 \n" ] } ], "source": [ "from numpy import mat,shape\n", "#Basic Gauss Elimination\n", "\n", "\n", "A = mat([[ 3, 6, 1],[ 2, 4 , 3],[ 1, 3, 2 ]])\n", "B = [16, 13, 9]\n", "S=shape(A)\n", "ar1,ac1 = S[0],S[1]\n", "Aug = mat([[3, 6 ,1 ,16],[2, 4, 3, 13],[1, 3, 2, 9]])\n", "for i in range(1,ar1):\n", " Aug[i,:] = Aug[i,:] - (Aug[i,0]/Aug[0,0])*Aug[0,:] \n", "\n", "print Aug\n", "print 'since Aug(2,2) = 0 elimination is not possible,so reordering the matrix'\n", "temp=A[2,:]\n", "A[2,:]=A[1,:]\n", "A[1,:]=temp\n", "print Aug\n", "print 'Elimination is complete and by back substitution the solution is\\n'\n", "print 'x3 = 1, x2 = 2 , x1 = 1 '" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Example No. 7_03 Pg No. 220" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Gauss Elimination using partial pivoting : ans is\n", "[ 9. -1. -10.]\n", "\n", " The solution can be obtained by back substitution \n", " x1 = 4 \n", " x2 = 0 \n", " x3 = -4 \n", "\n" ] } ], "source": [ "from numpy import array , zeros, dot, diag\n", "\n", "A = array([[ 2., 2., 1.],[4., 2., 3.],[ 1., -1., 1.]])\n", "B = array([[6.],[4.],[0.]])\n", "\n", "## Gauss Elimination using partial pivoting\n", "\n", "def GEPP(A, b):\n", " '''\n", " Gaussian elimination with partial pivoting.\n", " % input: A is an n x n nonsingular matrix\n", " % b is an n x 1 vector\n", " % output: x is the solution of Ax=b.\n", " % post-condition: A and b have been modified. \n", " '''\n", " n = len(A)\n", " if b.size != n:\n", " raise ValueError(\"Invalid argument: incompatible sizes between A & b.\", b.size, n)\n", " # k represents the current pivot row. Since GE traverses the matrix in the upper \n", " # right triangle, we also use k for indicating the k-th diagonal column index.\n", " for k in xrange(n-1):\n", " #Choose largest pivot element below (and including) k\n", " maxindex = abs(A[k:,k]).argmax() + k\n", " if A[maxindex, k] == 0:\n", " raise ValueError(\"Matrix is singular.\")\n", " #Swap rows\n", " if maxindex != k:\n", " A[[k,maxindex]] = A[[maxindex, k]]\n", " b[[k,maxindex]] = b[[maxindex, k]]\n", " for row in xrange(k+1, n):\n", " multiplier = A[row][k]/A[k][k]\n", " #the only one in this column since the rest are zero\n", " A[row][k] = multiplier\n", " for col in xrange(k + 1, n):\n", " A[row][col] = A[row][col] - multiplier*A[k][col]\n", " #Equation solution column\n", " b[row] = b[row] - multiplier*b[k]\n", " #print A ;print b\n", " x = zeros(n)\n", " k = n-1\n", " x[k] = b[k]/A[k,k]\n", " while k >= 0:\n", " x[k] = (b[k] - dot(A[k,k+1:],x[k+1:]))/A[k,k]\n", " k = k-1\n", " return x\n", "Aug=GEPP(A,B)\n", "print 'Gauss Elimination using partial pivoting : ans is\\n',Aug\n", "\n", "#Back Substitution\n", "def forward_elimination(A, b, n):\n", " \"\"\"\n", " Calculates the forward part of Gaussian elimination.\n", " \"\"\"\n", " for row in range(0, n-1):\n", " for i in range(row+1, n):\n", " factor = A[i,row] / A[row,row]\n", " for j in range(row, n):\n", " A[i,j] = A[i,j] - factor * A[row,j]\n", "\n", " b[i] = b[i] - factor * b[row]\n", "\n", " \n", " return A, b\n", "\n", "def back_substitution(a, b, n):\n", " \"\"\"\"\n", " Does back substitution, returns the Gauss result.\n", " \"\"\"\n", " x = zeros((n,1))\n", " x[n-1] = b[n-1] / a[n-1, n-1]\n", " for row in range(n-2, -1, -1):\n", " sums = b[row]\n", " for j in range(row+1, n):\n", " sums = sums - a[row,j] * x[j]\n", " x[row] = sums / a[row,row]\n", " return x\n", "\n", "def gauss(A, b):\n", " \"\"\"\n", " This function performs Gauss elimination without pivoting.\n", " \"\"\"\n", " n = A.shape[0]\n", "\n", " # Check for zero diagonal elements\n", " if any(diag(A)==0):\n", " raise ZeroDivisionError(('Division by zero will occur; '\n", " 'pivoting currently not supported'))\n", "\n", " A, b = forward_elimination(A, b, n)\n", " return back_substitution(A, b, n)\n", "\n", "x = gauss(A, B)\n", "print ('\\n The solution can be obtained by back substitution \\n x1 = %i \\n x2 = %i \\n x3 = %i \\n'%(x[0], x[1], x[2]))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Example No. 7_04 Pg No. 228" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Aug = \n", "[[ 2 4 -6 -8]\n", " [ 1 3 1 10]\n", " [ 2 -4 -2 -12]]\n", "Aug = \n", "[[ 1 2 -3 -4]\n", " [ 1 3 1 10]\n", " [ 2 -4 -2 -12]]\n", "Aug = \n", "[[ 1 2 -3 -4]\n", " [ 0 1 4 10]\n", " [ 0 -8 4 -12]]\n", "Aug = \n", "[[ 1 2 -3 -4]\n", " [ 0 1 4 10]\n", " [ 0 -8 4 -12]]\n", "Aug = \n", "[[ 1 0 -11 -4]\n", " [ 0 1 4 10]\n", " [ 0 0 36 -12]]\n", "Aug = \n", "[[ 1 0 -11 -4]\n", " [ 0 1 4 10]\n", " [ 0 0 1 -1]]\n", "Aug = \n", "[[ 1 0 0 -4]\n", " [ 0 1 0 10]\n", " [ 0 0 1 -1]]\n" ] } ], "source": [ "from numpy import mat, shape\n", "#Gauss Jordan Elimination\n", "\n", "A = mat([[ 2, 4, -6],[1, 3, 1],[2, -4, -2 ]])\n", "B = mat([[ -8],[ 10],[ -12 ]])\n", "S=shape(A)\n", "ar, ac =S[0],S[1]\n", "Aug = mat([[ 2, 4, -6, -8],[ 1, 3, 1, 10],[ 2, -4, -2, -12 ]])\n", "print 'Aug = \\n',Aug\n", "\n", "\n", "for i in range(0,ar):\n", " Aug[i,i:ar+1] = Aug[i,i:ar+1]/Aug[i,i]\n", " print 'Aug = \\n',Aug\n", " for k in range(0,ar):\n", " if k != i:\n", " Aug[k,i:ar] = Aug[k,i:ar] - Aug[k,i]*Aug[i,i:ar]\n", " \n", "\n", " print 'Aug = \\n',Aug" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Example No. 7_05 Pg No. 234" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "U = \n", "[[ 3. 2. 1. ]\n", " [ 0. 1.66666667 1.33333333]\n", " [ 0. 0. 1.6 ]]\n", "L = \n", "[[ 1. 0. 0. ]\n", " [ 0.66666667 1. 0. ]\n", " [ 0.33333333 0.8 1. ]]\n", "\n", " z(1) = 10.000000 \n", "\n", "\n", " z(2) = 14.000000 \n", "\n", "\n", " z(3) = 10.666667 \n", "\n", "\n", " x(3) = 6.666667 \n", "\n", "\n", " x(2) = 8.400000 \n", "\n", "\n", " x(1) = 1.111111 \n", "\n" ] } ], "source": [ "from numpy import array, arange, zeros\n", "from scipy.linalg import lu\n", "from __future__ import division\n", "\n", "#DoLittle LU Decomposition\n", "\n", "A = array([[ 3, 2, 1],[2 , 3 , 2],[1, 2, 3 ]])\n", "B = array([[ 10],[14],[14 ]])\n", "n,n = A.shape\n", "L=lu(A)[1]\n", "U=lu(A)[2]\n", "print 'U = \\n',U\n", "print 'L = \\n',L\n", "\n", "z=zeros([3,1])\n", "#Forward Substitution\n", "for i in range(0,n):\n", " z[i,0] = B[i,0]\n", " for j in range(0,i-1):\n", " z[i,0] = z[i,0] - L[i,0]*z[j,0]; \n", " \n", " print '\\n z(%i) = %f \\n'%(i+1,z[i,0])\n", "\n", "x=zeros([3,1]) \n", "#Back Substitution\n", "for i in arange(n-1,-1,-1):\n", " x[i,0] = z[i,0]\n", " for j in arange(n-1,i+1,-1):\n", " x[i,0] = x[i,0] - U[i,j]*x[j,0]\n", " \n", " x[i,0] = x[i,0]/U[i,i]\n", " print '\\n x(%i) = %f \\n'%(i+1,x[i,0])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Example No. 7_06 Pg No. 242" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "U =\n", "[[ 1. 2. 3. ]\n", " [ 0. 2.82842712 7.77817459]\n", " [ 0. 0. 8.54400375]]\n" ] } ], "source": [ "from numpy import sqrt, mat,shape, zeros\n", "#Cholesky's Factorisation\n", "\n", "A = mat([[ 1, 2, 3],[2, 8 ,22],[3, 22, 82 ]])\n", "n= shape(A)[0]\n", "U=zeros([n,n])\n", "\n", "for i in range(0,n):\n", " for j in range(0,n):\n", " if i == j:\n", " U[i,i] = A[i,i]\n", " for k in range(0,i-1):\n", " U[i,i] = U[i,i]-U[k,i]**2 \n", " \n", " U[i,i] = sqrt(U[i,i])\n", " elif i < j:\n", " U[i,j] = A[i,j]\n", " for k in range(0,i-1):\n", " U[i,j] = U[i,j] - U[k,i]*U[k,j]\n", " \n", " U[i,j] = U[i,j]/U[i,i]\n", " \n", "print 'U =\\n',U" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Example No. 7_07 Pg No. 245" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", " x(1) = 5.000000 \n", " x(2) = 5.000000 \n", "\n", "\n", " x(1) = -15.000000 \n", " x(2) = -15.000000 \n", "\n", "x=\n", "[[ 20.]\n", " [-15.]]\n", "r=\n", "[[ -2.66453526e-12]\n", " [ 1.00000000e-02]]\n" ] } ], "source": [ "from numpy import mat\n", "#Ill-Conditioned Systems\n", "\n", "A = mat([[ 2, 1],[2.001, 1]])\n", "B = mat([[ 25],[25.01 ]])\n", "x=[0,0]\n", "x[0] = (25 - 25.01)/(2 - 2.001)\n", "x[1] = (25.01*2 - 25*2.001)/(2*1 - 2.001*1)\n", "print '\\n x(1) = %f \\n x(2) = %f \\n'%(x[1],x[1])\n", "x[0] = (25 - 25.01)/(2 - 2.0005)#\n", "x[1] = (25.01*2 - 25*2.0005)/(2*1 - 2.0005*1)#\n", "print '\\n x(1) = %f \\n x(2) = %f \\n'%(x[1],x[1])\n", "x=mat([[x[0]],[x[1]]])\n", "r = A*(x)-B\n", "print 'x=\\n',x\n", "print 'r=\\n',r" ] } ], "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.9" } }, "nbformat": 4, "nbformat_minor": 0 }