summaryrefslogtreecommitdiff
path: root/Numerical_Methods_by_E._Balaguruswamy/chapter9.ipynb
diff options
context:
space:
mode:
authorTrupti Kini2017-01-12 23:33:05 +0600
committerTrupti Kini2017-01-12 23:33:05 +0600
commite2ec570a69371634ec7a6b27f9891e732c3428da (patch)
tree9a1604378caf041459fa446926a3ba5cd0f31eed /Numerical_Methods_by_E._Balaguruswamy/chapter9.ipynb
parent882a643ff3426d410a375df9c18ae2f90e045728 (diff)
downloadPython-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.ipynb798
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
+}