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