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