summaryrefslogtreecommitdiff
path: root/Numerical_Methods_by_E._Balaguruswamy/chapter7.ipynb
diff options
context:
space:
mode:
Diffstat (limited to 'Numerical_Methods_by_E._Balaguruswamy/chapter7.ipynb')
-rw-r--r--Numerical_Methods_by_E._Balaguruswamy/chapter7.ipynb524
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
+}