{
"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
}