diff options
Diffstat (limited to 'Numerical_Methods_by_E._Balaguruswamy/chapter7.ipynb')
-rw-r--r-- | Numerical_Methods_by_E._Balaguruswamy/chapter7.ipynb | 524 |
1 files changed, 524 insertions, 0 deletions
diff --git a/Numerical_Methods_by_E._Balaguruswamy/chapter7.ipynb b/Numerical_Methods_by_E._Balaguruswamy/chapter7.ipynb new file mode 100644 index 00000000..3dff4606 --- /dev/null +++ b/Numerical_Methods_by_E._Balaguruswamy/chapter7.ipynb @@ -0,0 +1,524 @@ +{ + "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 +} |