diff options
author | Trupti Kini | 2017-01-12 23:33:05 +0600 |
---|---|---|
committer | Trupti Kini | 2017-01-12 23:33:05 +0600 |
commit | e2ec570a69371634ec7a6b27f9891e732c3428da (patch) | |
tree | 9a1604378caf041459fa446926a3ba5cd0f31eed /Numerical_Methods_by_E._Balaguruswamy/chapter9.ipynb | |
parent | 882a643ff3426d410a375df9c18ae2f90e045728 (diff) | |
download | Python-Textbook-Companions-e2ec570a69371634ec7a6b27f9891e732c3428da.tar.gz Python-Textbook-Companions-e2ec570a69371634ec7a6b27f9891e732c3428da.tar.bz2 Python-Textbook-Companions-e2ec570a69371634ec7a6b27f9891e732c3428da.zip |
Added(A)/Deleted(D) following books
A Basic_mechanical_engineering_by_Basant_Agrawal_,_C.M_Agrawal/README.txt
A Numerical_Methods_by_E._Balaguruswamy/chapter10.ipynb
A Numerical_Methods_by_E._Balaguruswamy/chapter11.ipynb
A Numerical_Methods_by_E._Balaguruswamy/chapter12.ipynb
A Numerical_Methods_by_E._Balaguruswamy/chapter13.ipynb
A Numerical_Methods_by_E._Balaguruswamy/chapter14.ipynb
A Numerical_Methods_by_E._Balaguruswamy/chapter15.ipynb
A Numerical_Methods_by_E._Balaguruswamy/chapter3.ipynb
A Numerical_Methods_by_E._Balaguruswamy/chapter4.ipynb
A Numerical_Methods_by_E._Balaguruswamy/chapter6.ipynb
A Numerical_Methods_by_E._Balaguruswamy/chapter7.ipynb
A Numerical_Methods_by_E._Balaguruswamy/chapter8.ipynb
A Numerical_Methods_by_E._Balaguruswamy/chapter9.ipynb
A Numerical_Methods_by_E._Balaguruswamy/screenshots/greatest-precision-4.png
A Numerical_Methods_by_E._Balaguruswamy/screenshots/rounding-off-4.png
A Numerical_Methods_by_E._Balaguruswamy/screenshots/truncation-error-4.png
A Physics_BSc(Paper_4)_by_Sanjeeva_Rao,_Bhikshmaiah,_Ramakrishna_Reddy,_Ananta_Ramaiah/C_4JhcI7F.ipynb
A Physics_BSc(Paper_4)_by_Sanjeeva_Rao,_Bhikshmaiah,_Ramakrishna_Reddy,_Ananta_Ramaiah/C_5CfOfQx.ipynb
A Physics_BSc(Paper_4)_by_Sanjeeva_Rao,_Bhikshmaiah,_Ramakrishna_Reddy,_Ananta_Ramaiah/C_B30VPml.ipynb
A Physics_BSc(Paper_4)_by_Sanjeeva_Rao,_Bhikshmaiah,_Ramakrishna_Reddy,_Ananta_Ramaiah/C_BE8QjoS.ipynb
A Physics_BSc(Paper_4)_by_Sanjeeva_Rao,_Bhikshmaiah,_Ramakrishna_Reddy,_Ananta_Ramaiah/C_Gb015bo.ipynb
A Physics_BSc(Paper_4)_by_Sanjeeva_Rao,_Bhikshmaiah,_Ramakrishna_Reddy,_Ananta_Ramaiah/C_KQ8ycMr.ipynb
A Physics_BSc(Paper_4)_by_Sanjeeva_Rao,_Bhikshmaiah,_Ramakrishna_Reddy,_Ananta_Ramaiah/C_LIZWeY4.ipynb
A Physics_BSc(Paper_4)_by_Sanjeeva_Rao,_Bhikshmaiah,_Ramakrishna_Reddy,_Ananta_Ramaiah/C_Opo6g3L.ipynb
A Physics_BSc(Paper_4)_by_Sanjeeva_Rao,_Bhikshmaiah,_Ramakrishna_Reddy,_Ananta_Ramaiah/C_Pq40WOu.ipynb
A Physics_BSc(Paper_4)_by_Sanjeeva_Rao,_Bhikshmaiah,_Ramakrishna_Reddy,_Ananta_Ramaiah/C_QLXi0uM.ipynb
A Physics_BSc(Paper_4)_by_Sanjeeva_Rao,_Bhikshmaiah,_Ramakrishna_Reddy,_Ananta_Ramaiah/C_TdOIeIQ.ipynb
A Physics_BSc(Paper_4)_by_Sanjeeva_Rao,_Bhikshmaiah,_Ramakrishna_Reddy,_Ananta_Ramaiah/C_UiM06tF.ipynb
A Physics_BSc(Paper_4)_by_Sanjeeva_Rao,_Bhikshmaiah,_Ramakrishna_Reddy,_Ananta_Ramaiah/C_VQabJzR.ipynb
A Physics_BSc(Paper_4)_by_Sanjeeva_Rao,_Bhikshmaiah,_Ramakrishna_Reddy,_Ananta_Ramaiah/C_XKPPUxj.ipynb
A Physics_BSc(Paper_4)_by_Sanjeeva_Rao,_Bhikshmaiah,_Ramakrishna_Reddy,_Ananta_Ramaiah/C_dfFlLnm.ipynb
A Physics_BSc(Paper_4)_by_Sanjeeva_Rao,_Bhikshmaiah,_Ramakrishna_Reddy,_Ananta_Ramaiah/C_fegIkl6.ipynb
A Physics_BSc(Paper_4)_by_Sanjeeva_Rao,_Bhikshmaiah,_Ramakrishna_Reddy,_Ananta_Ramaiah/C_g1GxlUN.ipynb
A Physics_BSc(Paper_4)_by_Sanjeeva_Rao,_Bhikshmaiah,_Ramakrishna_Reddy,_Ananta_Ramaiah/C_gGvuusH.ipynb
A Physics_BSc(Paper_4)_by_Sanjeeva_Rao,_Bhikshmaiah,_Ramakrishna_Reddy,_Ananta_Ramaiah/C_gZeNd4b.ipynb
A Physics_BSc(Paper_4)_by_Sanjeeva_Rao,_Bhikshmaiah,_Ramakrishna_Reddy,_Ananta_Ramaiah/C_khAt4Y6.ipynb
A Physics_BSc(Paper_4)_by_Sanjeeva_Rao,_Bhikshmaiah,_Ramakrishna_Reddy,_Ananta_Ramaiah/C_oJQM5Mb.ipynb
A Physics_BSc(Paper_4)_by_Sanjeeva_Rao,_Bhikshmaiah,_Ramakrishna_Reddy,_Ananta_Ramaiah/C_p22tFeA.ipynb
A Physics_BSc(Paper_4)_by_Sanjeeva_Rao,_Bhikshmaiah,_Ramakrishna_Reddy,_Ananta_Ramaiah/C_qNHcrb8.ipynb
A Physics_BSc(Paper_4)_by_Sanjeeva_Rao,_Bhikshmaiah,_Ramakrishna_Reddy,_Ananta_Ramaiah/C_qeUxd0G.ipynb
A Physics_BSc(Paper_4)_by_Sanjeeva_Rao,_Bhikshmaiah,_Ramakrishna_Reddy,_Ananta_Ramaiah/C_r0nfWNs.ipynb
A Physics_BSc(Paper_4)_by_Sanjeeva_Rao,_Bhikshmaiah,_Ramakrishna_Reddy,_Ananta_Ramaiah/C_uWQUwaW.ipynb
A Physics_BSc(Paper_4)_by_Sanjeeva_Rao,_Bhikshmaiah,_Ramakrishna_Reddy,_Ananta_Ramaiah/C_vkMAnbH.ipynb
A Physics_BSc(Paper_4)_by_Sanjeeva_Rao,_Bhikshmaiah,_Ramakrishna_Reddy,_Ananta_Ramaiah/C_wP2wGMS.ipynb
A Physics_BSc(Paper_4)_by_Sanjeeva_Rao,_Bhikshmaiah,_Ramakrishna_Reddy,_Ananta_Ramaiah/C_x09tDSO.ipynb
A Physics_BSc(Paper_4)_by_Sanjeeva_Rao,_Bhikshmaiah,_Ramakrishna_Reddy,_Ananta_Ramaiah/C_yCW2orc.ipynb
Diffstat (limited to 'Numerical_Methods_by_E._Balaguruswamy/chapter9.ipynb')
-rw-r--r-- | Numerical_Methods_by_E._Balaguruswamy/chapter9.ipynb | 798 |
1 files changed, 798 insertions, 0 deletions
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 +} |