diff options
Diffstat (limited to 'Numerical_Methods_by_E._Balaguruswamy')
15 files changed, 7656 insertions, 0 deletions
diff --git a/Numerical_Methods_by_E._Balaguruswamy/chapter10.ipynb b/Numerical_Methods_by_E._Balaguruswamy/chapter10.ipynb new file mode 100644 index 00000000..7de04b99 --- /dev/null +++ b/Numerical_Methods_by_E._Balaguruswamy/chapter10.ipynb @@ -0,0 +1,282 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Chapter 10 - Curve fitting regression" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 10_01 Pg No. 326" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "b = [-1 -1 0 0 1]\n", + "a = [ 3 3 1 1 -2]\n", + "y= [-x + 3 -x + 3 1 1 x - 2]\n" + ] + } + ], + "source": [ + "from numpy import array\n", + "from sympy.abc import x\n", + "#Fitting a Straight Line\n", + "\n", + "X = range(1,6)\n", + "Y = array([[ 3, 4, 5 ,6 ,8 ]])\n", + "n = len(X)#\n", + "X=array(X)\n", + "b = ( n*sum(X*Y) - sum(X)*sum(Y) )/( n*sum(X*X) - (sum(X))**2 )\n", + "a = sum(Y)/n - b*sum(X)/n\n", + "print 'b = ',b\n", + "print 'a = ',a\n", + "y = a + b*x\n", + "print 'y=',y" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 10_02 Pg No. 331" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "b = 2.0\n", + "lna = -0.69314718056\n", + "a = 0.5\n", + "\n", + " The power function equation obtained is \n", + " y = 0.5x**2\n" + ] + } + ], + "source": [ + "from numpy import array,log,exp\n", + "from sympy.abc import x\n", + "\n", + "#Fitting a Power-Function model to given data\n", + "\n", + "X = array(range(1,6))\n", + "Y = [ 0.5, 2 ,4.5 ,8 ,12.5 ]\n", + "Xnew = log(X)\n", + "Ynew = log(Y)\n", + "n = len(Xnew)\n", + "b = ( n*sum(Xnew*Ynew) - sum(Xnew)*sum(Ynew) )/( n*sum(Xnew*Xnew) - ( sum(Xnew) )**2 )\n", + "lna = sum(Ynew)/n - b*sum(Xnew)/n\n", + "a = exp(lna)\n", + "print 'b = ',b\n", + "print 'lna = ',lna\n", + "print 'a = ',a\n", + "print '\\n The power function equation obtained is \\n y = %.4Gx**%.4G'%(a,b)#" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 10_03 Pg No. 332" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "b = 37.6294062985\n", + "a = 20.9234245534\n", + "The relationship between T and t is \n", + "T = 37.63*e**(t/4) + 20.92 \n", + "\n", + "The temperature at t = 6 is 189.566723485\n" + ] + } + ], + "source": [ + "from numpy import array,log,exp\n", + "time = array(range(1,5))\n", + "T = [ 70 ,83 ,100 ,124 ]\n", + "t = 6\n", + "Fx = exp(time/4.0)\n", + "n = len(Fx)\n", + "Y = T #\n", + "b = ( n*sum(Fx*Y) - sum(Fx)*sum(Y) )/( n*sum(Fx*Fx) - (sum(Fx))**2 )\n", + "a = sum(Y)/n - b*sum(Fx)/n\n", + "print 'b = ',b\n", + "print 'a = ',a\n", + "print 'The relationship between T and t is \\nT = %.4G*e**(t/4) + %.4G \\n'%(b,a)\n", + "#deff('T = T(t)'%('T = b*exp(t/4) + a '\n", + "def T(t):\n", + " tt=b*exp(t/4.0)+a\n", + " return tt\n", + " \n", + "T_6 = T(6)\n", + "print 'The temperature at t = 6 is',T_6" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 10_04 Pg NO. 335" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Using CA = B form , we get\n", + "B [ 6. 62. 190.]\n", + "C [[ 1. 1. 4.]\n", + " [ 1. 4. 10.]\n", + " [ 4. 10. 30.]]\n", + "A = [[ 20. 103.33333333 -190. ]\n", + " [ 10. 144.66666667 -190. ]\n", + " [ -6. -62. 95. ]]\n", + "Therefore the least sqaures polynomial is\n", + " y = 1J + 1J*x + 1J*x**2 \n", + "[ 20. 103.33333333 -190. ]\n", + "[ 10. 144.66666667 -190. ]\n", + "[ -6. -62. 95.]\n" + ] + } + ], + "source": [ + "from numpy import array,ones,identity\n", + "from numpy.linalg import inv\n", + "\n", + "#Curve Fitting\n", + "\n", + "x = array(range(1,5))\n", + "y = [6, 11, 18, 27 ]\n", + "n = len(x) #Number of data points\n", + "m = 2+1 #Number of unknowns\n", + "print 'Using CA = B form , we get'\n", + "C=identity(m)\n", + "B=ones(m)\n", + "for j in range(0,m):\n", + " for k in range(0,m):\n", + " C[j,k] = sum(x**(j+k-2))\n", + " \n", + " B[j] = sum( y*( x**( j-1 ) ) )\n", + "\n", + "print 'B',B\n", + "print 'C',C\n", + "A = inv(C)*B\n", + "print 'A = ',A\n", + "print 'Therefore the least sqaures polynomial is\\n y = 1J + 1J*x + 1J*x**2 \\n',(A[0])\n", + "print A[1]\n", + "print A[2]" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 10_05 Pg No. 342" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "C=\n", + "[[ 4. 10. 6.]\n", + " [ 10. 30. 20.]\n", + " [ 6. 20. 14.]]\n", + "B=\n", + "[ 84. 240. 156.]\n", + "\n", + " The regression plane is \n", + " y = 5 + 6*x + 0*z \n", + "\n" + ] + } + ], + "source": [ + "from numpy import array,ones,identity,arange,vstack,transpose\n", + "from scipy.sparse.linalg import lsqr\n", + "#Plane Fitting\n", + "\n", + "x = range(1,5)\n", + "z = range(0,4)\n", + "y = arange(12,31,6)\n", + "n = len(x) #Number of data points\n", + "m = 3 #Number of unknowns\n", + "G = vstack([ones(n),x,z])\n", + "H = transpose(G)\n", + "C = G.dot(H)\n", + "B = y.dot(H)\n", + "D = lsqr(C,B)\n", + "print 'C=\\n',C\n", + "print 'B=\\n',B\n", + "print '\\n The regression plane is \\n y = %d + %.f*x + %d*z \\n'%(D[0][0],D[0][1],D[0][2])\n" + ] + } + ], + "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 +} diff --git a/Numerical_Methods_by_E._Balaguruswamy/chapter11.ipynb b/Numerical_Methods_by_E._Balaguruswamy/chapter11.ipynb new file mode 100644 index 00000000..c717bdbd --- /dev/null +++ b/Numerical_Methods_by_E._Balaguruswamy/chapter11.ipynb @@ -0,0 +1,428 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Chapter 11 - Numerical differentiation" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 11_01 Pg No. 348" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "h\ty\terr\n", + "0.2 \t2.2 \t0.2\n", + "0.1 \t2.1 \t0.1\n", + "0.05 \t2.05 \t0.05\n", + "0.01 \t2.01 \t0.01\n" + ] + } + ], + "source": [ + "from __future__ import division\n", + "from scipy.misc import derivative\n", + "\n", + "#First order forward difference\n", + "\n", + "def f(x):\n", + " F=x**2\n", + " return F\n", + " \n", + "\n", + "def df(x,h):\n", + " DF=(f(x+h)-f(x))/h\n", + " return DF\n", + "\n", + "dfactual = derivative(f,1)\n", + "h = [0.2, 0.1, 0.05, 0.01 ]\n", + "print 'h\\ty\\terr'\n", + "for hh in h:\n", + " y = df(1,hh)\n", + " err = y - dfactual\n", + " print hh,'\\t',y,'\\t',err" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 11_02 Pg No. 350" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "h\ty\terr\n", + "0.2 \t2.0 \t-4.4408920985e-16\n", + "0.1 \t2.0 \t4.4408920985e-16\n", + "0.05 \t2.0 \t4.4408920985e-16\n", + "0.01 \t2.0 \t1.7763568394e-15\n" + ] + } + ], + "source": [ + "from __future__ import division\n", + "from scipy.misc import derivative\n", + "\n", + "#Three-Point Formula\n", + "\n", + "def f(x):\n", + " F = x**2\n", + " return F\n", + "def df(x,h):\n", + " DF = (f(x+h)-f(x-h))/(2*h)\n", + " return DF\n", + "dfactual = derivative(f,1)\n", + "h = [0.2, 0.1, 0.05, 0.01 ]\n", + "print 'h\\ty\\terr'\n", + "for hh in h:\n", + " y = df(1,hh)\n", + " err = y - dfactual\n", + " print hh,'\\t',y,'\\t',err\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 11_03 Pg No. 353" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "h\ty\t\terr\n", + "0.01 \t0.898257285429 \t-0.00218981692372\n", + "0.015 \t0.897151155627 \t-0.00329594672573\n", + "0.02 \t0.896037563392 \t-0.00440953896077\n", + "0.025 \t0.894916522709 \t-0.00553057964367\n", + "0.03 \t0.893788047675 \t-0.00665905467759\n", + "0.035 \t0.892652152498 \t-0.00779494985432\n", + "0.04 \t0.891508851498 \t-0.00893825085448\n", + "\n", + "M2 = 2.061e-08 & hopt = 4.706e-01\n" + ] + } + ], + "source": [ + "from __future__ import division\n", + "from scipy.misc import derivative\n", + "from numpy import arange,sin,cos,sqrt\n", + "x = 0.45#\n", + "def f(x):\n", + " F = sin(x)\n", + " return F\n", + "def df(x,h):\n", + " DF = (f(x+h) - f(x))/h\n", + " return DF\n", + "dfactual = cos(x)#\n", + "h = arange(0.01,0.005+0.04,0.005)\n", + "n = len(h)#\n", + "print 'h\\ty\\t\\terr'\n", + "for hh in h:\n", + " y = df(x,hh)\n", + " err = y - dfactual \n", + " print hh,'\\t',y,'\\t',err\n", + "\n", + "#using 16 significant digits so the bound for roundoff error is 0.5*10**(-16)\n", + "e = 0.5*10**(-16)\n", + "M2 = max(sin(x+h))#\n", + "hopt = 2*sqrt(e/M2)#\n", + "print '\\nM2 = %.3e & hopt = %.3e'%(hopt,M2)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 11_04 Pg No. 354" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + " d2fexact = -0.732 \n", + " err = -6.097e-06 \n", + " y = -0.732 \n" + ] + } + ], + "source": [ + "from __future__ import division\n", + "from numpy import arange,sin,cos,sqrt\n", + "\n", + "#Approximate second derivative\n", + "\n", + "x = 0.75#\n", + "h = 0.01#\n", + "def f(x):\n", + " F = cos(x)\n", + " return F\n", + "def d2f(x,h):\n", + " D2F = ( f(x+h) - 2*f(x) + f(x-h) )/h**2\n", + " return D2F\n", + "y = d2f(0.75,0.01)#\n", + "d2fexact = -cos(0.75)\n", + "err = d2fexact - y #\n", + "print ' d2fexact = %.3f \\n err = %.3e \\n y = %.3f '%(d2fexact,err,y)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 11_05 Pg No. 358" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "v(5) = 4.25\n", + "v(7) = 5.5\n", + "v(9) = 6.75\n" + ] + } + ], + "source": [ + "#Differentiation of tabulated data\n", + "\n", + "T = range(5,10)\n", + "s = [10, 14.5, 19.5, 25.5 , 32 ]\n", + "h = T[1]-T[0]\n", + "n = len(T)\n", + "def v(t):\n", + " if t in T and T.index(t) == 0:\n", + " V = (-3*s[T.index(t) ] + 4*s[T.index(t+h)] - s[T.index(t+2*h) ] ) /(2*h) #Three point forward difference formula\n", + " elif t in T and T.index(t) == n-1:\n", + " V = ( 3*s[T.index(t)] - 4*s[ T.index(t-h)] + s[T.index(t-2*h) ] )/(2*h) #Backward difference formula\n", + " else:\n", + " V = ( s[T.index(t+h)] - s[T.index(t-h)] )/(2*h) #Central difference formula\n", + " return V\n", + "\n", + "v_5 = v(5)\n", + "v_7 = v(7)\n", + "v_9 = v(9)\n", + "\n", + "print 'v(5) = ',v_5\n", + "print 'v(7) = ',v_7\n", + "print 'v(9) = ',v_9\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 11_06 Pg No. 359" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "a(7) = 0.5\n" + ] + } + ], + "source": [ + "T = range(5,10)\n", + "s = [10, 14.5 , 19.5 , 25.5, 32 ]\n", + "h = T[1]-T[0]\n", + "def a(t):\n", + " A = ( s[T.index(t+h) ] - 2*s[ T.index(t) ] + s[T.index(t-h) ] )/h**2\n", + " return A\n", + "a_7 = a(6)\n", + "\n", + "print 'a(7) = ',a_7" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 11_7 Pg No. 359" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "solving the above equations we get \n", + "\n", + " y1 = y(0.25) = -0.117562 \n", + " y2 = y(0.5) = -0.168475 \n", + " y3 = y(0.75) = -0.139088 \n", + " \n" + ] + } + ], + "source": [ + "from numpy import mat,divide,linalg\n", + "h = 0.25 #\n", + "#y''(x) = e**(x**2)\n", + "#y(0) = 0 , y(1) = 0\n", + "# y''(x) = y(x+h) - 2*y(x) + y(x-h)/h**2 = e**(x**2)\n", + "#( y(x + 0.25) - 2*y(x) + y(x-0.25) )/0.0625 = e**(x**2)\n", + "#y(x+0.25) - 2*y(x) + y(x - 0.25) = 0.0624*e**(x**2)\n", + "#y(0.5) - 2*y(0.25) + y(0) = 0.0665\n", + "#y(0.75) - 2*y(0.5) + y(0.25) = 0.0803\n", + "#y(1) - 2*y(0.75) + y(0.5) = 0.1097\n", + "#given y(0) = y(1) = 0\n", + "#\n", + "#0 + y2 - 2y1 = 0.06665\n", + "#y3 - 2*y2 + y1 = 0.0803\n", + "#-2*y3 + y2 + 0 = 0.1097\n", + "#Therefore\n", + "A = mat([[0, 1, -2],[1, -2, 1],[-2, 1, 0 ]])\n", + "B = mat([[ 0.06665],[0.0803] ,[0.1097 ]])\n", + "C = linalg.solve(A,B)\n", + "print 'solving the above equations we get \\n\\n y1 = y(0.25) = %f \\n y2 = y(0.5) = %f \\n y3 = y(0.75) = %f \\n '%(C[2],C[1],C[0])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 11_08 Pg No. 362" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "richardsons technique -\n", + "\n", + "df(0.5) = 1.64850499031\n", + "D(rh) = D(0.25) = 1.71828182846\n", + "D(0.5) = 1.66594919985\n", + "Exact df(0.5) = 1.93757920531\n", + "The result by richardsons technique is much better than other results\n", + "for r = 2\n", + "df(x) = 1.64518270284\n", + "D(rh) = 1.93757920531\n" + ] + } + ], + "source": [ + "from numpy import arange,exp\n", + "from scipy.misc import derivative\n", + "\n", + "x = arange(-0.5,0.25+1.5,0.25)\n", + "h = 0.5 ;\n", + "r = 1.0/2 ;\n", + "\n", + "def f(x):\n", + " F = exp(x)\n", + " return F\n", + "def D(x,h):\n", + " D = (f(x + h) - f(x-h) )/(2*h) \n", + " return D\n", + "def df(x,h,r):\n", + " Df = (D(x,r*h) - r**2*D(x,h))/(1-r**2)\n", + " return Df\n", + "\n", + "df_05 = df(0.5,0.5,1.0/2);\n", + "print 'richardsons technique -\\n\\ndf(0.5) = ',df_05\n", + "print 'D(rh) = D(0.25) = ',D(0.5,0.5)\n", + "print 'D(0.5) = ',D(0.5,0.25)\n", + "\n", + "dfexact = derivative(f,0.5)\n", + "print 'Exact df(0.5) = ',dfexact\n", + "print 'The result by richardsons technique is much better than other results'\n", + "\n", + "#r = 2\n", + "print 'for r = 2'\n", + "print 'df(x) = ',df(0.5,0.5,2)\n", + "print 'D(rh) = ',D(0.5,2*0.5)" + ] + } + ], + "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 +} diff --git a/Numerical_Methods_by_E._Balaguruswamy/chapter12.ipynb b/Numerical_Methods_by_E._Balaguruswamy/chapter12.ipynb new file mode 100644 index 00000000..9bd54b00 --- /dev/null +++ b/Numerical_Methods_by_E._Balaguruswamy/chapter12.ipynb @@ -0,0 +1,518 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Chapter 12 - Numerical integration" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 12_01 Pg No. 373" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "case(a)\n", + "It = 5\n", + "Ett = 1.0\n", + "Iexact = 4.75\n", + "True error = 0.25 Here Error bound is an overestimate of true error\n", + "case(b)\n", + "It = 1.59375\n", + "Ett = 0.09375\n", + "Iexact = 1.515625\n", + "True error = 0.078125 Here Error bound is an overestimate of true error\n" + ] + } + ], + "source": [ + "from scipy.misc import derivative\n", + "from mpmath import quad\n", + "#Trapezoidal Rule\n", + "\n", + "def f(x):\n", + " F = x**3 + 1\n", + " return F\n", + "\n", + "#case(a)\n", + "a = 1#\n", + "b = 2 #\n", + "h = b - a #\n", + "It = (b-a)*(f(a)+f(b))/2\n", + "d2f = derivative(f,2,n=2)\n", + "Ett = h**3*d2f/12.0\n", + "Iexact = quad(f,[1,2])\n", + "Trueerror = It - Iexact\n", + "print 'case(a)'\n", + "print 'It = ',It\n", + "print 'Ett = ',Ett\n", + "print 'Iexact = ',Iexact\n", + "print 'True error = ',Trueerror,\n", + "print 'Here Error bound is an overestimate of true error'\n", + "\n", + "#case(b)\n", + "a = 1#\n", + "b = 1.5 #\n", + "h = b - a #\n", + "It = (b-a)*(f(a)+f(b))/2\n", + "Ett = h**3*derivative(f,1.5,n=2)/12\n", + "Iexact = quad(f,[1,1.5])\n", + "Trueerror = It - Iexact\n", + "print 'case(b)'\n", + "print 'It = ',It\n", + "print 'Ett = ',Ett\n", + "print 'Iexact = ',Iexact\n", + "print 'True error = ',Trueerror,\n", + "print 'Here Error bound is an overestimate of true error'\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 12_02 Pg No. 376" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "intergral for case(a),Ia = 2.54308063482\n", + "intergral for case(b),Ib = 0.0\n", + "exact integral,Iexact = 2.3504\n", + "n = 4 case is better than n = 2 case\n" + ] + } + ], + "source": [ + "from math import exp\n", + "#Tapezoidal rule\n", + "\n", + "\n", + "def f(x):\n", + " F = exp(x)\n", + " return F\n", + "a = -1 #\n", + "b = 1 #\n", + "\n", + "#case(a)\n", + "n = 2\n", + "h = (b-a)/n \n", + "I = 0\n", + "for i in range(1,n+1):\n", + " I = I + f(a+(i-1)*h)+f(a+i*h)\n", + "\n", + "I = h*I/2 #\n", + "print 'intergral for case(a),Ia = ',I\n", + "\n", + "#case(b)\n", + "n = 4\n", + "h = (b-a)/n \n", + "I = 0\n", + "for i in range(1,n+1):\n", + " I = I + f(a+(i-1)*h)+f(a+i*h)#\n", + "\n", + "I = h*I/2 #\n", + "Iexact = 2.35040\n", + "print 'intergral for case(b),Ib = ',I\n", + "print 'exact integral,Iexact = ',Iexact\n", + "print 'n = 4 case is better than n = 2 case'" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 12_03 Pg No. 381" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "h = 1\n", + "Integral for case(a) , Is1 = 2.36205375654\n", + "h = 0.785398163397\n", + "Integral for case(b),Is1 = 1.14238405466\n" + ] + } + ], + "source": [ + "from math import pi,sin,exp,sqrt\n", + "#Simpon's 1/3 rule\n", + "\n", + "#case(a)\n", + "def f(x):\n", + " F = exp(x)\n", + " return F\n", + "a = -1#\n", + "b = 1#\n", + "h = (b-a)/2 \n", + "x1 = a+h\n", + "Is1 = h*( f(a) + f(b) + 4*f(x1) )/3 \n", + "print 'h = ',h\n", + "print 'Integral for case(a) , Is1 = ',Is1\n", + "\n", + "#case(b)\n", + "def f(x):\n", + " F = sqrt(sin(x))\n", + " return F\n", + "a = 0\n", + "b = pi/2\n", + "h = (b-a)/2 \n", + "x1 = a+h\n", + "Is1 = h*( f(a) + f(b) + 4*f(x1) )/3\n", + "print 'h = ',h\n", + "print 'Integral for case(b),Is1 = ',Is1" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 12_04 Pg No.382" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "h = 0.392699081699\n", + "Integral value for n = 4 is 1.17822754435\n", + "h = 0.261799387799\n", + "Integral value for n = 6 is 1.18728120346\n" + ] + } + ], + "source": [ + "from math import pi,sin,exp,sqrt\n", + "#Simpon's 1/3 rule\n", + "\n", + "def f(x):\n", + " F = sqrt( sin(x) )\n", + " return F\n", + "x0 = 0 #\n", + "xa = pi/2 #\n", + "\n", + "#case(a) n = 4\n", + "n = 4 #\n", + "h = ( xa-x0 )/n\n", + "I = 0 \n", + "for i in range(1,n/2+1):\n", + " I = I + f(x0 + (2*i-2)*h) + 4*f(x0 + (2*i-1)*h) + f(x0 + 2*i*h) #\n", + "I = h*I/3\n", + "print 'h = ',h\n", + "print 'Integral value for n = 4 is',I \n", + "\n", + "#case(b) n = 6\n", + "n = 6\n", + "h = ( xa-x0 )/n\n", + "I = 0 \n", + "for i in range(1,n/2+1):\n", + " I = I + f(x0 + (2*i-2)*h) + 4*f(x0 + (2*i-1)*h) + f(x0 + 2*i*h) #\n", + "\n", + "I = h*I/3\n", + "print 'h = ',h\n", + "print 'Integral value for n = 6 is',I" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 12_05 Pg No. 386" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Integral of x**3 +1 from 1 to 2 : 0\n", + "Integral of sqrt(sin(x)) from 0 to pi/2: 1.16104132669\n" + ] + } + ], + "source": [ + "from math import pi,sin,exp,sqrt\n", + "\n", + "#Simpson's 3/8 rule\n", + "\n", + "#case(a)\n", + "def f(x):\n", + " F = x**3 + 1\n", + " return F\n", + "a = 1 #\n", + "b = 2 #\n", + "h = (b-a)/3 \n", + "x1 = a + h \n", + "x2 = a + 2*h\n", + "Is2 = 3*h*( f(a) + 3*f(x1) + 3*f(x2) + f(b) )/8 #\n", + "print 'Integral of x**3 +1 from 1 to 2 : ',Is2\n", + "#case(b)\n", + "def f(x):\n", + " F = sqrt( sin(x) )\n", + " return F\n", + "a = 0 #\n", + "b = pi/2 #\n", + "h = (b-a)/3 \n", + "x1 = a + h \n", + "x2 = a + 2*h\n", + "Is2 = 3*h*( f(a) + 3*f(x1) + 3*f(x2) + f(b) )/8 #\n", + "print 'Integral of sqrt(sin(x)) from 0 to pi/2:',Is2" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 12_06 Pg No. 387" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Ib = 1.18061711033\n" + ] + } + ], + "source": [ + "from math import sin,sqrt,pi\n", + "#Booles's Five-Point formula\n", + "\n", + "def f(x):\n", + " F = sqrt( sin(x) )\n", + " return F\n", + "x0 = 0#\n", + "xb = pi/2 #\n", + "n = 4 #\n", + "h = (xb - x0)/n\n", + "Ib = 2*h*(7*f(x0) + 32*f(x0+h) + 12*f(x0 + 2*h) + 32*f(x0+3*h) + 7*f(x0+4*h))/45#\n", + "print 'Ib = ',Ib" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 12_07 Pg No. 391" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "R(0,0) = 0.75\n", + "\n", + "R(1,0) = 0.333333 \n", + "\n", + "\n", + "R(2,0) = 0.342857 \n", + "\n", + "\n", + "R(1,1) = 0.194444 \n", + "\n", + "\n", + "R(2,1) = 0.346032 \n", + "\n", + "\n", + "R(2,2) = 0.356138 \n", + "\n" + ] + } + ], + "source": [ + "from __future__ import division\n", + "from numpy import zeros\n", + "def f(x):\n", + " F = 1.0/x\n", + " return F\n", + "#since we can't have (0,0) element in matrix we start with (1,1)\n", + "a = 1 ;\n", + "b = 2 ;\n", + "h = b-a ;\n", + "R=zeros([3,3])\n", + "R[0,0] = h*(f(a)+f(b))/2 \n", + "print 'R(0,0) = ',R[0,0]\n", + "\n", + "h=[h]\n", + "for i in range(2,4):\n", + " h.append((b-a)/2**(i-1))\n", + " s = 0\n", + " for k in range(1,2**(i-2)+1):\n", + " s = s + f(a + (2*k - 1)*h[i-1]);\n", + " \n", + " R[i-1,0] = R[i-1,0]/2 + h[i-1]*s;\n", + " print '\\nR(%i,0) = %f \\n'%(i-1,R[i-1,0])\n", + "\n", + "for j in range(2,4):\n", + " for i in range(j,4):\n", + " R[i-1,j-1] = (4**(j-1)*R[i-1,j-2] - R[i-2,j-2])/(4**(j-1)-1);\n", + " print '\\nR(%i,%i) = %f \\n'%(i-1,j-1,R[i-1,j-1])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 12_08 Pg No. 397" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "x1 = -0.57735026919\n", + "x2 = 0.57735026919\n", + "I = 2.34269608791\n" + ] + } + ], + "source": [ + "from math import sqrt,exp \n", + "#Two Point Gauss -Legedre formula\n", + "\n", + "def f(x):\n", + " F = exp(x)\n", + " return F\n", + "x1 = -1/sqrt(3)\n", + "x2 = 1/sqrt(3)\n", + "I = f(x1) + f(x2)\n", + "print 'x1 = ',x1\n", + "print 'x2 = ',x2\n", + "print 'I = ',I, " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 12_09 Pg No. 398" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + " g(z) = exp(-(2.000000*z + 0.000000)/2) \n", + " C = 2.000000 \n", + " Ig = 4.685392 \n", + "\n" + ] + } + ], + "source": [ + "from math import sqrt,exp\n", + "#Gaussian two point formula\n", + "\n", + "a = -2 #\n", + "b = 2 #\n", + "def f(x):\n", + " F = exp(-x/2)\n", + " return F\n", + "A = (b-a)/2 \n", + "B = (a+b)/2\n", + "C = (b-a)/2\n", + "def g(z):\n", + " G = exp(-1*(A*z+B)/2)\n", + " return G\n", + "w1 = 1 #\n", + "w2 = 1 #\n", + "z1 = -1/sqrt(3)\n", + "z2 = 1/sqrt(3)\n", + "Ig = C*( w1*g(z1) + w2*g(z2) )\n", + "print ' g(z) = exp(-(%f*z + %f)/2) \\n C = %f \\n Ig = %f \\n'%(A,B,C,Ig)" + ] + } + ], + "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 +} diff --git a/Numerical_Methods_by_E._Balaguruswamy/chapter13.ipynb b/Numerical_Methods_by_E._Balaguruswamy/chapter13.ipynb new file mode 100644 index 00000000..2aa34503 --- /dev/null +++ b/Numerical_Methods_by_E._Balaguruswamy/chapter13.ipynb @@ -0,0 +1,911 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Chapter 13 - Numerical solution of ordinary differential equations" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 13_01 Pg No. 414" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + " y(0.25) = 5.3125\n", + " y(0.5) = 7.40685424949\n" + ] + } + ], + "source": [ + "from __future__ import division\n", + "from math import sqrt\n", + "#Taylor method\n", + "\n", + "def f(x,y):\n", + " F = x**2 + y**2\n", + " return F\n", + "def d2y(x,y):\n", + " D2Y = 2*x + 2*y*f(x,y)\n", + " return D2Y\n", + "def d3y(x,y):\n", + " D3Y = 2 + 2*y*d2y(x,y) + 2*f(x,y)**2\n", + " return D3Y\n", + "def y(x):\n", + " Y = 1 + f(0,1)*x + d2y(0,1)*x**2/2 + d3y(0,1)*sqrt(x)\n", + " return Y\n", + "print ' y(0.25) = ',y(0.25)\n", + "print ' y(0.5) = ',y(0.5)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 13_02 Pg No. 415" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + " Iteration-1\n", + "\n", + " dy(0) = 0.000000\n", + " d2y(0) = 0.000000\n", + " d3y(0) = 2.000000\n", + " d4y(0) = 0.000000\n", + " \n", + "y(0) = 0.002667\n", + "\n", + "\n", + " Iteration-2\n", + "\n", + " dy(0) = 0.040007\n", + " d2y(0) = 0.400213\n", + " d3y(0) = 2.005336\n", + " d4y(0) = 0.106763\n", + " \n", + "y(0) = 0.021353\n", + "\n", + "\n" + ] + } + ], + "source": [ + "from __future__ import division\n", + "from scipy.misc import factorial\n", + "#Recursive Taylor Method\n", + "\n", + "def f(x,y):\n", + " F = x**2 + y**2\n", + " return F\n", + "def d2y(x,y):\n", + " D2Y = 2*x + 2*y*f(x,y)\n", + " return D2Y\n", + "def d3y(x,y):\n", + " D3Y = 2 + 2*y*d2y(x,y) + 2*f(x,y)**2\n", + " return D3Y\n", + "def d4y(x,y):\n", + " D4Y = 6*f(x,y)*d2y(x,y) + 2*y*d3y(x,y)\n", + " return D4Y\n", + "h = 0.2 #\n", + "def y(x,y):\n", + " Y = y + f(x,y)*h + d2y(x,y)*h**2/2 + d3y(x,y)*h**3/6 + d4y(x,y)*h**4/factorial(4)\n", + " return Y\n", + "x0 = 0#\n", + "y0 = 0 #\n", + "y_=[]\n", + "for i in range(1,3):\n", + " y_.append(y(x0,y0))\n", + " print ' Iteration-%d\\n\\n dy(0) = %f\\n d2y(0) = %f\\n d3y(0) = %f\\n d4y(0) = %f\\n '%(i,f(x0,y0),d2y(x0,y0),d3y(x0,y0),d4y(x0,y0)) \n", + " x0 = x0 + i*h\n", + " y0 = y_[i-1]\n", + " print 'y(0) = %f\\n\\n'%(y_[i-1])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 13_03 Pg No. 417" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "for dy(x) = x**2 + y**2 the results are :\n", + "y(0.1) = 0.000333334920635\n", + "y(0.2) = 0.00266686984127\n", + "y(1) = 0.349206349206\n", + "for dy(x) = x*e**y the results are \n", + "y(0.1) = 0.0050125208594\n", + "y(0.2) = 0.0202013400268\n", + "y(1) = 0.6487212707\n" + ] + } + ], + "source": [ + "from __future__ import division\n", + "from math import exp \n", + "#Picard's Method\n", + "\n", + "#y'(x) = x**2 + y**2,y(0) = 0\n", + "#y(1) = y0 + integral(x**2 + y0**2,x0,x)\n", + "#y(1) = x**3/3\n", + "#y(2) = 0 + integral(xY2 + y1**2,x0,x)\n", + "# = integral(x**2 + x**6/9,0,x) = x**3/3 + x**7/63\n", + "#therefore y(x) = x**3 /3 + x**7/63\n", + "def y(x):\n", + " Y = x**3/3 + x**7/63 \n", + " return Y\n", + "print 'for dy(x) = x**2 + y**2 the results are :'\n", + "print 'y(0.1) = ',y(0.1)\n", + "print 'y(0.2) = ',y(0.2)\n", + "print 'y(1) = ',y(1)\n", + "\n", + "#y'(x) = x*e**y, y(0) = 0\n", + "#y0 = 0 , x0 = 0\n", + "#Y(1) = 0 + integral(x*e**0,0,x) = x**2/2\n", + "#y(2) = 0 + integral( x*e**( x**2/2 ) ,0,x) = e**(x**2/2)-1\n", + "#therefore y(x) = e**(x**2/2) - 1\n", + "def y(x):\n", + " Y = exp(x**2/2) - 1 \n", + " return Y\n", + "print 'for dy(x) = x*e**y the results are ' \n", + "print 'y(0.1) = ',y(0.1)\n", + "print 'y(0.2) = ',y(0.2)\n", + "print 'y(1) = ',y(1)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 13_04 Pg No. 417" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "for h = 0.500000\n", + "\n", + "y(1.500000) = 4.000000\n", + "\n", + "y(2.000000) = 7.875000\n", + "\n", + "\n", + "for h = 0.250000\n", + "\n", + "y(1.250000) = 3.000000\n", + "\n", + "y(1.500000) = 4.421875\n", + "\n", + "y(1.750000) = 6.359375\n", + "\n", + "y(2.000000) = 8.906250\n", + "\n" + ] + } + ], + "source": [ + "from __future__ import division\n", + "\n", + "#Euler's Method\n", + "\n", + "def dy(x):\n", + " DY = 3*x**2 + 1\n", + " return DY\n", + "x0 = 1\n", + "y = [0,2] #\n", + "#h = 0.5\n", + "h = 0.5\n", + "print 'for h = %f\\n'%(h)\n", + "for i in range(2,4):\n", + " y.append(y[i-1] + h*dy(x0+(i-2)*h))\n", + " print 'y(%f) = %f\\n'%(x0+(i-1)*h,y[i])\n", + "\n", + "#h = 0.25\n", + "h = 0.25\n", + "print '\\nfor h = %f\\n'%(h)\n", + "y = [0,2] #\n", + "for i in range(2,6):\n", + " y.append(y[i-1] + h*dy(x0+(i-2)*h))\n", + " print 'y(%f) = %f\\n'%(x0+(i-1)*h,y[i])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 13_05 Pg No. 422" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\n", + " Step 1 \n", + " x(1) = 1.500000\n", + " y(1.500000) = 4.000000\n", + "\n", + "\n", + " Et(1) = 0.875000\n", + "\n", + "\n", + " Step 2 \n", + " x(2) = 2.000000\n", + " y(2.000000) = 7.875000\n", + "\n", + "\n", + " Et(2) = 1.250000\n", + "\n", + "x Est y true y Et Globalerr\n", + "1.5 \t4.0 \t4.875 \t0.875 \t0.875\n", + "2.0 \t7.875 \t10.0 \t1.25 \t2.125\n" + ] + } + ], + "source": [ + "from numpy import nditer\n", + "from __future__ import division\n", + "#Error estimation in Euler's Method\n", + "\n", + "def dy(x):\n", + " DY = 3*x**2 + 1\n", + " return DY\n", + "def d2y(x):\n", + " D2Y = 6*x\n", + " return D2Y\n", + "def d3y(x):\n", + " D3Y = 6\n", + " return D3Y\n", + "def exacty(x):\n", + " Exacty = x**3 + x\n", + " return Exacty\n", + "x0 = 1\n", + "y = 2\n", + "h = 0.5\n", + "X=[];ESTY=[];TRUEY=[];ET=[];GERR=[];\n", + "for i in range(2,4):\n", + " x = x0 + (i-1)*h\n", + " y = y + h*dy(x0+(i-2)*h)\n", + " print '\\n Step %d \\n x(%d) = %f\\n y(%f) = %f\\n'%(i-1,i-1,x,x,y)\n", + " Et = d2y(x0+(i-2)*h)*h**2/2 + d3y(x0+(i-2)*h)*h**3/6\n", + " print '\\n Et(%d) = %f\\n'%(i-1,Et)\n", + " truey = exacty(x0+(i-1)*h)\n", + " gerr = truey - y\n", + " X.append(x)\n", + " ESTY.append(y)\n", + " TRUEY.append(truey)\n", + " ET.append(Et)\n", + " GERR.append(gerr)\n", + "\n", + "#table = [x y(2:3) truey Et gerr]\n", + "print 'x Est y true y Et Globalerr' \n", + "for a,b,c,d,e in nditer([X,ESTY,TRUEY,ET,GERR]):\n", + " print a,'\\t',b,'\\t',c,'\\t',d,'\\t',e" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 13_06 Pg No. 427" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "EULERS METHOD\n", + "y(1.250000) = 3.000000 \n", + " \n", + "y(1.500000) = 4.200000 \n", + " \n", + "y(1.750000) = 5.600000 \n", + " \n", + "y(2.000000) = 7.200000 \n", + " \n", + "HEUNS METHOD\n", + "\n", + "Iteration 1 \n", + " m1 = 4.000000\n", + " ye(1.250000) = 3.000000 \n", + " m2 = 4.800000 \n", + " y(1.250000) = 3.100000 \n", + "\n", + "\n", + "Iteration 2 \n", + " m1 = 4.960000\n", + " ye(1.500000) = 4.340000 \n", + " m2 = 5.786667 \n", + " y(1.500000) = 4.443333 \n", + "\n", + "\n", + "Iteration 3 \n", + " m1 = 5.924444\n", + " ye(1.750000) = 5.924444 \n", + " m2 = 6.770794 \n", + " y(1.750000) = 6.030238 \n", + "\n", + "\n", + "Iteration 4 \n", + " m1 = 6.891701\n", + " ye(2.000000) = 7.753163 \n", + " m2 = 7.753163 \n", + " y(2.000000) = 7.860846 \n", + "\n", + "x Eulers Heuns \t\tAnalytical\n", + "0.0 \t0.0 \t0.000000 \t0.0\n", + "1.0 \t2.0 \t2.000000 \t2.0\n", + "1.25 \t3.0 \t3.100000 \t2.2360679775\n", + "1.5 \t4.2 \t4.443333 \t2.44948974278\n", + "1.75 \t5.6 \t6.030238 \t2.64575131106\n", + "2.0 \t7.2 \t7.860846 \t2.82842712475\n" + ] + } + ], + "source": [ + "from numpy import nditer,sqrt\n", + "from __future__ import division\n", + "\n", + "#Heun's Method\n", + "\n", + "def f(x,y):\n", + " F = 2*y/x \n", + " return F\n", + "def exacty(x):\n", + " Exacty = 2*sqrt(x)\n", + " return Exacty\n", + "x = [0,1] #\n", + "y = [0,2] #\n", + "h = 0.25 #\n", + "#Euler's Method\n", + "print 'EULERS METHOD'\n", + "for i in range(2,6):\n", + " x.append(x[i-1] + h )\n", + " y.append(y[i-1] + h*f(x[i-1],y[i-1]))\n", + " print 'y(%f) = %f \\n '%(x[i],y[i])\n", + "\n", + "eulery = y\n", + "#Heun's Method\n", + "print 'HEUNS METHOD'\n", + "ye=[0,0]\n", + "y = [0,2] #\n", + "for i in range(2,6):\n", + " m1 = f(x[i-1],y[i-1]) #\n", + " ye.append(y[i-1] + h*f(x[(i-1)],y[(i-1)]))\n", + " m2 = f(x[(i)],ye[(i)]) \n", + " y.append(y[(i-1)] + h*(m1 + m2)/2)\n", + " print '\\nIteration %d \\n m1 = %f\\n ye(%f) = %f \\n m2 = %f \\n y(%f) = %f \\n'%(i-1,m1,x[(i)],ye[(i)],m2,x[(i)],y[(i)]) \n", + " \n", + " \n", + "truey = exacty(x) \n", + "print 'x Eulers Heuns \\t\\tAnalytical' \n", + "for a,b,c,d in nditer([x,eulery,y,truey]):\n", + " print a,'\\t',b,'\\t%.6f'%c,'\\t',d\n", + " " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 13_07 Pg NO. 433" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "y(1.250000) = 3.111111 \n", + " \n", + "y(1.500000) = 4.468687 \n", + " \n" + ] + } + ], + "source": [ + "from __future__ import division\n", + "\n", + "#Polygon Method\n", + "\n", + "def f(x,y):\n", + " F = 2*y/x\n", + " return F\n", + "x=[1]\n", + "y=[2]\n", + "h = 0.25 #\n", + "for i in range(1,3):\n", + " x.append(x[(i-1)] + h )\n", + " y.append(y[(i-1)] + h*f( x[(i-1)]+ h/2 , y[(i-1)] + h*f( x[(i-1)] , y[(i-1)] )/2 ))\n", + " print 'y(%f) = %f \\n '%(x[(i)],y[(i)])\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 13_08 Pg No. 439" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\n", + "Iteration - 1\n", + " m1 = 0.000000\n", + " m2 = 0.010000 \n", + " m3 = 0.010001 \n", + " m4 = 0.040004 \n", + " y(0.200000) = 0.002667 \n", + "\n", + "\n", + "Iteration - 2\n", + " m1 = 0.040007\n", + " m2 = 0.090044 \n", + " m3 = 0.090136 \n", + " m4 = 0.160428 \n", + " y(0.400000) = 0.021360 \n", + "\n" + ] + } + ], + "source": [ + "from __future__ import division\n", + "\n", + "#Classical Runge Kutta Method\n", + "\n", + "def f(x,y):\n", + " F = x**2 + y**2\n", + " return F\n", + "h = 0.2\n", + "x=[0]\n", + "y=[0]\n", + "\n", + "for i in range(0,2):\n", + " m1 = f( x[(i)] , y[(i)] ) #\n", + " m2 = f( x[i] + h/2 , y[(i)] + m1*h/2 ) #\n", + " m3 = f( x[(i)] + h/2 , y[(i)] + m2*h/2 ) #\n", + " m4 = f( x[(i)] + h , y[(i)] + m3*h ) #\n", + " x.append(x[(i)] + h )\n", + " y.append(y[(i)] + (m1 + 2*m2 + 2*m3 + m4)*h/6 )\n", + " \n", + " print '\\nIteration - %d\\n m1 = %f\\n m2 = %f \\n m3 = %f \\n m4 = %f \\n y(%f) = %f \\n'%(i+1,m1,m2,m3,m4,x[(i+1)],y[(i+1)])\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 13_09 Pg No. 444" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "for four decimal places\n", + "h = 0.0143668085653\n", + "for six decimal places\n", + "h = 0.00454318377739\n", + "Note-We can use h = 0.01 for four decimal places and h = 0.004 for six decimal places\n" + ] + } + ], + "source": [ + "from __future__ import division\n", + "\n", + "#Optimum step size\n", + "\n", + "x = 0.8 #\n", + "h1 = 0.05 #\n", + "y1 = 5.8410870 #\n", + "h2 = 0.025 #\n", + "y2 = 5.8479637 #\n", + "\n", + "#d = 4\n", + "h = ((h1**4 - h2**4)*10**(-4)/(2*(y2 - y1)))**(1/4)\n", + "print 'for four decimal places'\n", + "print 'h = ',h \n", + "\n", + "#d = 6\n", + "h = ((h1**4 - h2**4)*10**(-6)/(2*(y2 - y1)))**(1/4)\n", + "print 'for six decimal places'\n", + "print 'h = ' ,h\n", + "print 'Note-We can use h = 0.01 for four decimal places and h = 0.004 for six decimal places'" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 13_10 Pg NO. 446" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "yp4 = 8.00914285714\n", + "fp4 = 8.00914285714 yc4 = 8.00266666667\n", + "f4 = 8.00266666667\n", + "yc4 = 8.00212698413\n", + "Note- the exact solution is y(2) = 8\n" + ] + } + ], + "source": [ + "from __future__ import division\n", + "\n", + "#Milne-Simpson Predictor-Corrector method\n", + "\n", + "def f(x,y):\n", + " F = 2*y/x\n", + " return F\n", + "x0 = 1 #\n", + "y0 = 2 #\n", + "h = 0.25 #\n", + "#Assuming y1 ,y2 and y3(required for milne-simpson formula) are estimated using Fourth- order Runge kutta method\n", + "x1 = x0 + h \n", + "y1 = 3.13 #\n", + "x2 = x1 + h \n", + "y2 = 4.5 #\n", + "x3 = x2 + h\n", + "y3 = 6.13 #\n", + "#Milne Predictor formula\n", + "yp4 = y0 + 4*h*(2*f(x1,y1) - f(x2,y2) + 2*f(x3,y3))/3\n", + "x4 = x3 + h \n", + "fp4 = f(x4,yp4) #\n", + "print 'yp4 = ',yp4\n", + "print 'fp4 = ',fp4, \n", + "#Simpson Corrector formula\n", + "yc4 = y2 + h*( f(x2,y2) + 4*f(x3,y3) + fp4)/3 \n", + "f4 = f(x4,yc4)\n", + "print 'yc4 = ',yc4\n", + "print 'f4 = ',f4 \n", + "\n", + "yc4 = y2 + h*( f(x2,y2) + 4*f(x3,y3) + f4)/3 \n", + "print 'yc4 = ',yc4 \n", + "print 'Note- the exact solution is y(2) = 8'" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 13_11 Pg NO. 446" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Adams Predictor formula\n", + "yp4 = 8.01135714286\n", + "fp4 = 8.01135714286\n", + "\n", + "Adams Corrector formula\n", + "yc4 = 8.00727901786\n", + "f4 = 8.00727901786\n", + "\n", + "refined-yc4 = 8.00689669364\n" + ] + } + ], + "source": [ + "from __future__ import division\n", + "\n", + "#Adams-Bashforth-Moulton Method\n", + "\n", + "def f(x,y):\n", + " F = 2*y/x\n", + " return F\n", + "x0 = 1 #\n", + "y0 = 2 #\n", + "h = 0.25 #\n", + "x1 = x0 + h \n", + "y1 = 3.13 #\n", + "x2 = x1 + h \n", + "y2 = 4.5 #\n", + "x3 = x2 + h\n", + "y3 = 6.13 #\n", + "#Adams Predictor formula\n", + "yp4 = y3 + h*(55*f(x3,y3) - 59*f(x2,y2) + 37*f(x1,y1) - 9*f(x0,y0))/24\n", + "x4 = x3 + h \n", + "fp4 = f(x4,yp4) \n", + "print 'Adams Predictor formula'\n", + "print 'yp4 = ',yp4 \n", + "print 'fp4 = ',fp4 \n", + "#Adams Corrector formula\n", + "yc4 = y3 + h*( f(x1,y1) - 5*f(x2,y2) + 19*f(x3,y3) + 9*fp4 )/24 \n", + "f4 = f(x4,yc4)\n", + "print '\\nAdams Corrector formula'\n", + "print 'yc4 = ',yc4 \n", + "print 'f4 = ',f4 \n", + "\n", + "yc4 = y3 + h*( f(x1,y1) - 5*f(x2,y2) + 19*f(x3,y3) + 9*f4 )/24 \n", + "print '\\nrefined-yc4 = ',yc4" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 13_12 Pg No. 453" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\n", + " y4p = 0.557351\n", + " f4p = -0.310640 \n", + " y4c = 0.555396 \n", + " Modifier Etc = 0.000067 \n", + " Modified y4c = 0.555464 \n", + "\n", + "\n", + " y5p = 0.500837\n", + " f5p = -0.250837 \n", + " y5c = 0.499959 \n", + " Modifier Etc = 0.000030 \n", + " Modified y5c = 0.499989 \n", + "\n", + "error = 1.11332736336e-05\n" + ] + } + ], + "source": [ + "from __future__ import division\n", + "\n", + "#Milne-Simpson Method using modifier\n", + "\n", + "def f(y):\n", + " F = -y**2\n", + " return F\n", + "x = [ 1 , 1.2 , 1.4 , 1.6 ] \n", + "y = [ 1 , 0.8333333 , 0.7142857 , 0.625 ] \n", + "h = 0.2 #\n", + "\n", + "for i in range(0,2):\n", + " yp = y[(i)] + 4*h*( 2*f( y[(i+1)] ) - f( y[(i+2)] ) + 2*f( y[(i+3)] ) )/3\n", + " fp = f(yp) #\n", + " yc = y[( i+2)] + h*(f( y[(i+2)] ) + 4*f( y[(i+3)] ) + fp )/3 \n", + " Etc = -(yc - yp)/29\n", + " y.append(yc + Etc)\n", + " print '\\n y%dp = %f\\n f%dp = %f \\n y%dc = %f \\n Modifier Etc = %f \\n Modified y%dc = %f \\n'%(i+4,yp,i+4,fp,i+4,yc,Etc,i+4,y[(i+4)])\n", + "exactanswer = 0.5 #\n", + "err = exactanswer - y[5] #\n", + "print 'error = ',err" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 13_13 Pg No. 455" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "m1(1) = 0.000000\n", + " m1(2) = 1.000000\n", + " m2(1) = 0.200000\n", + " m2(2) = 1.100000\n", + " m(1) = 0.100000\n", + " m(2) = 1.050000\n", + " y1(0.1) = 1.010000\n", + " y2(0.1) = -0.895000\n", + "\n" + ] + } + ], + "source": [ + "from __future__ import division\n", + "\n", + "#System of differetial Equations\n", + "\n", + "def f1(x,y1,y2):\n", + " F1 = x + y1 + y2 \n", + " return F1\n", + "def f2(x,y1,y2):\n", + " F2 = 1 + y1 + y2 \n", + " return F2\n", + "\n", + "x0 = 0 #\n", + "y10 = 1 #\n", + "y20 = -1 #\n", + "h = 0.1 #\n", + "m1= [f1( x0 , y10 , y20 )]\n", + "m1.append(f2( x0 , y10 , y20 ))\n", + "m2=[f1( x0+h , y10 + h*m1[0] , y20 + h*m1[1] )]\n", + "m2.append(f2( x0+h , y10 + h*m1[0] , y20 + h*m1[1] ))\n", + "m= [(m1[0] + m2[0])/2 ]\n", + "m.append((m1[1] + m2[1])/2)\n", + "\n", + "y1_0_1 = y10 + h*m[0]\n", + "y2_0_1 = y20 + h*m[1]\n", + "\n", + "print 'm1(1) = %f\\n m1(2) = %f\\n m2(1) = %f\\n m2(2) = %f\\n m(1) = %f\\n m(2) = %f\\n y1(0.1) = %f\\n y2(0.1) = %f\\n'%(m1[0],m1[1],m2[0],m2[1],m[0],m[1],y1_0_1,y2_0_1) " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 13_14 Pg No. 457" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "m1(1) = 1.000000\n", + " m1(2) = -2.000000\n", + " m2(1) = 0.600000\n", + " m2(2) = 0.600000\n", + " m(1) = 0.800000\n", + " m(2) = -0.700000\n", + " y1(0.1) = 0.160000\n", + " y2(0.1) = 0.860000\n", + "\n" + ] + } + ], + "source": [ + "from __future__ import division\n", + "#Higher Order Differential Equations\n", + "\n", + "x0 = 0\n", + "y10 = 0\n", + "y20 = 1\n", + "h = 0.2\n", + "m1 = [y20] #\n", + "m1.append(6*x0 + 3*y10 - 2*y20)\n", + "m2= [y20 + h*m1[1]]\n", + "m2.append(6*(x0+h) + 3*(y10 + h*m1[0]) - 2*(y20 + h*m1[1]) )\n", + "m = [(m1[0] + m2[0])/2 ]\n", + "m.append((m1[1] + m2[1])/2)\n", + "\n", + "y1_0_2 = y10 + h*m[0]\n", + "y2_0_2 = y20 + h*m[1]\n", + "\n", + "print 'm1(1) = %f\\n m1(2) = %f\\n m2(1) = %f\\n m2(2) = %f\\n m(1) = %f\\n m(2) = %f\\n y1(0.1) = %f\\n y2(0.1) = %f\\n'%(m1[0],m1[1],m2[0],m2[1],m[0],m[1],y1_0_2,y2_0_2) " + ] + } + ], + "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 +} diff --git a/Numerical_Methods_by_E._Balaguruswamy/chapter14.ipynb b/Numerical_Methods_by_E._Balaguruswamy/chapter14.ipynb new file mode 100644 index 00000000..dbc3001f --- /dev/null +++ b/Numerical_Methods_by_E._Balaguruswamy/chapter14.ipynb @@ -0,0 +1,374 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Chapter 14 - Boundary value and eigen value problems" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 14_01 Pg No. 467" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "B1 = 7.75\n", + "Since B1 is less than B , let z(1) = y(1) = 4*(M2)\n", + "B2 = 9.75\n", + "Since B2 is larger than B ,let us have third estimate of z(1) = M3 \n", + "B3 = 9.0\n", + "The solution is [2, 4.375, 9.0]\n" + ] + } + ], + "source": [ + "#Shooting Method\n", + "\n", + "def heun(f,x0,y0,z0,h,xf):\n", + " x = [x0] #\n", + " global y\n", + " y = [y0] #\n", + " z = [z0] #\n", + " n = (xf - x0)/h\n", + " m1=[0,0];m2=[0,0];m=[0,0]\n", + " for i in range(0,int(n)):\n", + " m1[0] = z[i] \n", + " m1[1] = f(x[i],y[i])\n", + " m2[0] = z[i] + h*m1[1]\n", + " m2[1] = f(x[i]+h,y[i]+h*m1[0])\n", + " m[0] = (m1[0] + m2[0])/2 \n", + " m[1] = ( m1[1] + m2[1] )/2\n", + " x.append(x[(i)] + h)\n", + " y.append(y[(i)] + h*m[0])\n", + " z.append(z[(i)] + h*m[1])\n", + " \n", + " B = y[int(n)]\n", + " return B\n", + "\n", + "\n", + "def f(x,y):\n", + " F = 6*x\n", + " return F\n", + "\n", + "x0 = 1 #\n", + "y0 = 2 #\n", + "h = 0.5 #\n", + "z0 = 2\n", + "M1 = z0 \n", + "xf = 2\n", + "B = 9\n", + "B1 = heun(f,x0,y0,z0,h,xf)\n", + "print 'B1 = ',B1\n", + "\n", + "if B1 != B:\n", + " print 'Since B1 is less than B , let z(1) = y(1) = 4*(M2)'\n", + " z0 = 4\n", + " M2 = z0\n", + " B2 = heun(f,x0,y0,z0,h,xf)\n", + " print 'B2 = ',B2\n", + " if B2 != B:\n", + " print 'Since B2 is larger than B ,let us have third estimate of z(1) = M3 '\n", + " M3 = M2 - (B2 - B)*(M2 - M1)/(B2 - B1)\n", + " z0 = M3 \n", + " B3= heun(f,x0,y0,z0,h,xf)\n", + " print 'B3 = ',B3\n", + " print 'The solution is ',y" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 14_02 Pg No. 470" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + " \n", + " The solution is \n", + " y1 = y(0.25) = -0.100203 \n", + " y2 = y(0.5) = -0.137907 \n", + " y3 = y(0.75) = -0.109079 \n", + " \n" + ] + } + ], + "source": [ + "from numpy import array,exp,zeros,vstack,hstack,linalg\n", + "#Finite Difference Method\n", + "\n", + "def d2y(x):\n", + " D2Y = exp(x**2)\n", + " return D2Y\n", + "x_1 = 0#\n", + "y_0 = 0 #\n", + "y_1 = 0 #\n", + "h = 0.25\n", + "xf = 1\n", + "n = (xf-x_1)/h\n", + "A=zeros([int(n-1),int(n-1)])\n", + "B=zeros([int(n-1),1])\n", + "for i in range(0,int(n-1)):\n", + " A[i,:] = [1, -2, 1]\n", + " B[i,0] = exp((x_1 + i*h)**2)*h**2\n", + "\n", + "A[0,0] = 0 # #since we know y0 and y4\n", + "A[2,2] = 0 #\n", + "A[0,:] = hstack([ A[0,1:3], [0]]) #rearranging terms\n", + "A[2,0:3] = hstack([[ 0], A[2,0:2]]) \n", + "C=linalg.solve(A,B)\n", + "print ' \\n The solution is \\n y1 = y(0.25) = %f \\n y2 = y(0.5) = %f \\n y3 = y(0.75) = %f \\n '%(C[0],C[1],C[2]) " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 14_03 Pg No. 473" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\n", + " The roots are \n", + " lamda1 = 4.000000 \n", + " lamda2 = 6.000000 \n", + " \n", + "X1 = Matrix([[1.00000000000000], [1]])\n", + "X2 = Matrix([[2.00000000000000], [1]])\n" + ] + } + ], + "source": [ + "#Eigen Vectors\n", + "from numpy import eye\n", + "from sympy import symbols,Matrix,det,solve\n", + "A = [[8 ,-4],[ 2, 2 ] ]\n", + "A=Matrix(A)\n", + "lamd = symbols('lamd')\n", + "p = det(A - lamd*eye(2))\n", + "root = solve(p,lamd)\n", + "print '\\n The roots are \\n lamda1 = %f \\n lamda2 = %f \\n '%(root[0],root[1])\n", + "A1 = A - root[0]*eye(2)\n", + "X1 = Matrix([[-1*A1[0,1]/A1[0,0]],[1]])\n", + "print 'X1 = ',X1\n", + "A2 = A - root[1]*eye(2)\n", + "X2 = Matrix([[-1*A2[0,1]/A2[0,0]],[ 1]])\n", + "print 'X2 = ',X2" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 14_04 Pg No. 474" + ] + }, + { + "cell_type": "code", + "execution_count": 28, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\n", + "A2 = \n", + "[[-5. 0. 0.]\n", + " [ 1. -8. 9.]\n", + " [ 0. 4. -9.]]\n", + "\n", + "p2 = -11.000000\n", + "\n", + "\n", + "A3 = \n", + "[[ -6. 0. 0.]\n", + " [ 1. -6. 27.]\n", + " [ 0. 8. -6.]]\n", + "\n", + "p3 = -6.000000\n", + "\n" + ] + } + ], + "source": [ + "from numpy import array,shape,trace,poly\n", + "#Fadeev - Leverrier method\n", + "\n", + "A = [[ -1, 0, 0],[ 1, -2, 3],[ 0, 2, -3 ]]\n", + "A=array(A)\n", + "r,c= shape(A)[0],shape(A)[1]\n", + "A1 = A\n", + "p= [trace(A1)]\n", + "for i in range(1,r):\n", + " A1 = A*( A1 - p[(i-1)]*eye(3))\n", + " p.append(trace(A1)/(i+1))\n", + " print '\\nA%d = '%(i+1)\n", + " print A1\n", + " print '\\np%d = %f\\n'%((i+1),p[i])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 14_05 Pg No. 476" + ] + }, + { + "cell_type": "code", + "execution_count": 46, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\n", + " Eigen vector - 1 \n", + " for lamda1 = 0.000000 \n", + " X1 = \n", + "[-inf nan nan]\n", + "\n", + " Eigen vector - 2 \n", + " for lamda2 = 0.000000 \n", + " X2 = \n", + "[ 0. 0. 0.]\n", + "\n", + " Eigen vector - 3 \n", + " for lamda3 = 0.000000 \n", + " X3 = \n", + "[ 0. 0. -1.]\n" + ] + } + ], + "source": [ + "import numpy as np\n", + "np.seterr(divide='ignore', invalid='ignore')\n", + "from __future__ import division\n", + "from numpy.linalg import eig\n", + "from numpy import array,diag\n", + "#Eigen Vectors\n", + "\n", + "A = [[ -1, 0, 0],[1, -2, 3],[0, 2, -3]]\n", + "A=array(A)\n", + "evectors,evalues = eig(A)\n", + "evectors=diag(evectors)\n", + "for i in range(0,3):\n", + " print '\\n Eigen vector - %d \\n for lamda%d = %f \\n X%d = '%(i+1,i+1,evalues[0,0],i+1)\n", + " evectors[:,0] = evectors[:,0]/evectors[1,i]\n", + " print evectors[:,i]" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 14_06 Pg No. 478" + ] + }, + { + "cell_type": "code", + "execution_count": 48, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Iterations:\n", + " 0 1 2 3 4 5 6 7 \n", + "X = [[ 0. 1. 0.8 1. 0.97560976 1.\n", + " 0.99726027 1. ]\n", + " [ 1. 0.5 1. 0.92857143 1. 0.99180328\n", + " 1. 0.99908592]\n", + " [ 0. 0. 0. 0. 0. 0. 0.\n", + " 0. ]]\n", + "\n", + "Y = [[None 2.0 2.0 2.8 2.8571428571428577 2.975609756097561 2.9836065573770494\n", + " 2.9972602739726026]\n", + " [None 1.0 2.5 2.6 2.928571428571429 2.951219512195122 2.9918032786885247\n", + " 2.9945205479452053]\n", + " [None 0.0 0.0 0.0 0.0 0.0 0.0 0.0]]\n" + ] + } + ], + "source": [ + "from numpy import array, zeros,expand_dims, dot, hstack\n", + "A = array([[ 1, 2, 0],[2, 1 ,0],[0, 0, -1 ]])\n", + "X=zeros([3,8])\n", + "Y=zeros([3,7])\n", + "X[:,0] =[0,1,0]\n", + "\n", + "for i in range(1,8):\n", + " Y[:,i-1] = [xx for xx in A.dot(expand_dims(X[:,i-1], axis=1))]\n", + " \n", + " X[:,i] = Y[:,i-1]/max(Y[:,i-1])\n", + "\n", + "print 'Iterations:\\n 0 1 2 3 4 5 6 7 '\n", + "print 'X = ',X\n", + "Y=hstack([[[None],[None],[None]],Y])\n", + "print '\\nY = ',Y\n" + ] + } + ], + "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 +} diff --git a/Numerical_Methods_by_E._Balaguruswamy/chapter15.ipynb b/Numerical_Methods_by_E._Balaguruswamy/chapter15.ipynb new file mode 100644 index 00000000..0e991840 --- /dev/null +++ b/Numerical_Methods_by_E._Balaguruswamy/chapter15.ipynb @@ -0,0 +1,404 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Chapter 15 - Solution of partial differential equations" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 15_01 Pg No. 488" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\n", + " The solution of the system is \n", + " f1 = 75 \n", + " f2 = 50 \n", + " f3 = 50 \n", + " f4 = 24 \n", + " \n" + ] + } + ], + "source": [ + "from numpy import zeros,mat\n", + "from numpy.linalg import solve\n", + "#Elliptic Equations\n", + "\n", + "l = 15\n", + "h = 5\n", + "n = 1 + 15/5\n", + "f=zeros([4,4])\n", + "f[0,0:3]=100\n", + "f[0:3,0]=100\n", + "\n", + "#At point 1 : f2 + f3 - 4f1 + 100 + 100 = 0\n", + "#At point 2 : f1 + f4 - 4f2 + 100 + 0 = 0\n", + "#At point 3 : f1 + f4 - 4f3 + 100 + 0 = 0\n", + "#At point 4 : f2 + f3 - 4f4 + 0 + 0 = 0\n", + "#\n", + "#Final Equations are\n", + "# -4f1 + f2 + f3 + 0 = -200\n", + "# f1 - 4f2 + 0 + f4 = -100\n", + "# f1 + 0 - 4f3 + f4 = -100\n", + "# 0 + f2 + f3 - 4f4 = 0\n", + "A = mat([[ -4, 1 ,1 ,0],[1, -4, 0 ,1],[1, 0 ,-4 ,1],[0, 1, 1, -4 ]])\n", + "B = mat([[-200],[-100],[-100],[0]])\n", + "C = solve(A,B)\n", + "print '\\n The solution of the system is \\n f1 = %d \\n f2 = %d \\n f3 = %d \\n f4 = %d \\n '%(C[0],C[1],C[2],C[3])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 15_02 Pg No. 489" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "First row of below matrix represents iteration number:\n", + "[[ 0. 1. 2. 3. 4.]\n", + " [ 75. 75. 75. 75. 75.]\n", + " [ 0. 0. 0. 0. 0.]\n", + " [ 0. 0. 0. 0. 0.]\n", + " [ 0. 50. 50. 50. 50.]]\n" + ] + } + ], + "source": [ + "from numpy import zeros,mat\n", + "from numpy.linalg import solve\n", + "\n", + "#Liebmann's Iterative method\n", + "\n", + "f=zeros([4,4])\n", + "f[0,0:3]=100\n", + "f[0:3,0]=100\n", + "A=zeros([5,5])\n", + "A[0,0:5] = range(0,5)\n", + "for n in range(1,6):\n", + " for i in range(2,4):\n", + " for j in range(2,4):\n", + " if n == 1 and i == 2 and j == 2 :\n", + " f[i-1,j-1] = ( f[i,j] + f[i-2,j-2] + f[i-2,j] + f[i,j-2] )/4\n", + " else:\n", + " f[i,j] = ( f[i,j-1] + f[i-2,j-1] + f[i-1,j]+ f[i-1,j-2] )/4\n", + " A[1:5,n-1] = mat([f[1,1],f[1,2],f[2,1],f[2,2] ])\n", + "\n", + "\n", + "print 'First row of below matrix represents iteration number:\\n',A" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 15_03 Pg No. 490" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "The solution is \n", + " f1 = -5.500000 \n", + " f2 = -10.750000 \n", + " f3 = -3.250000 \n", + " f4 = -5.500000 \n", + " \n" + ] + } + ], + "source": [ + "from numpy import zeros,mat\n", + "from numpy.linalg import solve\n", + "\n", + "#Poisson's Equation\n", + "\n", + "#D2f = 2*x**2 * y**2\n", + "# f = 0\n", + "# h = 1 \n", + "#Point 1 : 0 + 0 + f2 + f3 - 4f1 = 2(1)**2 * 2**2\n", + "# f2 + f3 - 4f1 = 8\n", + "#Point 2 : 0 + 0 + f1 + f4 -4f2 = 2*(2)**2*2**2\n", + "# f1 - 4f2 = f4 = 32\n", + "#Point 3 : 0 + 0 + f1 + f4 - 4f4 = 2*(1**2)*1**2\n", + "# f1 -4f3 + f4 = 2\n", + "#Point 4 : 0 + 0 + f2 + f3 - 4f4 = 2* 2**2 * 1**2\n", + "# f2 + f3 - 4f4 = 8\n", + "#Rearranging the equations\n", + "# -4f1 + f2 + f3 = 8\n", + "# f1 - 4f2 + f4 = 32\n", + "# f1 - 4f3 + f4 = 2\n", + "# f2 + f3 - 4f4 = 8\n", + "A = mat([[ -4, 1, 1, 0],[1, -4 ,0 ,1],[1, 0, -4, 1],[0, 1, 1, -4]])\n", + "B = mat([[ 8],[32],[2],[8 ]])\n", + "C = solve(A,B)\n", + "print 'The solution is \\n f1 = %f \\n f2 = %f \\n f3 = %f \\n f4 = %f \\n '%( C[0],C[1],C[2],C[3])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 15_04 Pg No. 491" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\n", + " Iteration 1\n", + " f1 = -2.000000, f2 = -9.000000, f3 = -2.000000, f4 = -2.000000\n", + "\n", + "\n", + " Iteration 2\n", + " f1 = -5.000000, f2 = -11.000000, f3 = -3.000000, f4 = -5.000000\n", + "\n", + "\n", + " Iteration 3\n", + " f1 = -6.000000, f2 = -11.000000, f3 = -4.000000, f4 = -6.000000\n", + "\n", + "\n", + " Iteration 4\n", + " f1 = -6.000000, f2 = -11.000000, f3 = -4.000000, f4 = -6.000000\n", + "\n" + ] + } + ], + "source": [ + "#Gauss-Seidel Iteration\n", + "\n", + "f2 = 0\n", + "f3 = 0\n", + "for i in range(1,5):\n", + " f1 = (f2 + f3 - 8)/4\n", + " f4 = f1 \n", + " f2 = (f1 + f4 -32)/4\n", + " f3 = (f1 + f4 - 2)/4\n", + " print '\\n Iteration %d\\n f1 = %f, f2 = %f, f3 = %f, f4 = %f\\n'%(i,f1,f2,f3,f4)\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 15_05 Pg No. 494" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "The final results are \n", + "[[ 0. 150. 100. 50. 0. ]\n", + " [ 0. 50. 100. 50. 0. ]\n", + " [ 0. 50. 50. 50. 0. ]\n", + " [ 0. 25. 50. 25. 0. ]\n", + " [ 0. 25. 25. 25. 0. ]\n", + " [ 0. 12.5 25. 12.5 0. ]\n", + " [ 0. 12.5 12.5 12.5 0. ]]\n" + ] + } + ], + "source": [ + "from numpy import zeros,mat\n", + "\n", + "#Initial Value Problems\n", + "\n", + "h = 1 #\n", + "k = 2 #\n", + "tau = h**2/(2*k)\n", + "f=zeros([7,5])\n", + "for i in range(1,5):\n", + " f[0,i] = 50*( 4 - (i) )\n", + "\n", + "f[0:7,0] = 0\n", + "f[0:7,4] = 0 \n", + "for j in range(0,6):\n", + " for i in range(1,4):\n", + " f[j+1,i] = ( f[j,i-1] + f[j,i+1] )/2 \n", + "print 'The final results are \\n',f" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 15_06 Pg No. 497" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "The final solution using crank nicholson implicit method is :\n", + "[[ 0. 150. 100. 50. 0. ]\n", + " [ 0. 42.85714286 71.42857143 42.85714286 0. ]\n", + " [ 0. 26.53061224 34.69387755 26.53061224 0. ]\n", + " [ 0. 13.70262391 20.11661808 13.70262391 0. ]\n", + " [ 0. 7.70512287 10.70387339 7.70512287 0. ]]\n" + ] + } + ], + "source": [ + "from numpy import zeros,mat,linalg\n", + "#Crank-Nicholson Implicit Method\n", + "\n", + "h = 1 #\n", + "k = 2 #\n", + "tau = h**2/(2*k)\n", + "f=zeros([5,5])\n", + "for i in range(1,4):\n", + " f[0,i] = 50*( 4 - (i) )\n", + "f[0:5,0] = 0 \n", + "f[0:5,4] = 0 \n", + "A = mat([[4, -1, 0 ],[-1 , 4 , -1 ],[0 , -1 , 4]])\n", + "B=zeros([3,1])\n", + "for j in range(0,4):\n", + " for i in range(1,4):\n", + " B[i-1,0] = f[j,i-1] + f[j,i+1]\n", + " \n", + " C = linalg.solve(A,B)\n", + " f[j+1,1] = C[0]\n", + " f[j+1,2] = C[1]\n", + " f[j+1,3] = C[2]\n", + "\n", + "print 'The final solution using crank nicholson implicit method is :\\n',f" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 15_07 Pg No. 500" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "The values estimated are :\n", + "[[ 0. 4. 6. 6. 4. 0.]\n", + " [ 0. 3. 5. 5. 3. 0.]\n", + " [ 0. 1. 2. 2. 1. 0.]\n", + " [ 0. -1. -2. -2. -1. 0.]\n", + " [ 0. -3. -5. -5. -3. 0.]\n", + " [ 0. -4. -6. -6. -4. 0.]]\n" + ] + } + ], + "source": [ + "from __future__ import division\n", + "from numpy import sqrt\n", + "#Hyperbolic Equations\n", + "\n", + "h = 1\n", + "Tbyp = 4\n", + "tau = sqrt(h**2 /4)\n", + "r = int(1+(2.5 - 0)/tau)\n", + "c = int(1+(5 - 0)/h)\n", + "\n", + "f=zeros([6,6])\n", + "for i in range(1,int(c)-1):\n", + " f[0,i] = (i)*(5 - (i) )\n", + "\n", + "f[0:r-1,0] = 0\n", + "f[0:r-1,c-1] = 0\n", + "g=[]\n", + "for i in range(1,c-1):\n", + " g.append(0)\n", + " f[1,i] = (f[0,i+1] + f[0,i-1])/2 + tau*g[i-1] \n", + "\n", + "for j in range(1,r-1):\n", + " for i in range(1,c-1):\n", + " f[j+1,i] = -f[j-1,i] + f[j,i+1] + f[j,i-1]\n", + " \n", + "\n", + "print 'The values estimated are :\\n',f" + ] + } + ], + "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 +} diff --git a/Numerical_Methods_by_E._Balaguruswamy/chapter3.ipynb b/Numerical_Methods_by_E._Balaguruswamy/chapter3.ipynb new file mode 100644 index 00000000..d4901132 --- /dev/null +++ b/Numerical_Methods_by_E._Balaguruswamy/chapter3.ipynb @@ -0,0 +1,1170 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Chapter 3 - Computer codes and arithmetic" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 3_01 Pg No. 45" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Decimal value = 13.8125\n" + ] + } + ], + "source": [ + "#Binary to decimal\n", + "\n", + "b = '1101.1101'\n", + "def parse_bin(s):\n", + " t = s.split('.')\n", + " return int(t[0], 2) + int(t[1], 2) / 2.**len(t[1])\n", + "\n", + "d = parse_bin(b) # #Integral and fractional parts\n", + "print 'Decimal value = ',d" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 3_02 Pg No. 46" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Decimal value = 4783\n" + ] + } + ], + "source": [ + "#hexadecimal to decimal\n", + "\n", + "h = '12AF' \n", + "d=int(h,16)\n", + "print 'Decimal value = ',d" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 3_03 Pg No. 47" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Binary equivalent = 101011.011\n" + ] + } + ], + "source": [ + "#Decimal to Binary\n", + "\n", + "d = 43.375 #\n", + "\n", + "def parse_float(num):\n", + " exponent=0\n", + " shifted_num=num\n", + " while shifted_num != int(shifted_num): \n", + " shifted_num*=2\n", + " exponent+=1\n", + " if exponent==0:\n", + " return '{0:0b}'.format(int(shifted_num))\n", + " binary='{0:0{1}b}'.format(int(shifted_num),exponent+1)\n", + " integer_part=binary[:-exponent]\n", + " fractional_part=binary[-exponent:].rstrip('0')\n", + " return '{0}.{1}'.format(integer_part,fractional_part)\n", + "\n", + "b = parse_float(d) \n", + "print 'Binary equivalent = ',b" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 3_04 Pg No. 48" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Octal number = 243\n" + ] + } + ], + "source": [ + "#Decimal to Octal\n", + "\n", + "d = 163 #\n", + "def base10toN(num, base):\n", + " \"\"\"Change ``num'' to given base\n", + " Upto base 36 is supported.\"\"\"\n", + "\n", + " converted_string, modstring = \"\", \"\"\n", + " currentnum = num\n", + " if not 1 < base < 37:\n", + " raise ValueError(\"base must be between 2 and 36\")\n", + " if not num:\n", + " return '0'\n", + " while currentnum:\n", + " mod = currentnum % base\n", + " currentnum = currentnum // base\n", + " converted_string = chr(48 + mod + 7*(mod > 10)) + converted_string\n", + " return converted_string\n", + "\n", + "\n", + "\n", + "Oct = base10toN(d,8)\n", + "print 'Octal number = ',Oct" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 3_05 Pg No. 48" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Binary equivalent = 0.10100110011001100110011001100110011001100110011001101\n" + ] + } + ], + "source": [ + "#Decimal to binary\n", + "\n", + "d = 0.65\n", + "\n", + "def parse_float(num):\n", + " exponent=0\n", + " shifted_num=num\n", + " while shifted_num != int(shifted_num): \n", + " shifted_num*=2\n", + " exponent+=1\n", + " if exponent==0:\n", + " return '{0:0b}'.format(int(shifted_num))\n", + " binary='{0:0{1}b}'.format(int(shifted_num),exponent+1)\n", + " integer_part=binary[:-exponent]\n", + " fractional_part=binary[-exponent:].rstrip('0')\n", + " return '{0}.{1}'.format(integer_part,fractional_part)\n", + "\n", + "\n", + "b=parse_float(d) # binary equibvalent\n", + "print 'Binary equivalent = ',b" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 3_06 Pg No. 49 " + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Hexadecimal equivalent of octal number 243 is : a3\n" + ] + } + ], + "source": [ + "#Octal to Hexadecimal \n", + "\n", + "Oct = '243' \n", + "dec=int(Oct,8)\n", + "h = hex(dec)[2:]\n", + "\n", + "print 'Hexadecimal equivalent of octal number 243 is :',h" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 3_07 Pg No. 49" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Octal number equivalent of Hexadecimal number 39.B8 is: 57.71875\n" + ] + } + ], + "source": [ + "#Hexadecimal to Octal \n", + "\n", + "\n", + "h = '39.B8' #\n", + " \n", + " \n", + "def float_to_binary(num):\n", + " exponent=0\n", + " shifted_num=num\n", + " while shifted_num != int(shifted_num): \n", + " shifted_num*=2\n", + " exponent+=1\n", + " if exponent==0:\n", + " return '{0:0b}'.format(int(shifted_num))\n", + " binary='{0:0{1}b}'.format(int(shifted_num),exponent+1)\n", + " integer_part=binary[:-exponent]\n", + " fractional_part=binary[-exponent:].rstrip('0')\n", + " return '{0}.{1}'.format(integer_part,fractional_part)\n", + "\n", + "def floathex_to_binary(floathex):\n", + " num = float.fromhex(floathex)\n", + " return float_to_binary(num)\n", + "\n", + "b= floathex_to_binary(h)\n", + "def parse_bin(s):\n", + " t = s.split('.')\n", + " return int(t[0], 2) + int(t[1], 2) / 2.**len(t[1])\n", + "\n", + "d = parse_bin(b) # #Integral and fractional parts\n", + "print 'Octal number equivalent of Hexadecimal number 39.B8 is:',d\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 3_08 Pg No. 50" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Binary equivalent of -13 is : 10011\n" + ] + } + ], + "source": [ + "## -ve Integer to binary\n", + "\n", + "negint = -13\n", + "def to_twoscomplement(bits, value):\n", + " if value < 0:\n", + " value = ( 1<<bits ) + value\n", + " formatstring = '{:0%ib}' % bits\n", + " return formatstring.format(value)\n", + "\n", + "compl_2 = to_twoscomplement(5, negint)\n", + "\n", + "print 'Binary equivalent of -13 is :',compl_2" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 3_09 Pg No. 51" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Binary equivalent of -32768 in a 16 bit word is: 1000000000000000\n" + ] + } + ], + "source": [ + "#Binary representation\n", + "\n", + "n = -32768\n", + "\n", + "def to_twoscomplement(bits, value):\n", + " if value < 0:\n", + " value = ( 1<<bits ) + value\n", + " formatstring = '{:0%ib}' % bits\n", + " return formatstring.format(value)\n", + "\n", + "binfinal = to_twoscomplement(16, n)\n", + "\n", + "\n", + "print 'Binary equivalent of -32768 in a 16 bit word is:',binfinal " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 3_10 Pg No. 52" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\n", + " 0.00596 is expressed as 0.596000*10**-2 \n", + "\n", + "\n", + " 65.7452 is expressed as 0.657452*10**2 \n", + "\n", + "\n", + " -486.8 is expressed as -0.486800*10**3 \n", + "\n" + ] + } + ], + "source": [ + "#Floating Point Notation\n", + "\n", + "def float_notation(n):\n", + " m = n #\n", + " for i in range(1,17):\n", + " if abs(m) >= 1:\n", + " m = n/10**i\n", + " e = i\n", + " elif abs(m) < 0.1:\n", + " m = n*10**i\n", + " e = -i\n", + " else:\n", + " if i == 1:\n", + " e = 0\n", + " \n", + " break \n", + " return [m,e]\n", + " \n", + "\n", + "[m,e] = float_notation(0.00596)\n", + "print '\\n 0.00596 is expressed as %f*10**%d \\n'%(m,e)\n", + "[m,e] = float_notation(65.7452)\n", + "print '\\n 65.7452 is expressed as %f*10**%d \\n'%(m,e)\n", + "[m,e] = float_notation(-486.8)\n", + "print '\\n -486.8 is expressed as %f*10**%d \\n'%(m,e)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 3_11 Pg No. 53" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "25+12 = 37\n", + "25-12 = 13\n", + "12-25 = -13\n", + "25*12 = 300\n", + "25/12 = 2\n", + "12/25 = 0\n" + ] + } + ], + "source": [ + "#Interger Arithmetic\n", + "\n", + "print '25+12 =',int(25 + 12)\n", + "print '25-12 =',int(25 - 12) \n", + "print '12-25 =',int(12 - 25)\n", + "print '25*12 =',int(25*12)\n", + "print '25/12 =',int(25/12)\n", + "print '12/25 =',int(12/25)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 3_12 Pg No. 53" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "(a+b)/c = 4 \n", + " a/c + b/c = 3\n", + "The results are not identical.This is because the remainder of an integer division is always truncated\n" + ] + } + ], + "source": [ + "#Integer Arithmetic\n", + "\n", + "a = 5 #\n", + "b = 7 #\n", + "c = 3 #\n", + "Lhs = int((a + b)/c)\n", + "Rhs = int(a/c) + int(b/c)\n", + "print '(a+b)/c = ',Lhs,'\\n a/c + b/c = ',Rhs\n", + "if Lhs != Rhs:\n", + " print 'The results are not identical.This is because the remainder of an integer division is always truncated'\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 3_13 Pg No. 54" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Ez= 5 \n", + " fy= 0.000964572 \n", + " fz= 0.587315572\n", + "\n", + " z = 0.587316 E5 \n", + "\n" + ] + } + ], + "source": [ + "#Floating Point Arithmetic\n", + "\n", + "fx = 0.586351 #\n", + "Ex = 5 #\n", + "fy = 0.964572 #\n", + "Ey = 2 #\n", + "v=[Ex,Ey]\n", + "Ez= max(v)\n", + "n=v.index(Ez)+1\n", + "if n == 1:\n", + " fy = fy*10**(Ey-Ex)\n", + " fz = fx + fy\n", + " if fz > 1:\n", + " fz = fz*10**(-1) \n", + " Ez = Ez + 1\n", + " \n", + " print 'Ez=',Ez,'\\n fy=',fy,'\\n fz=',fz\n", + "else:\n", + " fx = fx*10**(Ex - Ey)\n", + " fz = fx + fy\n", + " if fz > 1:\n", + " fz = fz*10**(-1)\n", + " Ez = Ez + 1\n", + " \n", + " print 'Ez=',Ez,'\\n fy=',fy,'\\n fz=',fz\n", + "\n", + "print '\\n z = %f E%d \\n'%(fz,Ez)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 3_14 Pg No. 54" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Ez= 5 \n", + " fy= 0.635742 \n", + " fz= 0.1371558\n", + "\n", + " z = 0.137156 E5 \n", + "\n" + ] + } + ], + "source": [ + "#Floating Point Arithmetic\n", + "\n", + "fx = 0.735816 #\n", + "Ex = 4 #\n", + "fy = 0.635742 #\n", + "Ey = 4 #\n", + "v=[Ex,Ey]\n", + "Ez= max(v)\n", + "n=v.index(Ez)+1\n", + "if n == 1:\n", + " fy = fy*10**(Ey-Ex)\n", + " fz = fx + fy\n", + " if fz > 1:\n", + " fz = fz*10**(-1)\n", + " Ez = Ez + 1\n", + " \n", + " print 'Ez=',Ez,'\\n fy=',fy,'\\n fz=',fz\n", + "else:\n", + " fx = fx*10**(Ex - Ey)\n", + " fz = fx + fy\n", + " if fz > 1:\n", + " fz = fz*10**(-1) \n", + " Ez = Ez + 1\n", + " \n", + " print 'Ez=',Ez,'\\n fy=',fy,'\\n fz=',fz\n", + "\n", + "print '\\n z = %f E%d \\n'%(fz,Ez)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 3_15 Pg No. 54" + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Ez = -3 fz = 0.005082\n", + "\n", + " z = 0.005082 E-3 \n", + "\n", + "\n", + " z = 0.005082 E-3 (normalised) \n", + "\n" + ] + } + ], + "source": [ + "#Floating Point Arithmetic\n", + "\n", + "fx = 0.999658 #\n", + "Ex = -3 #\n", + "fy = 0.994576 #\n", + "Ey = -3 #\n", + "Ez = max(Ex,Ey)\n", + "fy = fy*10**(Ey-Ex)\n", + "fz = fx - fy\n", + "print 'Ez =',Ez,'fz =',fz\n", + "print '\\n z = %f E%d \\n'%(fz,Ez)\n", + "if fz < 0.1 :\n", + " fz = fz*10**6 #Since we are using 6 significant digits\n", + " n = len(str(fz))\n", + " fz = fz/10**n\n", + " Ez = Ez + n - 6\n", + " print '\\n z = %f E%d (normalised) \\n'%(fz,Ez)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 3_16 Pg No. 55" + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\n", + " fz = 0.080000 \n", + " Ez = 2 \n", + " z = 0.080000 E2 \n", + "\n", + "\n", + " z = 0.800000 E1 (normalised) \n", + "\n" + ] + } + ], + "source": [ + "#Floating Point Arithmetic\n", + "\n", + "fx = 0.200000 #\n", + "Ex = 4 #\n", + "fy = 0.400000 #\n", + "Ey = -2 #\n", + "fz = fx*fy\n", + "Ez = Ex + Ey \n", + "print '\\n fz = %f \\n Ez = %d \\n z = %f E%d \\n'%(fz,Ez,fz,Ez)\n", + "if fz < 0.1:\n", + " fz = fz*10\n", + " Ez = Ez - 1\n", + " print '\\n z = %f E%d (normalised) \\n'%(fz,Ez)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 3_17 Pg No. 55" + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\n", + " fz = 4.382715 \n", + " Ez = -2 \n", + " z = 4.382715 E-2 \n", + "\n", + "\n", + " z = 0.438271 E-1 (normalised) \n", + "\n" + ] + } + ], + "source": [ + "#Floating Point Arithmetic\n", + "\n", + "fx = 0.876543 #\n", + "Ex = -5 #\n", + "fy = 0.200000 #\n", + "Ey = -3 #\n", + "fz = fx/fy\n", + "Ez = Ex - Ey \n", + "print '\\n fz = %f \\n Ez = %d \\n z = %f E%d \\n'%(fz,Ez,fz,Ez)\n", + "\n", + "if fz > 1:\n", + " fz = fz/10\n", + " Ez = Ez + 1\n", + " print '\\n z = %f E%d (normalised) \\n'%(fz,Ez)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 3_18 Pg No. 56" + ] + }, + { + "cell_type": "code", + "execution_count": 18, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Ez = 2 \n", + " fx = 50000000.0\n", + "\n", + " fz = 5000000.010000 \n", + " z = 5000000.010000 E2 \n", + "\n" + ] + } + ], + "source": [ + "#Floating Point Arithmetic\n", + "\n", + "fx = 0.500000 #\n", + "Ex = 1 #\n", + "fy = 0.100000 #\n", + "Ey = -7 #\n", + "Ez= max([Ex,Ey])\n", + "n=[Ex,Ey].index(Ez)\n", + "if n == 1:\n", + " fy = fy*10**(Ey-Ex)\n", + " fz = fx + fy\n", + " if fz > 1:\n", + " fz = fz*10**(-1) \n", + " Ez = Ez + 1\n", + " \n", + " print 'Ez =',Ez,'\\n fy =',fy\n", + "else:\n", + " fx = fx*10**(Ex - Ey)\n", + " fz = fx + fy\n", + " if fz > 1:\n", + " fz = fz*10**(-1)\n", + " Ez = Ez + 1\n", + " \n", + " print 'Ez =',Ez,'\\n fx =',fx\n", + "\n", + "print '\\n fz = %f \\n z = %f E%d \\n'%(fz,fz,Ez)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 3_19 Pg No. 56" + ] + }, + { + "cell_type": "code", + "execution_count": 19, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\n", + " fz = 0.175000 \n", + " Ez = 110 \n", + " z = 0.175000 E110 \n", + "\n" + ] + } + ], + "source": [ + "#Floating Point Arithmetic\n", + "\n", + "fx = 0.350000 #\n", + "Ex = 40 #\n", + "fy = 0.500000 #\n", + "Ey = 70 #\n", + "fz = fx*fy\n", + "Ez = Ex + Ey \n", + "print '\\n fz = %f \\n Ez = %d \\n z = %f E%d \\n'%(fz,Ez,fz,Ez)\n", + "if fz < 0.1:\n", + " fz = fz*10\n", + " Ez = Ez - 1\n", + " print '\\n z = %f E%d (normalised) \\n'%(fz,Ez)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 3_20 Pg No. 56" + ] + }, + { + "cell_type": "code", + "execution_count": 20, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\n", + " fz = 4.375000 \n", + " Ez = -113 \n", + " z = 4.375000 E-113 \n", + "\n", + "\n", + " z = 0.437500 E-112 (normalised) \n", + "\n" + ] + } + ], + "source": [ + "#Floating Point Arithmetic\n", + "\n", + "fx = 0.875000 #\n", + "Ex = -18 #\n", + "fy = 0.200000 #\n", + "Ey = 95 #\n", + "fz = fx/fy\n", + "Ez = Ex - Ey \n", + "print '\\n fz = %f \\n Ez = %d \\n z = %f E%d \\n'%(fz,Ez,fz,Ez)\n", + "\n", + "if fz > 1:\n", + " fz = fz/10\n", + " Ez = Ez + 1\n", + " print '\\n z = %f E%d (normalised) \\n'%(fz,Ez)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 3_21 Pg No. 57" + ] + }, + { + "cell_type": "code", + "execution_count": 21, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Ez = 0 \n", + " fz = 2e-06\n", + "\n", + " z = 0.000002 E0 \n", + "\n", + "\n", + " z = 0.002000 E-3 (normalised) \n", + "\n" + ] + } + ], + "source": [ + "#Floating Point Arithmetic\n", + "\n", + "fx = 0.500000 #\n", + "Ex = 0 #\n", + "fy = 0.499998 #\n", + "Ey = 0 #\n", + "Ez = 0 #\n", + "fz = fx - fy\n", + "print 'Ez =',Ez,'\\n fz =',fz\n", + "print '\\n z = %f E%d \\n'%(fz,Ez)\n", + "if fz < 0.1 :\n", + " fz = fz*10**6\n", + " n = len(str(fz))\n", + " fz = fz/10**n\n", + " Ez = Ez + n - 6\n", + " print '\\n z = %f E%d (normalised) \\n'%(fz,Ez)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 3_22 Pg No. 57" + ] + }, + { + "cell_type": "code", + "execution_count": 22, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + " fxy = 2.480183\n", + " Exy = 1 \n", + " fxy_z = 24.553832\n", + " Exy_z = 1 \n", + " fyz = -1.000000 \n", + " Eyz = -2 \n", + " fx_yz = -1.000000 \n", + " Ex_yz = 0 \n", + "\n", + " (x+y) + z != x + (y+z)\n" + ] + } + ], + "source": [ + "#Laws of Arithmetic\n", + "\n", + "def add_sub(fx,Ex,fy,Ey): #addition and subtraction fuction\n", + " if fx*fy >= 0:\n", + " #Addition\n", + " Ez = max([Ex,Ey])\n", + " n=[Ex,Ey].index(Ez)\n", + " if n == 1:\n", + " fy = fy*10**(Ey-Ex)\n", + " fz = fx + fy\n", + " if fz > 1:\n", + " fz = fz*10**(-1)\n", + " Ez = Ez + 1\n", + " \n", + " else:\n", + " fx = fx*10**(Ex - Ey)\n", + " fz = fx + fy\n", + " if fz > 1:\n", + " fz = fz*10**(-1) \n", + " Ez = Ez + 1\n", + " \n", + " \n", + " \n", + " else:\n", + " #Subtraction\n", + " Ez = max([Ex,Ey])\n", + " n=[Ex,Ey].index(Ez)\n", + " if n == 1:\n", + " fy = fy*10**(Ey-Ex)\n", + " fz = fx + fy\n", + " if abs(fz) < 0.1:\n", + " fz = fz*10**6\n", + " fz = floor(fz)\n", + " nfz = len(str(abs(fz)))\n", + " fz = fz/10**nfz\n", + " Ez = nfz - 6 \n", + " \n", + " else:\n", + " fx = fx*10**(Ex - Ey)\n", + " fz = fx + fy\n", + " if fz < 0.1:\n", + " fz = fz*10**6\n", + " fz = int(fz)\n", + " nfz = len(str(abs(fz)))\n", + " fz = fz/10**nfz\n", + " Ez = nfz - 6\n", + " \n", + " \n", + " \n", + " return [fz,Ez]\n", + "\n", + "fx = 0.456732\n", + "Ex = -2\n", + "fy = 0.243451\n", + "Ey = 0\n", + "fz = -0.24800\n", + "Ez = 0\n", + "\n", + "[fxy,Exy] = add_sub(fx,Ex,fy,Ey)\n", + "[fxy_z,Exy_z] = add_sub(fxy,Exy,fz,Ez)\n", + "[fyz,Eyz] = add_sub(fy,Ey,fz,Ez)\n", + "[fx_yz,Ex_yz] = add_sub(fx,Ex,fyz,Eyz)\n", + "print ' fxy = %f\\n Exy = %d \\n fxy_z = %f\\n Exy_z = %d \\n fyz = %f \\n Eyz = %d \\n fx_yz = %f \\n Ex_yz = %d \\n'%(fxy,Exy,fxy_z,Exy_z,fyz,Eyz,fx_yz,Ex_yz)\n", + "\n", + "if fxy_z != fx_yz | Exy_z != Ex_yz:\n", + " print ' (x+y) + z != x + (y+z)' " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 3_23 Pg No. 58" + ] + }, + { + "cell_type": "code", + "execution_count": 23, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "In book they have considered the maximum exponent can be only 99, since 110 is greater than 99 the result is erroneous\n", + "but in scilab the this value is much larger than 110 so we get a correct result \n", + "xy_z = 6e+78\n", + "x_yz = 6e+78\n" + ] + } + ], + "source": [ + "#Associative law\n", + "x = 0.400000*10**40\n", + "y = 0.500000*10**70\n", + "z = 0.300000*10**(-30)\n", + "print 'In book they have considered the maximum exponent can be only 99, since 110 is greater than 99 the result is erroneous'\n", + "print 'but in scilab the this value is much larger than 110 so we get a correct result '\n", + "print 'xy_z = ',(x*y)*z\n", + "print 'x_yz = ',x*(y*z)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 3_24 Pg No. 58" + ] + }, + { + "cell_type": "code", + "execution_count": 24, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "x_yz = 4e-06\n", + "xy_xz = 0.0\n" + ] + } + ], + "source": [ + "from math import floor\n", + "#Distributive law\n", + "x = 0.400000*10**1 #\n", + "fx = 0.400000\n", + "Ex = 1\n", + "y = 0.200001*10**0 #\n", + "z = 0.200000*10**0 #\n", + "x_yz = x*(y-z)\n", + "x_yz = x_yz*10**6\n", + "x_yz = floor(x_yz) #considering only six significant digits\n", + "n = len(str(x_yz))\n", + "fx_yz = x_yz/10**n\n", + "Ex_yz = n - 6\n", + "x_yz = fx_yz *10**Ex_yz\n", + "print 'x_yz = ',x_yz\n", + "\n", + "fxy = fx*y\n", + "fxy = fxy*10**6\n", + "fxy = floor(fxy) #considering only six significant digits\n", + "n = len(str(fxy))\n", + "fxy = fxy/10**n\n", + "Exy = n - 6\n", + "xy = fxy * 10**Exy\n", + "\n", + "fxz = fx*z\n", + "fxz = fxz*10**6\n", + "fxz = floor(fxz) #considering only six significant digits\n", + "n = len(str(fxz))\n", + "fxz = fxz/10**n\n", + "Exz = n - 6\n", + "xz = fxz * 10**Exz\n", + "\n", + "xy_xz = xy - xz\n", + "print 'xy_xz = ',xy_xz" + ] + } + ], + "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 +} diff --git a/Numerical_Methods_by_E._Balaguruswamy/chapter4.ipynb b/Numerical_Methods_by_E._Balaguruswamy/chapter4.ipynb new file mode 100644 index 00000000..832b55d6 --- /dev/null +++ b/Numerical_Methods_by_E._Balaguruswamy/chapter4.ipynb @@ -0,0 +1,891 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Chapter 4 - Approximations and errors in computing" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 4_01 Pg No. 63" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\n", + " 4.3201 has a precision of 10**-4\n", + "\n", + "\n", + " 4.32 has a precision of 10**-2\n", + "\n", + "\n", + " 4.320106 has a precision of 10**-6\n", + "\n", + "\n", + " The number with highest precision is 4.320106\n", + "\n" + ] + } + ], + "source": [ + "#Greatest precision\n", + "\n", + "a = '4.3201'\n", + "b = '4.32'\n", + "c = '4.320106'\n", + "na = len(a)-(a.index('.')+1)\n", + "print '\\n %s has a precision of 10**-%d\\n'%(a,na)\n", + "nb = len(b)-(b.index('.')+1)\n", + "print '\\n %s has a precision of 10**-%d\\n'%(b,nb)\n", + "nc = len(c)-c.index('.')-1\n", + "print '\\n %s has a precision of 10**-%d\\n'%(c,nc)\n", + "e = max(na,nb,nc)\n", + "if e ==na:\n", + " print '\\n The number with highest precision is %s\\n'%(a)\n", + "elif e == nb:\n", + " print '\\n The number with highest precision is %s\\n'%(b)\n", + "else:\n", + " print '\\n The number with highest precision is %s\\n'%(c)\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 4_02 Pg No. 63" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "95.763 has 5 significant digits\n", + "\n", + "0.008472 has 7 significant digits.The leading or higher order zeros are only place holders\n", + "\n", + "0.0456000 has 8 significant digits\n", + "\n", + "36.0 has 3 significant digits\n", + "\n", + "3600.00 has 6 significant digits\n", + "\n" + ] + } + ], + "source": [ + "#Accuracy of numbers\n", + "\n", + "def sd(x):\n", + " nd = x.index('.') #position of point\n", + " num = [ord(c) for c in x] #str2code(x)\n", + " if (nd)==None and num(length(x)) == 0:\n", + " print 'Accuracy is not specified\\n'\n", + " n = 0 #\n", + " else:\n", + " if num[0]>= 1 and (nd==None):\n", + " n = len(x)\n", + " elif num[1] >= 1 and not((nd==None)):\n", + " n = len(x) - 1\n", + " else:\n", + " for i in range(0,len(x)):\n", + " if num[i] >= 1 and num[i] <= 9:\n", + " break;\n", + " \n", + " \n", + " n = len(x)- i + 1\n", + " return n\n", + " \n", + " \n", + "\n", + "a = '95.763'\n", + "na = sd(a)\n", + "print '%s has %d significant digits\\n'%(a,na)\n", + "b = '0.008472'\n", + "nb = sd(b)\n", + "print '%s has %d significant digits.The leading or higher order zeros are only place holders\\n'%(b,nb)\n", + "c = '0.0456000'\n", + "nc = sd(c)\n", + "print '%s has %d significant digits\\n'%(c,nc)\n", + "d = '36.0'\n", + "nd = sd(d)\n", + "print '%s has %d significant digits\\n'%(d,nd)\n", + "e = '3600.0'\n", + "sd(e)\n", + "f = '3600.00'\n", + "nf = sd(f)\n", + "print '%s has %d significant digits\\n'%(f,nf)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 4_03 Pg No. 64" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\n", + " 0.1_10 = 0.[0.0, 0.0, 0.0, 1.0, 1.0, 0.0, 0.0, 1.0] \n", + " 0.4_10 = 0.[0.0, 1.0, 1.0, 0.0, 0.0, 1.0, 1.0, 0.0] \n", + " \n", + "sum = 0.9921875\n", + "Note : The answer should be 0.5, but it is not so.This is due to the error in conversion from decimal to binary form.\n" + ] + } + ], + "source": [ + "from __future__ import division\n", + "from math import floor\n", + "from numpy import arange,zeros\n", + "a = 0.1\n", + "b = 0.4\n", + "afrac=[];bfrac=[];\n", + "for i in range(0,8):\n", + " afrac.append(floor(a*2))\n", + " a = a*2 - floor(a*2)\n", + " bfrac.append(floor(b*2))\n", + " b = b*2 - floor(b*2)\n", + "\n", + "afrac_s = '0' + '.' +(str(afrac)) #string form binary equivalent of a i.e 0.1\n", + "bfrac_s = '0' + '.' +(str(bfrac))\n", + "print '\\n 0.1_10 = %s \\n 0.4_10 = %s \\n '%( afrac_s , bfrac_s)\n", + " \n", + "summ=zeros(8)\n", + "for j in arange(7,0,-1):\n", + " summ[j] = afrac[j] + bfrac[j]\n", + " if summ[j] > 1:\n", + " summ[j] = summ[j]-2\n", + " afrac[(j-1)] = afrac[(j-1)] + 1\n", + "summ_dec = 0\n", + "for k in arange(7,0,-1):\n", + " summ_dec = summ_dec + summ[k]\n", + " summ_dec = summ_dec*1.0/2 \n", + "\n", + "print 'sum =',summ_dec\n", + "print 'Note : The answer should be 0.5, but it is not so.This is due to the error in conversion from decimal to binary form.' " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 4_04 Pg No. 66" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\n", + " Chooping Method : \n", + " Approximate x = 0.7526*10**3 \n", + " Error = 0.0835 \n", + " \n", + "\n", + " Symmetric Rounding :\n", + " Approximate x = 0.7527*10**3 \n", + " Error = -0.0165 \n", + " \n" + ] + } + ], + "source": [ + "#Rounding-Off\n", + "\n", + "fx = 0.7526\n", + "E =3\n", + "gx = 0.835\n", + "d = E - (-1)\n", + "#Chopping Method\n", + "Approx_x = fx*10**E\n", + "Err = gx*10**(E-d)\n", + "print '\\n Chooping Method : \\n Approximate x = %.4f*10**%d \\n Error = %.4f \\n '%(fx,E,Err)\n", + "#Symmetric Method\n", + "if gx >= 0.5:\n", + " Err = (gx -1)*10**(-1)\n", + " Approx_x = (fx + 10**(-d))*10**E\n", + "else:\n", + " Approx_x = fx*10**E\n", + " Err = gx * 10**(E-d)\n", + "\n", + "print '\\n Symmetric Rounding :\\n Approximate x = %.4f*10**%d \\n Error = %.4f \\n '%(fx + 10**(-d),E,Err)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 4_05 Pg No. 68" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\n", + " a) When first three terms are used \n", + " Truncation error = 1.402756E-03 \n", + " \n", + "\n", + " b) When first four terms are used \n", + " Truncation error = 6.942222E-05 \n", + " \n", + "\n", + " c) When first five terms are used \n", + " Truncation error = 2.755556E-06 \n", + " \n" + ] + } + ], + "source": [ + "from scipy.misc import factorial\n", + "#Truncation Error\n", + "\n", + "x = 1.0/5\n", + "#When first three terms are used\n", + "Trunc_err = x**3/factorial(3) + x**4/factorial(4) + x**5/factorial(5) + x**6/factorial(6)\n", + "print '\\n a) When first three terms are used \\n Truncation error = %.6E \\n '%(Trunc_err)\n", + "\n", + "#When four terms are used\n", + "Trunc_err = x**4/factorial(4) + x**5/factorial(5) + x**6/factorial(6)\n", + "print '\\n b) When first four terms are used \\n Truncation error = %.6E \\n '%(Trunc_err)\n", + "\n", + "#When Five terms are used\n", + "Trunc_err = x**5/factorial(5) + x**6/factorial(6)\n", + "print '\\n c) When first five terms are used \\n Truncation error = %.6E \\n '%(Trunc_err)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 4_06 Pg No. 68" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\n", + " a) When first three terms are used \n", + " Truncation error = -1.269244E-03 \n", + " \n", + "\n", + " b) When first four terms are used \n", + " Truncation error = 6.408889E-05 \n", + " \n", + "\n", + " c) When first five terms are used \n", + " Truncation error = -2.577778E-06 \n", + " \n" + ] + } + ], + "source": [ + "from scipy.misc import factorial\n", + "\n", + "#Truncation Error\n", + "\n", + "x = -1.0/5\n", + "#When first three terms are used\n", + "Trunc_err = x**3/factorial(3) + x**4/factorial(4) + x**5/factorial(5) + x**6/factorial(6)\n", + "print '\\n a) When first three terms are used \\n Truncation error = %.6E \\n '%(Trunc_err)\n", + "\n", + "#When four terms are used\n", + "Trunc_err = x**4/factorial(4) + x**5/factorial(5) + x**6/factorial(6)\n", + "print '\\n b) When first four terms are used \\n Truncation error = %.6E \\n '%(Trunc_err)\n", + "\n", + "#When Five terms are used\n", + "Trunc_err = x**5/factorial(5) + x**6/factorial(6)\n", + "print '\\n c) When first five terms are used \\n Truncation error = %.6E \\n '%(Trunc_err)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 4_07 Pg No. 71" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\n", + " For Building : \n", + " Absolute error, e1 = 5 \n", + " Relative error , e1_r = 0.17 percent \n", + " \n", + "\n", + " For Beam : \n", + " Absolute error, e2 = 5 \n", + " Relative error , e2_r = 17 percent \n", + " \n" + ] + } + ], + "source": [ + "#Absolute and Relative Errors\n", + "\n", + "h_bu_t = 2945 #\n", + "h_bu_a = 2950 #\n", + "h_be_t = 30 #\n", + "h_be_a = 35 #\n", + "e1 = abs(h_bu_t - h_bu_a)\n", + "e1_r = e1/h_bu_t\n", + "e2 = abs(h_be_t - h_be_a)\n", + "e2_r = e2/h_be_t\n", + "print '\\n For Building : \\n Absolute error, e1 = %d \\n Relative error , e1_r = %.2f percent \\n '%(e1,e1_r*100) \n", + "print '\\n For Beam : \\n Absolute error, e2 = %d \\n Relative error , e2_r = %.2G percent \\n '%(e2,e2_r*100)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 4_08 Pg No. 72" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "q = 7.9 \n", + " We can say that the computer can store numbers with 7 significant decimal digits \n", + " \n" + ] + } + ], + "source": [ + "from math import log10\n", + "#Machine Epsilon\n", + "\n", + "def Q(p):\n", + " q = 1 + (p-1)*log10(2)\n", + " return q\n", + "p = 24\n", + "q = Q(p)\n", + "print 'q = %.1f \\n We can say that the computer can store numbers with %d significant decimal digits \\n '%(q,q)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 4_09 Pg No. 75" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\n", + " |er_x| <= 0.05 o/o\n", + " |er_y| <= 0.05o/o \n", + " ex = 0.617 \n", + " ey = 0.616 \n", + " |ez| = 1.233 \n", + " |er_z| = 61.65o/o \n", + "\n" + ] + } + ], + "source": [ + "#Propagation of Error\n", + "\n", + "x = 0.1234*10**4\n", + "y = 0.1232*10**4\n", + "d = 4\n", + "er_x = 10**(-d + 1)/2\n", + "er_y = 10**(-d + 1)/2\n", + "ex = x*er_x\n", + "ey = y*er_y\n", + "ez = abs(ex) + abs(ey)\n", + "er_z = abs(ez)/abs(x-y)\n", + "\n", + "print '\\n |er_x| <= %.2f o/o\\n |er_y| <= %.2fo/o \\n ex = %.3f \\n ey = %.3f \\n |ez| = %.3f \\n |er_z| = %.2fo/o \\n'%(er_x *100,er_y*100,ex,ey,ez,er_z*100)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 4_10 Pg No. 77" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\n", + " ex = 0.01175 \n", + " ey = 0.03370 \n", + " ez = 0.01725 \n", + " exy = 0.15839 \n", + " ew = 0.17564 \n", + "\n" + ] + } + ], + "source": [ + "#Errors in Sequence of Computations\n", + "\n", + "x_a = 2.35 #\n", + "y_a = 6.74 #\n", + "z_a = 3.45 #\n", + "ex = abs(x_a)*10**(-3+1)/2\n", + "ey = abs(y_a)*10**(-3+1)/2\n", + "ez = abs(z_a)*10**(-3+1)/2\n", + "exy = abs(x_a)*ey + abs(y_a)*ex\n", + "ew = abs(exy) + abs(ez)\n", + "print '\\n ex = %.5f \\n ey = %.5f \\n ez = %.5f \\n exy = %.5f \\n ew = %.5f \\n'%(ex,ey,ez,exy,ew)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 4_11 Pg No. 77" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "w = 10000.0\n", + "True w = 10444\n", + "ew = 444.0\n", + "er,w = 0.0425124473382\n" + ] + } + ], + "source": [ + "from math import floor\n", + "#Addition of Chain of Numbers\n", + "\n", + "x = 9678 #\n", + "y = 678 #\n", + "z = 78 #\n", + "d = 4 # #length of mantissa\n", + "fx = x/10**4\n", + "fy = y/10**4\n", + "fu = fx + fy\n", + "Eu = 4\n", + "if fu >= 1:\n", + " fu = fu/10\n", + " Eu = Eu + 1\n", + "\n", + "#since length of mantissa is only four we need to maintain only four places in decimal, so\n", + "fu = floor(fu*10**4)/10**4\n", + "u = fu * 10**Eu\n", + "w = u + z\n", + "n = len(str(w))\n", + "w = floor(w/10**(n-4))*10**(n-4) #To maintain length of mantissa = 4\n", + "print 'w =',w\n", + "True_w = 10444\n", + "ew = True_w - w\n", + "er_w = (True_w - w)/True_w\n", + "print 'True w =',True_w\n", + "print 'ew =',ew\n", + "print 'er,w =',er_w " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 4_12 Pg No. 77" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "w = 10430.0\n", + "True w = 10444\n", + "ew = 14.0\n", + "er,w = 0.00134048257373\n" + ] + } + ], + "source": [ + "#Addition of chain Numbers\n", + "\n", + "x = 9678 #\n", + "y = 678 #\n", + "z = 78 #\n", + "d = 4 # #length of mantissa\n", + "n = max(len( str(y) ) , len(str(z)))\n", + "fy = y/10**n\n", + "fz = z/10**n\n", + "fu = fy + fz\n", + "Eu = n\n", + "if fu >= 1:\n", + " fu = fu/10\n", + " Eu = Eu + 1\n", + "\n", + "u = fu * 10**Eu\n", + "n = max(len( str(x) ) , len(str(u)))\n", + "fu = u/10**4\n", + "fx = x/10**4\n", + "fw = fu + fx\n", + "Ew = 4\n", + "if fw >= 1:\n", + " fw = fw/10\n", + " Ew = Ew + 1\n", + "\n", + "#since length of mantissa is only four we need to maintain only four places in decimal, so\n", + "fw = floor(fw*10**4)/10**4\n", + "w = fw*10**Ew\n", + "print 'w =',w\n", + "True_w = 10444\n", + "ew = True_w - w\n", + "er_w = (True_w - w)/True_w\n", + "print 'True w =',True_w\n", + "print 'ew =',ew \n", + "print 'er,w =',er_w" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 4_14 Pg No. 79" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\n", + " ex = 5E-04 \n", + " df(xa) = 1.25 \n", + " ef = 6.26E-04 \n", + " er,f = 1.04E-04 \n", + "\n" + ] + } + ], + "source": [ + "from math import sqrt\n", + "from scipy.misc import derivative\n", + "#Absolute & Relative Errors\n", + "\n", + "xa = 4.000\n", + "def f(x):\n", + " f = sqrt(x) + x\n", + " return f\n", + "#Assuming x is correct to 4 significant digits\n", + "ex = 0.5 * 10**(-4 + 1)\n", + "df_xa = derivative(f,4)\n", + "ef = ex * df_xa\n", + "er_f = ef/f(xa)\n", + "print '\\n ex = %.0E \\n df(xa) = %.2f \\n ef = %.2E \\n er,f = %.2E \\n'%( ex,df_xa,ef,er_f)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 4_15 Pg No. 80" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "ef = 0.07\n" + ] + } + ], + "source": [ + "#Error Evaluation\n", + "\n", + "x = 3.00 #\n", + "y = 4.00 #\n", + "def f(x,y):\n", + " f = x**2 + y**2\n", + " return f\n", + "def df_x(x):\n", + " df_x = 2*x\n", + " return df_x\n", + "def df_y(y):\n", + " df_y = 2*y\n", + " return df_y\n", + "ex = 0.005\n", + "ey = 0.005\n", + "ef = df_x(x)*ex + df_y(y)*ey\n", + "print 'ef =',ef" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 4_16 Pg No. 82" + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "x = 400.0\n", + "y = 807.0\n", + "Changing m2 from 2.01 to 2.005\n", + "\n", + " x = 800 \n", + " y = 1607 \n", + " From the above results we can see that for small change in m2 results in almost 100 percent change in the values of x and y.Therefore, the problem is absolutely ill-conditioned \n", + "\n" + ] + } + ], + "source": [ + "#Condition and Stability\n", + "\n", + "C1 = 7.00 #\n", + "C2 = 3.00 #\n", + "m1 = 2.00 #\n", + "m2 = 2.01 #\n", + "x = (C1 - C2)/(m2 - m1)\n", + "y = m1*((C1 - C2)/(m2 - m1)) + C1\n", + "print 'x =',x\n", + "print 'y =',y\n", + "print 'Changing m2 from 2.01 to 2.005'\n", + "m2 = 2.005\n", + "x = (C1 - C2)/(m2 - m1)\n", + "y = m1*((C1 - C2)/(m2 - m1)) + C1\n", + "print '\\n x = %d \\n y = %d \\n From the above results we can see that for small change in m2 results in almost 100 percent change in the values of x and y.Therefore, the problem is absolutely ill-conditioned \\n'%(x,y)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 4_18 Pg No. 84" + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "z = sqrt(x) - sqrt(y) = 0.0\n", + "z = ( x-y )/( sqrt(x) + sqrt(y) ) = 0.022\n" + ] + } + ], + "source": [ + "from math import sqrt, floor\n", + "#Difference of Square roots\n", + "\n", + "x = 497.0 #\n", + "y = 496.0 #\n", + "sqrt_x = sqrt(497)\n", + "sqrt_y = sqrt(496)\n", + "nx = len( str( floor( sqrt_x ) ) )\n", + "ny = len( str( floor( sqrt_y ) ) )\n", + "sqrt_x = floor(sqrt_x*10**(4-nx))/10**(4-nx)\n", + "sqrt_y = floor(sqrt_y*10**(4-ny))/10**(4-ny)\n", + "z1 = sqrt_x - sqrt_y\n", + "print 'z = sqrt(x) - sqrt(y) =',z1\n", + "z2 = ( x -y)/(sqrt_x + sqrt_y)\n", + "if z2 < 0.1:\n", + " z2 = z2*10**4\n", + " nz = len(str(floor(z2)))\n", + " z2 = floor(z2*10**(4-nz))/10**(8-nz)\n", + "\n", + "print 'z = ( x-y )/( sqrt(x) + sqrt(y) ) =',z2" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example 4_21 Pg No. 85" + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Truncation Error = 0.000505399929331\n", + "calculated e**x using roundoff = -0.000459999999793\n", + "actual e**x = 4.53999297625e-05\n", + "Here we can see the difference between calculated e**x and actual e**x this is due to trucation error (which is greater than final value of e**x ), so the roundoff error totally dominates the solution\n" + ] + } + ], + "source": [ + "from __future__ import division\n", + "from math import exp,floor\n", + "x = -10\n", + "T_act=[1]\n", + "T_trc=[1]\n", + "e_x_cal = 1\n", + "TE=[]\n", + "for i in range(0,100):\n", + " T_act.append(T_act[i]*x/(i+1))\n", + " T_trc.append(floor(T_act[i+1]*10**5)/10**5)\n", + " TE.append(abs(T_act[i+1]-T_trc[i+1]))\n", + " e_x_cal = e_x_cal + T_trc[i+1]\n", + "e_x_act = exp(-10)\n", + "Sum=0\n", + "for s in TE:\n", + " Sum+=s\n", + "print 'Truncation Error =',Sum\n", + "print 'calculated e**x using roundoff =',e_x_cal\n", + "print 'actual e**x = ',e_x_act\n", + "print 'Here we can see the difference between calculated e**x and actual e**x this is due to trucation error (which is greater than final value of e**x ), so the roundoff error totally dominates the solution'" + ] + } + ], + "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 +} diff --git a/Numerical_Methods_by_E._Balaguruswamy/chapter6.ipynb b/Numerical_Methods_by_E._Balaguruswamy/chapter6.ipynb new file mode 100644 index 00000000..76664ab7 --- /dev/null +++ b/Numerical_Methods_by_E._Balaguruswamy/chapter6.ipynb @@ -0,0 +1,1039 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Chapter 6 - Roots of non linear equations" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 6_01 Pg No. 126" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "The largest possible root is x1 = [[4]]\n", + "No root can be larger than the value = [[4]]\n", + "\n", + "All real roots lie in the interval (-3.741657,3.741657)\n", + "\n", + "We can use these two points as initial guesses for the bracketing methods and one of them for open end methods\n" + ] + } + ], + "source": [ + "from numpy import mat,shape\n", + "from math import sqrt\n", + "\n", + "#Possible Initial guess values for roots\n", + "\n", + "A = mat([[ 2],[-8] ,[2],[12]]) # Coefficients of x terms in the decreasing order of power\n", + "n = shape(A)#\n", + "x1 = -A[1]/A[0]#\n", + "print 'The largest possible root is x1 =',x1\n", + "print 'No root can be larger than the value =',x1\n", + "\n", + "x = sqrt((A[1]/A[0])**2 - 2*(A[2]/A[0])**2)\n", + "\n", + "print '\\nAll real roots lie in the interval (-%f,%f)\\n'%(x,x)\n", + "print 'We can use these two points as initial guesses for the bracketing methods and one of them for open end methods'" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 6_03 Pg No. " + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "p(4) = 1\n", + "\n", + " p(3)= -2\n", + "\n", + "\n", + " p(2)= 1\n", + "\n", + "\n", + " f(2) = p(1) = 0\n" + ] + } + ], + "source": [ + "from numpy import mat,shape\n", + "from math import sqrt\n", + "#Evaluating Polynomial using Horner's rule\n", + "\n", + "#Coefficients of x terms in the increasing order of power\n", + "A = mat([[6],[1],[-4],[1]])\n", + "x = 2\n", + "N=shape(A)\n", + "n,c = N[0],N[1]\n", + "p=[0,0,0]\n", + "p.append(A[n-1])\n", + "print 'p(4) =',p[n-1][0,0]\n", + "for i in range(1,n-1):\n", + " p[n-i]= p[n-i+1-1]*x + A[n-i-1]\n", + " print '\\n p(%d)= %d\\n'%(n-i,p[n-i])\n", + "\n", + "print '\\n f(%d) = p(1) = %d'%(x,p[0])\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 6_04 Pg No. 132" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "First finding the interval that contains a root,this can be done by using Eq 6.10\n", + "\n", + " Both the roots lie in the interval (-6,6) \n", + "\n", + "\n", + " The root lies in the interval (5,5)\n", + "\n" + ] + } + ], + "source": [ + "from numpy import poly1d,polyval, sqrt\n", + "\n", + "#Root of a Equation Using Bisection Method\n", + "\n", + "#Coefficients in increasing order of power of x starting from 0\n", + "A = [-10 ,-4, 1]#\n", + "print 'First finding the interval that contains a root,this can be done by using Eq 6.10'\n", + "xmax = sqrt((A[1]/A[2])**2 - 2*(A[0]/A[2]))\n", + "print '\\n Both the roots lie in the interval (-%d,%d) \\n'%(xmax,xmax)\n", + "x = range(-6,7)\n", + "p= poly1d(A)# p = poly(A,'x'%('c'\n", + "\n", + "fx = [p(xx) for xx in x]#\n", + "for i in range(0,12):\n", + " if fx[i]*fx[i] < 0:\n", + " break \n", + " \n", + "\n", + "print '\\n The root lies in the interval (%d,%d)\\n'%(x[i],x[i])\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 6_05 Pg No. 139" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Iteration No. 1 \n", + "\n", + "x0 = 1.000000 \n", + "\n", + "Iteration No. 2 \n", + "\n", + "x0 = 1.000000 \n", + "\n", + "Iteration No. 3 \n", + "\n", + "x0 = 1.000000 \n", + "\n", + "Iteration No. 4 \n", + "\n", + "x0 = 1.000000 \n", + "\n", + "Iteration No. 5 \n", + "\n", + "x0 = 1.000000 \n", + "\n", + "Iteration No. 6 \n", + "\n", + "x0 = 1.000000 \n", + "\n", + "Iteration No. 7 \n", + "\n", + "x0 = 1.000000 \n", + "\n", + "Iteration No. 8 \n", + "\n", + "x0 = 1.000000 \n", + "\n", + "Iteration No. 9 \n", + "\n", + "x0 = 1.000000 \n", + "\n", + "Iteration No. 10 \n", + "\n", + "x0 = 1.000000 \n", + "\n", + "Iteration No. 11 \n", + "\n", + "x0 = 1.000000 \n", + "\n", + "Iteration No. 12 \n", + "\n", + "x0 = 1.000000 \n", + "\n", + "Iteration No. 13 \n", + "\n", + "x0 = 1.000000 \n", + "\n", + "Iteration No. 14 \n", + "\n", + "x0 = 1.000000 \n", + "\n", + "Iteration No. 15 \n", + "\n", + "x0 = 1.000000 \n", + "\n" + ] + } + ], + "source": [ + "from numpy import poly1d,polyval\n", + "\n", + "#False Position Method\n", + "\n", + "#Coefficients of polynomial in increasing order of power of x\n", + "A = [-2 , -1, 1]\n", + "x1 = 1 #\n", + "x2 = 3 #\n", + "fx = poly1d(A)\n", + "for i in range(1,16):\n", + " print 'Iteration No. %d \\n'%(i)\n", + " fx1 = fx(x1)\n", + " fx2 = fx(x2)\n", + " x0 = x1 - fx1*(x2-x1)/(fx2-fx1) \n", + " print 'x0 = %f \\n'%(x0)#\n", + " fx0 = fx(x0)#\n", + " if fx1*fx0 < 0:\n", + " x2 = x0 \n", + " else:\n", + " x1 = x0 #\n", + " " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 6_07 Pg No. 147" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "x2 = 1.000000\n", + "\n", + "Since f(1.000000) = 0, the root closer to the point x = 0 is 1.000000 \n", + "\n" + ] + } + ], + "source": [ + "from numpy import poly1d,polyval,polyder\n", + "#Root of the Equation using Newton Raphson Method\n", + "\n", + "#Coefficients of polynomial in increasing order of power of x\n", + "A = [ 2, -3, 1]#\n", + "fx = poly1d(A)\n", + "dfx = polyder(fx)\n", + "\n", + "x=[0]\n", + "f=[]\n", + "df=[]\n", + "for i in range(1,11):\n", + " f.append(fx(x[i-1]))\n", + " if f[i-1] != 0:\n", + " df.append(dfx(x[i-1]))\n", + " x.append(x[i-1] - f[i-1]/df[i-1])\n", + " print 'x%d = %f\\n'%(i+1,x[i])# \n", + " else:\n", + " print 'Since f(%f) = 0, the root closer to the point x = 0 is %f \\n'%(x[i-1],x[i-1] )\n", + " break\n", + " " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 6_08 Pg No. 151" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "x2 = 3.342105\n", + "\n", + "x3 = 2.248631\n", + "\n", + "x4 = 1.535268\n", + "\n", + "x5 = 1.079139\n", + "\n", + "x6 = 0.797330\n", + "\n", + "x7 = 0.632716\n", + "\n", + "From the results we can see that number of correct digits approximately doubles with each iteration\n" + ] + } + ], + "source": [ + "from numpy import poly1d,polyval,polyder\n", + "\n", + "#Root of the Equation using Newton Raphson Method\n", + "\n", + "#Coefficients of polynomial in increasing order of power of x\n", + "A = [ 6, 1 , -4 , 1 ]\n", + "fx = poly1d(A)\n", + "dfx = polyder(fx)\n", + "f=[];df=[]\n", + "x = [5.0] #\n", + "for i in range(1,7):\n", + " f.append(fx(x[i-1]))\n", + " if f[i-1] != 0:\n", + " df.append(dfx(x[i-1]))\n", + " x.append(x[i-1] - f[i-1]/df[i-1])\n", + " print 'x%d = %f\\n'%(i+1,x[i])\n", + " \n", + "\n", + "print 'From the results we can see that number of correct digits approximately doubles with each iteration'" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 6_09 Pg No. 153" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\n", + " For Iteration No. 1\n", + "\n", + "\n", + " x1 = 4.000000\n", + " x2 = 2.000000 \n", + " fx1 = -175.000000 \n", + " fx2 = -47.000000 \n", + " x3 = 1.265625 \n", + "\n", + "\n", + " For Iteration No. 2\n", + "\n", + "\n", + " x1 = 2.000000\n", + " x2 = 1.265625 \n", + " fx1 = -47.000000 \n", + " fx2 = -20.080566 \n", + " x3 = 0.717818 \n", + "\n", + "\n", + " For Iteration No. 3\n", + "\n", + "\n", + " x1 = 1.265625\n", + " x2 = 0.717818 \n", + " fx1 = -20.080566 \n", + " fx2 = -7.023891 \n", + " x3 = 0.423122 \n", + "\n", + "\n", + " For Iteration No. 4\n", + "\n", + "\n", + " x1 = 0.717818\n", + " x2 = 0.423122 \n", + " fx1 = -7.023891 \n", + " fx2 = -2.482815 \n", + " x3 = 0.261999 \n", + "\n", + "\n", + " For Iteration No. 5\n", + "\n", + "\n", + " x1 = 0.423122\n", + " x2 = 0.261999 \n", + " fx1 = -2.482815 \n", + " fx2 = -0.734430 \n", + " x3 = 0.194317 \n", + "\n", + "\n", + " For Iteration No. 6\n", + "\n", + "\n", + " x1 = 0.261999\n", + " x2 = 0.194317 \n", + " fx1 = -0.734430 \n", + " fx2 = -0.154860 \n", + " x3 = 0.176233 \n", + "\n", + "This can be still continued further for accuracy\n" + ] + } + ], + "source": [ + "from __future__ import division\n", + "from numpy import poly1d\n", + "#Root of the equation using SECANT Method\n", + "\n", + "#Coefficients of polynomial in increasing order of power of x\n", + "A = [ -10, -4, 1]\n", + "x1 = 4 #\n", + "x2 = 2 #\n", + "fx = poly1d(A)\n", + "for i in range(1,7):\n", + " print '\\n For Iteration No. %d\\n'%(i)\n", + " fx1 = fx(x1)\n", + " fx2 = fx(x2)\n", + " x3 = x2 - fx2*(x2-x1)/(fx2-fx1) #\n", + " print '\\n x1 = %f\\n x2 = %f \\n fx1 = %f \\n fx2 = %f \\n x3 = %f \\n'%(x1,x2,fx1,fx2,x3) #\n", + " x1 = x2#\n", + " x2 = x3#\n", + "\n", + "print 'This can be still continued further for accuracy'" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 6_11 Pg No. 161" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\n", + " x1 = -1.000000\n", + "\n", + "\n", + " x2 = 1.000000\n", + "\n", + "\n", + " x3 = 1.000000\n", + "\n", + "\n", + "1.000000 is root of the equation,since x3 - x2 = 0 \n", + "\n", + "\n", + "x1 = 1.000000\n", + "\n", + "\n", + "x2 = 1.000000\n", + "\n", + "\n", + " 1.000000 is root of the equation,since x2 - x1 = 0\n" + ] + } + ], + "source": [ + "from numpy import poly1d\n", + "from __future__ import division\n", + "#Fixed point method\n", + "\n", + "#Coefficients of polynomial in increasing order of power of x\n", + "A = [ -2, 1, 1 ]#\n", + "B = [ 2, 0, -1 ]#\n", + "gx = poly1d(B)\n", + "x = [0] ##initial guess x0 = 0\n", + "for i in range(2,11):\n", + " x.append (gx(x[i-2]))\n", + " print '\\n x%d = %f\\n'%(i-1,x[i-1])\n", + " if (x[i-1]-x[(i-2)]) == 0:\n", + " print '\\n%f is root of the equation,since x%d - x%d = 0 \\n'%(x[i-1],i-1,i-2)\n", + " break\n", + " \n", + "\n", + "#Changing initial guess x0 = -1\n", + "x[0] = -1 #\n", + "for i in range(2,11):\n", + " x[i-1]= gx(x[i-2])\n", + " print '\\nx%d = %f\\n'%(i-1,x[i-1])\n", + " if (x[i-1]-x[i-2]) == 0:\n", + " print '\\n %f is root of the equation,since x%d - x%d = 0'%(x[i-1],i-1,i-2)\n", + " break" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 6_12 Pg No. 162" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\n", + " x0 = 1.000000 \n", + "\n", + " x1 = 2.500000 \n", + "\n", + " x2 = 1.666667 \n", + "\n", + " x3 = 1.250000 \n", + "\n", + " x4 = 1.000000 \n", + "\n", + "\n", + " x0 = 0.000000 \n", + "\n", + " x1 = 1.000000 \n", + "\n", + " x2 = 7.000000 \n", + "\n", + " x3 = 15.000000 \n", + "\n", + " x4 = 25.000000 \n", + "\n", + "\n", + " x0 = 1.000000 \n", + "\n", + " x1 = 2.250000 \n", + "\n", + " x2 = 2.333333 \n", + "\n", + " x3 = 2.625000 \n", + "\n", + " x4 = 3.000000 \n", + "\n", + " x5 = 3.416667 \n", + "\n", + " x6 = 3.857143 \n", + "\n" + ] + } + ], + "source": [ + "#Fixed point method\n", + "\n", + "A = [ -5, 0, 1 ]#\n", + "def g(x):\n", + " x = 5.0/x\n", + " return x\n", + "x = [1] #\n", + "print '\\n x0 = %f \\n'%(x[0])\n", + "for i in range(2,6):\n", + " x.append(g(i))\n", + " print ' x%d = %f \\n'%(i-1,x[i-1])\n", + " \n", + "\n", + "#Defining g(x) in different way\n", + "def g(x):\n", + " x = x**2 + x - 5\n", + " return x\n", + "x=[0]\n", + "print '\\n x0 = %f \\n'%(x[0])\n", + "for i in range(2,6):\n", + " x.append(g(i))\n", + " print ' x%d = %f \\n'%(i-1,x[i-1])\n", + "\n", + "#Third form of g(x)\n", + "def g(x):\n", + " x = (x + 5/x)/2\n", + " return x\n", + "x=[1]\n", + "print '\\n x0 = %f \\n'%(x[0])\n", + "for i in range(2,8):\n", + " x.append(g(i))\n", + " print ' x%d = %f \\n'%(i-1,x[i-1])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 6_13 Pg No. 169" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + " x**2 - y**2 = 3 \n", + " x**2 + x*y \n", + "\n", + "\n", + " x0 = 1.000000 \n", + " y0 = 1.000000 \n", + "\n", + "\n", + " x1 = 2.500000 \n", + " y1 = 5.000000 \n", + "\n", + "\n", + " x2 = 2.750000 \n", + " y2 = 1.000000 \n", + "\n", + "\n", + " x3 = 3.500000 \n", + " y3 = -1.000000 \n", + "\n" + ] + } + ], + "source": [ + "#Solving System of non-linear equations using FIXED POINT METHOD\n", + "\n", + "print ' x**2 - y**2 = 3 \\n x**2 + x*y \\n'\n", + "def f(x,y):\n", + " x = y + 3/(x+y)\n", + " return x\n", + "def g(x):\n", + " y = (6-x**2)/x\n", + " return y\n", + "x=[1]\n", + "y=[1]\n", + "print '\\n x0 = %f \\n y0 = %f \\n'%(x[0],y[0])\n", + "for i in range(2,5):\n", + " x.append(f((i-1),(i-1)))\n", + " y.append(g((i-1)))\n", + " print '\\n x%d = %f \\n y%d = %f \\n'%(i-1,x[i-1],i-1,y[i-1])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 6_14 Pg No. 172" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "x**2 + x*y = 6 \n", + " x**2 - y**2 = 3 \n", + "\n", + "\n", + " x1 = 0.875000 \n", + " y1 = -0.625000 \n", + "\n", + "\n", + " x2 = -0.437500 \n", + " y2 = -2.687500 \n", + "\n" + ] + } + ], + "source": [ + "#Solving System of Non-linear equations using Newton Raphson Method\n", + "\n", + "\n", + "print 'x**2 + x*y = 6 \\n x**2 - y**2 = 3 \\n'#\n", + "def F(x,y):\n", + " f = x**2 + x*y - 6\n", + " return f\n", + "def G(x,y):\n", + " g = x**2 - y**2 -3\n", + " return g\n", + "def dFx(x,y):\n", + " f1 = 2*x + y\n", + " return f1\n", + "def dFy(x,y):\n", + " f2 = y\n", + " return f2\n", + "def dGx(x,y):\n", + " g1 = 2*x\n", + " return g1\n", + "def dGy(x,y):\n", + " g2 = -2*y\n", + " return g2\n", + "x=[1]\n", + "y=[1]\n", + "\n", + "for i in range(2,4):\n", + " Fval = F(i,i)\n", + " Gval = G(i,i)\n", + " f1 = dFx(i-1,i-1)\n", + " f2 = dFy(i-1,i-1)\n", + " g1 = dGx(i-1,i-1)\n", + " g2 = dGy(i-1,i-1)\n", + " D = f1*g2 - f2*g1 \n", + " \n", + " x.append(x[i-2] - (Fval*g2 - Gval*f2)/D )\n", + " y.append(y[i-2] - (Gval*f1 - Fval*g1)/D )\n", + " print '\\n x%d = %f \\n y%d = %f \\n'%(i-1,x[i-1],i-1,y[i-1]) " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 6_15 Pg No. 176" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "b3 = 0.000000\n", + "\n", + "b2 = 0.000000\n", + "\n", + "b1 = 0.000000\n", + "\n", + "Thus the polynomial is\n", + " 2\n", + "15 x - 7 x + 1\n" + ] + } + ], + "source": [ + "from __future__ import division\n", + "from numpy import poly1d, arange\n", + "#Synthetic Division\n", + "\n", + "a = [-9 ,15 ,-7, 1]\n", + "b=[0,0,0,0]\n", + "for i in arange(3,0,-1):\n", + " b[i]= a[i] + b[i]*3\n", + " print 'b%d = %f\\n'%(i,b[i-1])\n", + "print 'Thus the polynomial is' \n", + "print poly1d(b)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 6_16 Pg No. 187" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "b3 = 1.000000 \n", + "\n", + "b2 = 1.800000 \n", + "\n", + "b1 = 3.240000 \n", + "\n", + "b0 = 14.032000 \n", + "\n", + "c3 = 1.000000 \n", + "\n", + "c2 = 1.800000 \n", + "\n", + "c1 = 3.240000 \n", + "\n", + "c0 = 14.032000 \n", + "\n", + "\n", + " D = 4.240000 \n", + " du = -0.694340 \n", + " dv = -5.250566 \n", + " u = 1.105660\n", + " v = -1.694340 \n", + "\n" + ] + } + ], + "source": [ + "#Quadratic factor of a polynomial using Bairstow's Method\n", + "\n", + "a = [ 10, 1 ,0 ,1]#\n", + "n = len(a)#\n", + "u = 1.8 #\n", + "v = -1 #\n", + "b=[];c=[]\n", + "for nn in range(n):\n", + " b.append(0)\n", + " c.append(0)\n", + "b[n-1] = a[n-1]\n", + "b[n-2] = a[n-2] + u*b[n-1]\n", + "c[n-1] = 0 \n", + "c[n-2] = b[n-1]\n", + "\n", + "\n", + "for i in range(n-2,0,-1):\n", + " b[i-1]= a[i-1]+ u*b[i] + v*b[i+1]\n", + " c[i-1]= b[i] + u*c[i] + v*c[i+1] \n", + "\n", + "\n", + "for i in range(n,0,-1):\n", + " print 'b%d = %f \\n'%(i-1,b[i-1])\n", + "\n", + "for i in range(n,0,-1):\n", + " print 'c%d = %f \\n'%(i-1,b[i-1])\n", + "\n", + "\n", + "D = c[1]*c[1] - c[0]*c[2]\n", + "du = -1*(b[1]*c[1] - c[0]*c[2])/D \n", + "dv = -1*(b[0]*c[1] - b[1]*c[0])/D \n", + "u = u + du #\n", + "v = v + du #\n", + "print '\\n D = %f \\n du = %f \\n dv = %f \\n u = %f\\n v = %f \\n'%(D,du,dv,u,v)\n", + " " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 6_17 Pg No. 197" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\n", + " x1 = 0.000000\n", + " x2 = 1.000000\n", + " x3 = 2.000000\n", + " f1 = -20.000000\n", + " f2 = -7.000000\n", + " f3 = 16.000000\n", + " h1 = -2.000000\n", + " h2 = -1.000000\n", + " d1 = -36.000000\n", + " d2 = -23.000000\n", + " a0 = 16.000000\n", + " a1 = 28.000000\n", + " a2 = 5.000000\n", + " h4 = -0.645934\n", + " x4 = 1.354066\n", + " \n", + "\n", + " x1 = 1.000000\n", + " x2 = 2.000000\n", + " x3 = 1.354066\n", + " f1 = -7.000000\n", + " f2 = 16.000000\n", + " f3 = -0.309679\n", + " h1 = -0.354066\n", + " h2 = 0.645934\n", + " d1 = -6.690321\n", + " d2 = 16.309679\n", + " a0 = -0.309679\n", + " a1 = 21.145451\n", + " a2 = 6.354066\n", + " h4 = 0.014581\n", + " x4 = 1.368647\n", + " \n", + "\n", + " x1 = 2.000000\n", + " x2 = 1.354066\n", + " x3 = 1.368647\n", + " f1 = 16.000000\n", + " f2 = -0.309679\n", + " f3 = -0.003394\n", + " h1 = 0.631353\n", + " h2 = -0.014581\n", + " d1 = 16.003394\n", + " d2 = -0.306286\n", + " a0 = -0.003394\n", + " a1 = 21.103381\n", + " a2 = 6.722713\n", + " h4 = 0.000161\n", + " x4 = 1.368808\n", + " \n", + "\n", + " x1 = 1.354066\n", + " x2 = 1.368647\n", + " x3 = 1.368808\n", + " f1 = -0.309679\n", + " f2 = -0.003394\n", + " f3 = -0.000001\n", + " h1 = -0.014742\n", + " h2 = -0.000161\n", + " d1 = -0.309678\n", + " d2 = -0.003392\n", + " a0 = -0.000001\n", + " a1 = 21.096136\n", + " a2 = 6.091521\n", + " h4 = 0.000000\n", + " x4 = 1.368808\n", + " \n", + "root of the polynomial is x4 = 1.368808\n" + ] + } + ], + "source": [ + "from math import sqrt\n", + "#Solving Leonard's equation using MULLER'S Method\n", + "\n", + "def f(x):\n", + " y = x**3 + 2*x**2 + 10*x - 20\n", + " return y\n", + "x1 = 0 #\n", + "x2 = 1 #\n", + "x3 = 2 #\n", + "for i in range(1,11):\n", + " f1 = f(x1)\n", + " f2 = f(x2)\n", + " f3 = f(x3)\n", + " h1 = x1-x3 #\n", + " h2 = x2-x3 #\n", + " d1 = f1 - f3 #\n", + " d2 = f2 - f3 #\n", + " D = h1*h2*(h1-h2)#\n", + " a0 = f3 #\n", + " a1 = (d2*h1**2 - d1*h2**2)/D #\n", + " a2 = (d1*h2 - d2*h1)/D #\n", + " if abs(-2*a0/( a1 + sqrt( a1**2 - 4*a0*a2 ) )) < abs( -2*a0/( a1 - sqrt( a1**2 - 4*a0*a2 ) )):\n", + " h4 = -2*a0/(a1 + sqrt(a1**2 - 4*a0*a2))#\n", + " else:\n", + " h4 = -2*a0/(a1 - sqrt(a1**2 - 4*a0*a2))\n", + " \n", + " x4 = x3 + h4 #\n", + " print '\\n x1 = %f\\n x2 = %f\\n x3 = %f\\n f1 = %f\\n f2 = %f\\n f3 = %f\\n h1 = %f\\n h2 = %f\\n d1 = %f\\n d2 = %f\\n a0 = %f\\n a1 = %f\\n a2 = %f\\n h4 = %f\\n x4 = %f\\n '%(x1,x2,x3,f1,f2,f3,h1,h2,d1,d2,a0,a1,a2,h4,x4) #\n", + " relerr = abs((x4-x3)/x4)#\n", + " if relerr <= 0.00001:\n", + " print 'root of the polynomial is x4 = %f'%(x4)\n", + " break\n", + " \n", + " x1 = x2 #\n", + " x2 = x3 #\n", + " x3 = x4 #" + ] + } + ], + "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 +} 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 +} diff --git a/Numerical_Methods_by_E._Balaguruswamy/chapter8.ipynb b/Numerical_Methods_by_E._Balaguruswamy/chapter8.ipynb new file mode 100644 index 00000000..12e74ff2 --- /dev/null +++ b/Numerical_Methods_by_E._Balaguruswamy/chapter8.ipynb @@ -0,0 +1,317 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Chapter 8 - Iterative solutions of linear equations" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 8_01 Page No. 254" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "x1 = (5 - x2 - x3)/2 \n", + "x2 = (15 - 3x1 - 2x3)/5 \n", + "x3 = (8 - 2x1 - x2)/4\n", + "\n", + " Iteration Number : 1\n", + "\n", + " \n", + " x1 = 2.500000\n", + " x2 = 3.000000\n", + " x3 = 2.000000\n", + "\n", + "\n", + " Iteration Number : 2\n", + "\n", + " \n", + " x1 = 0.000000\n", + " x2 = 0.700000\n", + " x3 = 0.000000\n", + "\n", + "\n", + " Iteration Number : 3\n", + "\n", + " \n", + " x1 = 2.150000\n", + " x2 = 3.000000\n", + " x3 = 1.825000\n", + "\n", + "\n", + " Iteration Number : 4\n", + "\n", + " \n", + " x1 = 0.087500\n", + " x2 = 0.980000\n", + " x3 = 0.175000\n", + "\n" + ] + } + ], + "source": [ + "from __future__ import division\n", + "from numpy import mat\n", + "#Gauss Jacobi\n", + "\n", + "A = mat([[ 2, 1 , 1] ,[3, 5, 2],[2, 1, 4]])\n", + "B = mat([[ 5],[ 15],[ 8]])\n", + "x1old = 0 ;x2old = 0 ; x3old = 0 #intial assumption of x1,x2 & x3\n", + "\n", + "print 'x1 = (5 - x2 - x3)/2 '\n", + "print 'x2 = (15 - 3x1 - 2x3)/5 '\n", + "print 'x3 = (8 - 2x1 - x2)/4'\n", + "\n", + "for i in range(1,5):\n", + " print '\\n Iteration Number : %d\\n'%(i)\n", + " \n", + " x1 = (5 - x2old - x3old)/2 #\n", + " x2 = (15 - 3*x1old - 2*x3old)/5 # \n", + " x3 = (8 - 2*x1old - x2old)/4 #\n", + " \n", + " print ' \\n x1 = %f\\n x2 = %f\\n x3 = %f\\n'%(x1,x2,x3)\n", + " \n", + " x1old = x1#\n", + " x2old = x2#\n", + " x3old = x3#\n", + " \n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 8_02 Page No.261" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "(x1 = 5 - x2 - x3)/2 \n", + "(x2 = 15 - 3x1 - 2x3)/5 \n", + "(x3 = 8 - 2x1 - x2)/4\n", + "\n", + " Iteration Number : 1\n", + " \n", + " x1 = 2.500000\n", + " x2 = 1.500000\n", + " x3 = 0.375000\n", + "\n", + "\n", + " Iteration Number : 2\n", + " \n", + " x1 = 1.562500\n", + " x2 = 1.912500\n", + " x3 = 0.740625\n", + "\n" + ] + } + ], + "source": [ + "from __future__ import division\n", + "from numpy import mat\n", + "#Gauss Seidel\n", + "\n", + "A = mat([[ 2, 1 , 1] ,[3, 5, 2],[2, 1, 4]])\n", + "B = mat([[ 5],[ 15],[ 8]])\n", + "\n", + "x1old = 0 ;x2old = 0 ; x3old = 0 #intial assumption\n", + "\n", + "print '(x1 = 5 - x2 - x3)/2 '\n", + "print '(x2 = 15 - 3x1 - 2x3)/5 '\n", + "print '(x3 = 8 - 2x1 - x2)/4'\n", + " \n", + "for i in range(1,3):\n", + " \n", + " print '\\n Iteration Number : %d'%(i)\n", + " \n", + " x1 = (5 - x2old - x3old)/2 #\n", + " x1old = x1# \n", + " x2 = (15 - 3*x1old - 2*x3old)/5 #\n", + " x2old = x2# \n", + " x3 = (8 - 2*x1old - x2old)/4 #\n", + " x3old = x3#\n", + " \n", + " print ' \\n x1 = %f\\n x2 = %f\\n x3 = %f\\n'%(x1,x2,x3)\n", + " \n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 8_03 page no. 269" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Using a matrix to display the results after each iteration, first row represents initial assumption\n", + "x1 = (5-x2)/3\n", + "x2 = (x1 - 5)/3\n", + "The system converges to the solution ( 1.999949 , -1.000017 ) in 4 iterations\n", + "\n", + "Final Result is as : \n", + "0.0000000000 \t0.0000000000 \t1.9999491947 \t1.00001693509\n", + "-1.1111111111 \t1.6666666667 \t0.3332825281 \t0.111094176023\n", + "-0.9876543210 \t2.0370370370 \t0.0370878423 \t0.0123626141002\n", + "-1.0013717421 \t1.9958847737 \t0.0040644211 \t0.00135480702467\n", + "-0.9998475842 \t2.0004572474 \t0.0005080526 \t0.000169350878084\n" + ] + } + ], + "source": [ + "from __future__ import division\n", + "from numpy import array,zeros,ones,nditer,vstack,hstack\n", + "\n", + "#Gauss Seidel\n", + "\n", + "\n", + "A = array([[ 3, 1],[ 1 ,-3]])\n", + "B = array([[ 5],[5 ]])\n", + "\n", + "X=zeros([6,2])\n", + "print ('Using a matrix to display the results after each iteration, first row represents initial assumption')\n", + "X[0,0] = 0; X[0,1] = 0 ;#initial assumption\n", + "\n", + "maxit = 1000;#Maximum number of iterations\n", + "err = 0.0003 ;\n", + "\n", + "print('x1 = (5-x2)/3');\n", + "print('x2 = (x1 - 5)/3');\n", + "\n", + "for i in range(1,maxit):\n", + " \n", + " X[i,0] = (5 - X[i-1,1])/3 ;\n", + " X[i,1] = (X[i,0] - 5)/3 ;\n", + " \n", + " #Error Calculations\n", + " err1 =abs((X[i,0] - X[i-1,0])/X[i,0]) \n", + " err2 =abs((X[i,1]- X[i-1,1])/X[i,1])\n", + " \n", + " #Terminating Condition \n", + " if err >= err1 and err >= err2:\n", + " print 'The system converges to the solution ( %f , %f ) in %d iterations\\n'%(X[i,0],X[i,1],i-1) \n", + " break\n", + " \n", + " \n", + "\n", + "#calcution of true error i.e. difference between final result and results from each iteration\n", + "trueerr1 = abs(X[:,0] - X[i,0]*ones([i+1,1])) ;\n", + "trueerr2 = abs(X[:,1] - X[i,1]*ones([i+1,1])) ;\n", + "\n", + "#displaying final results\n", + "print 'Final Result is as : '\n", + "for i in range(0,5):\n", + " print '%.10f'%X[i,:][1],'\\t%.10f'%X[i,:][0],'\\t%.10f'%trueerr1[0,i],'\\t',trueerr2[0,i]" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 8_04 Page No.261" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "x1 = 5 + 3*x2 \n", + "x2 = 5 - 3*x1 \n", + "\n", + " Iteration : 1 x1 = 5 and x2 = -10\n", + "\n", + "\n", + " Iteration : 2 x1 = -25 and x2 = 80\n", + "\n", + "\n", + " Iteration : 3 x1 = 245 and x2 = -730\n", + "\n", + "It is clear that the process do not converge towards the solution, rather it diverges.\n" + ] + } + ], + "source": [ + "from numpy import mat\n", + "#Gauss Seidel\n", + "\n", + "A = mat([[ 1, -3],[3, 1 ]])\n", + "B = mat([[ 5],[5] ])\n", + "x1old = 0 #intial assumption\n", + "x2old = 0 #intial assumption\n", + "\n", + "print 'x1 = 5 + 3*x2 '\n", + "print 'x2 = 5 - 3*x1 '\n", + " \n", + "for i in range(1,4):\n", + " x1 = 5 + 3*x2old \n", + " x1old = x1\n", + " x2 = 5 - 3*x1old \n", + " x2old = x2\n", + " \n", + " print '\\n Iteration : %d x1 = %d and x2 = %d\\n'%(i,x1,x2)\n", + " \n", + "print 'It is clear that the process do not converge towards the solution, rather it diverges.'" + ] + } + ], + "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 +} diff --git a/Numerical_Methods_by_E._Balaguruswamy/chapter9.ipynb b/Numerical_Methods_by_E._Balaguruswamy/chapter9.ipynb new file mode 100644 index 00000000..539b10b0 --- /dev/null +++ b/Numerical_Methods_by_E._Balaguruswamy/chapter9.ipynb @@ -0,0 +1,798 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Chapter 9 - Curve fitting interpolation" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 9_01 Pg No.277" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "solving linear equations \n", + " a0 + 100a1 = 3/7 \n", + " a0 + 101a1 = -4/7 \n", + " we get\n", + " a0 = 100 \n", + " a1 = -1\n", + "\n", + " p(100) = 0 \n", + " p(101) = -1\n", + "\n" + ] + } + ], + "source": [ + "from numpy.linalg import solve\n", + "from numpy import mat\n", + "from sympy import Symbol\n", + "print 'solving linear equations \\n a0 + 100a1 = 3/7 \\n a0 + 101a1 = -4/7 \\n we get'#\n", + "C = mat([[ 1, 100],[1 ,101] ])\n", + "p = mat([[ 3/7],[-4/7] ])\n", + "a = solve(C,p )\n", + "print ' a0 = %.f \\n a1 = %.f'%(a[0],a[1])\n", + "\n", + "x = Symbol('x')\n", + "def horner(a,x):\n", + " px = a[0] + a[1]*x\n", + " return px\n", + "p100 = horner(a,100.0)\n", + "p101 = horner(a,101.0)\n", + "print '\\n p(100) = %.f \\n p(101) = %.f\\n'%(p100,p101)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 9_02 Page No. 278" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\n", + " a0 = 0.428571 \n", + " a1 = -1.428571 \n", + "\n", + "\n", + " p(100) = 0.428571 \n", + " p(101) = -1.000000\n", + "\n" + ] + } + ], + "source": [ + "from numpy.linalg import solve\n", + "from numpy import mat\n", + "from sympy import Symbol\n", + "\n", + "C = mat([[1, 100-100],[1, 101-100]])\n", + "p = mat([[ 3.0/7],[-4/7] ])\n", + "a = solve(C,p)\n", + "print '\\n a0 = %f \\n a1 = %f \\n'%(a[0],a[1])\n", + "\n", + "x = Symbol('x')\n", + "def horner(a,x):\n", + " px = a[0]+ a[1]*(x - 100) \n", + " return px\n", + "p100 = horner(a,100)\n", + "p101 = horner(a,101)\n", + "print '\\n p(100) = %f \\n p(101) = %f\\n'%(p100,p101)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 9_03 Page No. 280" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "2.5 lies between points 2 and 3\n", + "f(2.5) = 1.57315\n", + "error1 = 0.00795\n", + "The correct answer is 1.5811.The difference between results is due to use of a linear model to a nonlinear use\n", + "repeating the procedure using x1 = 2 and x2 = 4\n", + "error2 = 0.02045\n", + "f(2.5) = 1.56065\n", + "NOTE- The increase in error due to the increase in the interval between the interpolating data\n" + ] + } + ], + "source": [ + "x = range(0,6)\n", + "f = [0,1, 1.4142, 1.7321, 2, 2.2361]\n", + "n = 2.5\n", + "for i in range(1,6):\n", + " if n <= x[(i)]:\n", + " break;\n", + " \n", + "\n", + "print '%.1f lies between points %d and %d'%(n,x[(i-1)],x[(i)])\n", + "f2_5 = f[(i-1)] + ( n - x[(i-1)] )*( f[(i)] - f[(i-1)] )/( x[(i)] - x[(i-1)] )\n", + "err1 = 1.5811 - f2_5\n", + "print 'f(2.5) = ',f2_5\n", + "print 'error1 = ',err1\n", + "print 'The correct answer is 1.5811.The difference between results is due to use of a linear model to a nonlinear use'\n", + "print 'repeating the procedure using x1 = 2 and x2 = 4'\n", + "x1 = 2\n", + "x2 = 4\n", + "f2_5 = f[(x1)] + ( 2.5 - x1 )*( f[(x2)] - f[(x1)] )/( x2 - x1 )\n", + "err2 = 1.5811 - f2_5\n", + "print 'error2 = ',err2\n", + "print 'f(2.5) = ',f2_5\n", + "print 'NOTE- The increase in error due to the increase in the interval between the interpolating data'" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 9_04 Pg No. 282" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "For x = 2.5 we have,\n", + " L0(2.5) = 1.500000 \n", + " L1(2.5) = 0.250000 \n", + " L2(2.5) = 1.000000 \n", + " p(2.5) = 8.893756 \n", + "\n", + "The error is -7.312617\n" + ] + } + ], + "source": [ + "from sympy import Symbol\n", + "from math import sqrt\n", + "#Lagrange Interpolation- Second order\n", + "\n", + "X = [0, 1, 2 ,3 ,4 ,5]\n", + "Fx = [0, 1, 1.4142 ,1.7321 ,2, 2.2361]\n", + "X = X[2:5]\n", + "Fx = Fx[2:5]\n", + "x0 = 2.5 \n", + "x = Symbol('x')\n", + "p = 0\n", + "L=[0]\n", + "def horner(X,p,Fx,x,m):\n", + " for i in range(1,3):\n", + " L.append(1)\n", + " for j in range(1,3):\n", + " if j == i:\n", + " continue #\n", + " else:\n", + " L[(i)] = L[(i)]*( x - X[(j)] )/( X[(i)] - X[(j)] ) \n", + " \n", + " p = p + Fx[(i)]*L[(i)] \n", + " return [L[m],p]\n", + "\n", + "x=2.5\n", + "L0 = horner(X,p,Fx,x,1)[0]\n", + "L1 = horner(X,p,Fx,x,2)[0]\n", + "L2 = horner(X,p,Fx,x,3)[0]\n", + "p2_5 = horner(X,p,Fx,x,3)[1]\n", + "print 'For x = 2.5 we have,\\n L0(2.5) = %f \\n L1(2.5) = %f \\n L2(2.5) = %f \\n p(2.5) = %f \\n'%(L0,L1,L2,p2_5)\n", + "\n", + "err = sqrt(2.5) - p2_5 #\n", + "print 'The error is %f'%(err)# " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 9_05 Pg No. 283" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "The Lagrange basis polynomials are:\n", + "(-x + 1)*(x - 3)*(x - 2)/6\n", + "x*(x - 3)*(x - 2)/2\n", + "-x*(x - 3)*(x - 1)/2\n", + "x*(x - 2)*(x - 1)/6\n", + "The interpolation polynomial is \n", + "0.85915*x*(x - 3)*(x - 2) - 3.19455*x*(x - 3)*(x - 1) + 3.18091666666667*x*(x - 2)*(x - 1)\n", + "The interpolation value at x = 1.5 is : 4.36756875000000 \n", + "e**1.5 = 3.367569\n" + ] + } + ], + "source": [ + "from sympy import Symbol\n", + "#Lagrange Interpolation\n", + "\n", + "i = [ 0, 1, 2, 3 ]\n", + "X = [ 0 ,1 ,2 ,3 ]\n", + "Fx = [ 0 ,1.7183 ,6.3891, 19.0855 ]\n", + "x = Symbol('x')\n", + "n = 3 #order of lagrange polynomial \n", + "p = 0\n", + "L=[]\n", + "for i in range(0,n+1):\n", + " L.append(1)\n", + " for j in range(0,n+1):\n", + " if j == i:\n", + " continue \n", + " else:\n", + " L[i] = L[i]*( x - X[j] )/( X[i] - X[j] ) \n", + " \n", + " \n", + " p = p + Fx[i]*L[i]\n", + "\n", + "print \"The Lagrange basis polynomials are:\"\n", + "for i in range(0,4):\n", + " print str(L[i])\n", + "\n", + "\n", + "print \"The interpolation polynomial is \"\n", + "print str(p)\n", + "\n", + "print 'The interpolation value at x = 1.5 is :', \n", + "\n", + "p1_5 = p.subs(x,1.5)\n", + "e1_5 = p1_5 + 1 #\n", + "print e1_5,'\\ne**1.5 = %f'%(p1_5)# \n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 9_06 Pg No. 288" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "The coefficients of the polynomial are,\n", + " a0 = 0 \n", + " a1 = 0.301 \n", + " a2 = -0.06245 \n", + "\n", + "Polynomial is : 0.301*x - 0.06245*(1.0*x - 2)*(1.0*x - 1) - 0.301\n", + "p(2.5) = 0.4047 \n", + "\n" + ] + } + ], + "source": [ + "from numpy import zeros, prod, ones, array\n", + "from sympy.abc import x\n", + "i = [ 0, 1 ,2 ,3]\n", + "X = range(1,5)\n", + "Fx = [ 0, 0.3010, 0.4771, 0.6021] \n", + "X = range(1,4)\n", + "Fx = Fx[0:3]\n", + "#x = poly(0,'x');\n", + "#A = Fx'\n", + "A=zeros([3,3])\n", + "A[:,0]=Fx\n", + "for i in range(2,4):\n", + " for j in range(1,4-i+1):\n", + " A[j-1,i-1] = ( A[j+1-1,i-1-1]-A[j-1,i-1-1] )/(X[j+i-1-1]-X[j-1]) ;\n", + "\n", + "print 'The coefficients of the polynomial are,\\n a0 = %.4G \\n a1 = %.4G \\n a2 = %.4G \\n'%(A[0,0],A[0,1],A[0,2])\n", + "p = A[0,0]\n", + "\n", + "for i in range(2,4):\n", + " p = p +A[0,i-1]* prod(array([x*xx for xx in ones([1,i-1])]) - array(X[0:i-1]))\n", + "print 'Polynomial is : ',str(p)\n", + "x=2.5\n", + "p=A[0,0]\n", + "for i in range(2,4):\n", + " p = p +A[0,i-1]* prod(array([x*xx for xx in ones([1,i-1])]) - array(X[0:i-1]))\n", + "print 'p(2.5) = %.4G \\n'%p" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 9_07 Pg No. 291" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + " p0(1.5) = 0.000000 \n", + " p1(1.5) = 3.500000 \n", + " p2(1.5) = 2.000000 \n", + " p3(1.5) = 2.375000 \n", + " p4(1.5) = 2.375000 \n", + "\n", + "The function value at x = 1.5 is : 2.375\n" + ] + } + ], + "source": [ + "from numpy import zeros, prod, ones, array\n", + "#Newton Divided Difference Interpolation\n", + "\n", + "i = range(0,5)\n", + "X = range(1,6)\n", + "Fx = [ 0, 7 ,26 ,63 ,124]\n", + "#x = Symbol('x')\n", + "A=zeros([5,7])\n", + "A[:,0]=i\n", + "A[:,1]=X\n", + "A[:,2]=Fx\n", + "\n", + "for i in range(4,8):\n", + " for j in range(1,9-i):\n", + " A[j-1,i-1] = ( A[j,i-2]-A[j-1,i-2] )/(X[j+i-4]-X[j-1]) \n", + " \n", + "\n", + "p = A[0,2]\n", + "p1_5 = [p,0,0,0,0,0,0,0] \n", + "x=1.5\n", + "for i in range(4,8):\n", + " p = p +A[0,i-1]* prod(array([x*xx for xx in ones([1,i-3])]) - array(X[0:i-3]))\n", + " p1_5[i-3] = p#horner(p,1.5)#\n", + "\n", + "print ' p0(1.5) = %f \\n p1(1.5) = %f \\n p2(1.5) = %f \\n p3(1.5) = %f \\n p4(1.5) = %f \\n'%(p1_5[0],p1_5[1],p1_5[2],p1_5[3],p1_5[4])\n", + "print 'The function value at x = 1.5 is : ',p1_5[4] " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 9_08 Pg No. 297" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "A = \n", + "[[ 1.00000000e+01 1.73600000e-01 1.68400000e-01 -1.04000000e-02\n", + " -4.80000000e-03 4.00000000e-04]\n", + " [ 2.00000000e+01 3.42000000e-01 1.58000000e-01 -1.52000000e-02\n", + " -4.40000000e-03 0.00000000e+00]\n", + " [ 3.00000000e+01 5.00000000e-01 1.42800000e-01 -1.96000000e-02\n", + " 0.00000000e+00 0.00000000e+00]\n", + " [ 4.00000000e+01 6.42800000e-01 1.23200000e-01 0.00000000e+00\n", + " 0.00000000e+00 0.00000000e+00]\n", + " [ 5.00000000e+01 7.66000000e-01 0.00000000e+00 0.00000000e+00\n", + " 0.00000000e+00 0.00000000e+00]]\n", + "\n", + " p1(s) = 0.3472 \n", + " p2(s) = 0.3472 \n", + " p3(s) = 0.3472 \n", + " p4(s) = 0.3472 \n", + "\n", + "Thus sin(25) = 0.3472 \n", + " \n" + ] + } + ], + "source": [ + "from numpy import diff, zeros, prod, array, ones\n", + "from scipy.misc import factorial\n", + "#Newton-Gregory forward difference formula\n", + "\n", + "X = [ 10, 20 ,30, 40, 50]\n", + "Fx = [ 0.1736, 0.3420 ,0.5000 ,0.6428, 0.7660]\n", + "#x = poly(0,'x'#\n", + "\n", + "A=zeros([5,6])\n", + "A[:,0]=X\n", + "A[:,1]=Fx\n", + "\n", + "\n", + "for i in range(3,7):\n", + " A[0:7-i,i-1] = diff(A[0:8-i,i-2])\n", + " \n", + "print 'A = \\n',A\n", + "\n", + "x0 = X[0]\n", + "h = X[1] - X[0] #\n", + "x1 = 25\n", + "s = (x1 - x0)/h #\n", + "p = [Fx[0]] \n", + "\n", + "for j in range(0,4):\n", + " p.append(p[j] + prod( array([s*xx for xx in ones([1,j+1])])-array([range(0,j+1)]) )*A[0,j+1]/factorial(j+1))\n", + "\n", + "print '\\n p1(s) = %.4G \\n p2(s) = %.4G \\n p3(s) = %.4G \\n p4(s) = %.4G \\n'%(p[1],p[2],p[3],p[4]) \n", + "print 'Thus sin(%d) = %.4G \\n '%(x1,p[4]) " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 9_09 Pg No. 299" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "A=\n", + "[[ 1.00000000e+01 1.73600000e-01 0.00000000e+00 0.00000000e+00\n", + " 0.00000000e+00 0.00000000e+00]\n", + " [ 2.00000000e+01 3.42000000e-01 1.68400000e-01 0.00000000e+00\n", + " 0.00000000e+00 0.00000000e+00]\n", + " [ 3.00000000e+01 5.00000000e-01 1.58000000e-01 -1.04000000e-02\n", + " 0.00000000e+00 0.00000000e+00]\n", + " [ 4.00000000e+01 6.42800000e-01 1.42800000e-01 -1.52000000e-02\n", + " -4.80000000e-03 0.00000000e+00]\n", + " [ 5.00000000e+01 7.66000000e-01 1.23200000e-01 -1.96000000e-02\n", + " -4.40000000e-03 4.00000000e-04]]\n", + "\n", + " s = -3\n", + "\n", + "\n", + " p1(s) = 0.3964 \n", + " p2(s) = 0.2788 \n", + " p3(s) = 0.3228 \n", + " p4(s) = 0.3288 \n", + "\n", + " Thus sin(25) = 0.3288 \n", + " \n" + ] + } + ], + "source": [ + "from numpy import diff, prod, array, ones, zeros\n", + "from scipy.misc import factorial\n", + "#Newton Backward difference formula\n", + "\n", + "X = [ 10, 20, 30 ,40, 50]\n", + "Fx = [ 0.1736, 0.3420, 0.5000, 0.6428, 0.7660]\n", + "#x = poly(0,'x'#\n", + "#A = [X' Fx']#\n", + "A=zeros([5,6])\n", + "A[:,0]=X\n", + "A[:,1]=Fx\n", + "\n", + "\n", + "for i in range(2,6):\n", + " A[i-1:5,i] = diff(A[i-2:5,i-1])\n", + "\n", + "print 'A=\\n',A\n", + "\n", + "xn = X[4]\n", + "h = 10 #\n", + "xuk = 25#\n", + "s = (xuk - xn)/h #\n", + "print '\\n s =',s\n", + "p = [Fx[4]]\n", + "\n", + "for j in range(1,5):\n", + " p.append(p[j-1] + prod(array([s*xx for xx in ones([1,j])]-array([range(0,j)])))*A[4,j+1]/factorial(j) )\n", + " \n", + "print '\\n\\n p1(s) = %.4G \\n p2(s) = %.4G \\n p3(s) = %.4G \\n p4(s) = %.4G \\n'%(p[1],p[2],p[3],p[4]) \n", + "print ' Thus sin(%d) = %.4G \\n '%(xuk,p[4]) " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 9_10 Pg No. 301" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "case 1:\n", + "The piecewise polynomials are continuous and f(x) is a linear spline\n", + "case 2:\n", + "The 1th derivative of polynomial is not continuours\n", + "case 3\n", + "The polynomial is continuous and its derivatives from 1 to 1 are continuous, f(x) is a 2th degree polynomial\n" + ] + } + ], + "source": [ + "from sympy import symbols, degree\n", + "from sympy.polys.polyfuncs import horner\n", + "\n", + "x = symbols('x')\n", + "def isitspline(f1,f2,f3,x0,x1,x2,x3):\n", + " n1 = degree(f1)\n", + " n2 = degree(f2)\n", + " n3 = degree(f3)\n", + " n = max(n1,n2,n3)\n", + " f1_x1 = f1.subs(x,x1)\n", + " f2_x1 = f2.subs(x,x1)\n", + " f2_x2 = f2.subs(x,x2)\n", + " f3_x2 = f3.subs(x,x2)\n", + " if n ==1 and f1_x1 == f2_x1 and f2_x2 == f3_x2:\n", + " print 'The piecewise polynomials are continuous and f(x) is a linear spline'\n", + " elif f1_x1 == f2_x1 and f2_x2 == f3_x2:\n", + " for i in range(1,n):\n", + " df1 = f1.diff()\n", + " df2 = f2.diff()\n", + " df3 = f3.diff()\n", + " df1_x1 = df1.subs(x,x1)\n", + " df2_x1 = df2.subs(x,x1)\n", + " df2_x2 = df2.subs(x,x2)\n", + " df3_x2 = df3.subs(x,x2)\n", + " f1 = df1\n", + " f2 = df2\n", + " f3 = df3\n", + " if df1_x1 != df2_x1 or df2_x2 != df3_x2:\n", + " print 'The %dth derivative of polynomial is not continuours'%i\n", + " break;\n", + " \n", + " \n", + " if i == n-1 and df1_x1 == df2_x1 and df2_x2 == df3_x2:\n", + " print 'The polynomial is continuous and its derivatives from 1 to %i are continuous, f(x) is a %ith degree polynomial'%(i,i+1)\n", + " \n", + " else:\n", + " print 'The polynomial is not continuous'\n", + " \n", + " \n", + "n = 4 \n", + "x0 = -1 \n", + "x1 = 0\n", + "x2 = 1\n", + "x3 = 2\n", + "f1 = x+1 ;\n", + "f2 = 2*x + 1 ;\n", + "f3 = 4 - x ;\n", + "print 'case 1:'\n", + "isitspline(f1,f2,f3,x0,x1,x2,x3)\n", + "n = 4\n", + "x0 = 0 \n", + "x1= 1 \n", + "x2 = 2 \n", + "x3 = 3\n", + "f1 = x**2 + 1 ;\n", + "f2 = 2*x**2 ;\n", + "f3 = 5*x - 2 ;\n", + "print 'case 2:'\n", + "isitspline(f1,f2,f3,x0,x1,x2,x3)\n", + "n = 4\n", + "x0 = 0\n", + "x1 = 1\n", + "x2 = 2\n", + "x3 = 3\n", + "f1 = x\n", + "f2 = x**2 - x + 1\n", + "f3 = 3*x - 3\n", + "print 'case 3'\n", + "isitspline(f1,f2,f3,x0,x1,x2,x3)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 9_11 Pg No. 306" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "h= [5 7]\n", + "a1 = -0.0142857142857\n", + "s(x) = 0.497619047619048*x - 0.0119047619047619*(1.0*x - 4)**3 + 0.00952380952380905\n", + "s(7) : 3.17142857142857\n" + ] + } + ], + "source": [ + "from numpy import array,diff, zeros, ones\n", + "from sympy import symbols\n", + "from __future__ import division\n", + "X = [ 4, 9, 16]\n", + "Fx = [ 2, 3, 4]\n", + "n = len(X)\n", + "h = diff(X)\n", + "print 'h=',h\n", + "x = symbols('x')\n", + "#A(1) = 0;\n", + "#A(n) = 0;\n", + "A=zeros(n)\n", + "#Since we do not know only a1 = A(2) we just have one equation which can be solved directly without solving tridiagonal matrix\n", + "A[1] = 6*( ( Fx[2] - Fx[1] )/h[1] - ( Fx[1] - Fx[0] )/h[0] )/( 2*( h[0] + h[1] ) )\n", + "print 'a1 = ',A[1]\n", + "xuk = 7;\n", + "for i in range(1,n):\n", + " if xuk <= X[i]:\n", + " break;\n", + " \n", + "\n", + "#u = x*ones([1,2]) - X[i-1:i+1]\n", + "u = array([x*xx for xx in ones([1,2])]) - array(X[i-1:i+1])\n", + "s = ( A[1]*( u[0][i-1]**3 - ( h[i-1]**2 )*u[0][i-1])/6*h[i-1] ) + ( Fx[i]*u[0][i-1]- Fx[i-1]*u[0][i] )/h[i-1]\n", + "print 's(x) =',s\n", + "s_7 = s.subs(x,xuk);\n", + "print 's(7) :',s_7" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example No. 9_12 Pg No. 313" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "h= [1 1 1]\n", + "[[4 1]\n", + " [1 4]]\n", + "D= [ 0.5004 0.1998]\n", + "A= [ 0. 0.12012 0.01992 0. ]\n", + "s(x) = -0.0666*x - 0.02002*(1.0*x - 3)**2 + 0.00332*(1.0*x - 2)**3 + 0.44648\n", + "s(2.5): 0.275390000000000\n" + ] + } + ], + "source": [ + "from numpy import diff, diag, transpose, zeros,array, ones\n", + "from sympy import symbols\n", + "from numpy.linalg import solve\n", + "#Cubic Spline Interpolation\n", + "\n", + "X = range(1,5)\n", + "Fx = [ 0.5, 0.3333, 0.25, 0.2]\n", + "n = len(X)\n", + "h = diff(X)\n", + "print 'h=',h\n", + "x = symbols('x')\n", + "A=zeros(n)\n", + "#Forming Tridiagonal Matrix\n", + "#take make diagonal below main diagonal be 1 , main diagonal is 2 and diagonal above main diagonal is 3\n", + "diag1 = h[1:n-2]\n", + "diag2 = 2*(h[0:n-2]+h[1:n-1])\n", + "diag3 = h[1:n-2]\n", + "TridiagMat = diag(diag1,-1)+diag(diag2)+diag(diag3,1)\n", + "print TridiagMat\n", + "D = diff(Fx)#\n", + "D = 6*diff(D/h)\n", + "print 'D=',D\n", + "A[1:n-1] = solve(array(TridiagMat),array(D))\n", + "print 'A=',A\n", + "xuk = 2.5;\n", + "for i in range(1,n):\n", + " if xuk <= X[i]:\n", + " break;\n", + " \n", + "\n", + "u = array([x*xx for xx in ones([1,2])]) - array(X[i-1:i+1])\n", + "s = ( A[i-1]*( h[i]**2*u[0][1] - u[0][1]**2 )/( 6*h[i] ) ) + ( A[i]*( u[0][0]**3 - ( h[i-1]**2 )*u[0][0])/6*h[i-1] ) + ( Fx[i]*u[0][0]- Fx[i-1]*u[0][1] )/h[i-1];\n", + "print 's(x) = ',s\n", + "s2_5 = s.subs(x,xuk)\n", + "print 's(2.5):',s2_5" + ] + } + ], + "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 +} diff --git a/Numerical_Methods_by_E._Balaguruswamy/screenshots/greatest-precision-4.png b/Numerical_Methods_by_E._Balaguruswamy/screenshots/greatest-precision-4.png Binary files differnew file mode 100644 index 00000000..27825c68 --- /dev/null +++ b/Numerical_Methods_by_E._Balaguruswamy/screenshots/greatest-precision-4.png diff --git a/Numerical_Methods_by_E._Balaguruswamy/screenshots/rounding-off-4.png b/Numerical_Methods_by_E._Balaguruswamy/screenshots/rounding-off-4.png Binary files differnew file mode 100644 index 00000000..ff48734e --- /dev/null +++ b/Numerical_Methods_by_E._Balaguruswamy/screenshots/rounding-off-4.png diff --git a/Numerical_Methods_by_E._Balaguruswamy/screenshots/truncation-error-4.png b/Numerical_Methods_by_E._Balaguruswamy/screenshots/truncation-error-4.png Binary files differnew file mode 100644 index 00000000..fbdf6b2c --- /dev/null +++ b/Numerical_Methods_by_E._Balaguruswamy/screenshots/truncation-error-4.png |
