Chapter 9 - Curve fitting interpolation

Example No. 9_01 Pg No.277

In [1]:
from numpy.linalg import solve
from numpy import mat
from sympy import Symbol
print 'solving linear equations \n a0 + 100a1 = 3/7 \n a0 + 101a1 = -4/7 \n we get'#
C = mat([[ 1, 100],[1 ,101] ])
p = mat([[ 3/7],[-4/7] ])
a = solve(C,p )
print ' a0 = %.f \n a1 = %.f'%(a[0],a[1])

x = Symbol('x')
def horner(a,x):
    px = a[0] + a[1]*x
    return px
p100 = horner(a,100.0)
p101 = horner(a,101.0)
print '\n p(100) = %.f \n p(101) = %.f\n'%(p100,p101)

solving linear equations 
 a0 + 100a1 = 3/7 
 a0 + 101a1 = -4/7 
 we get
 a0 = 100 
 a1 = -1

 p(100) = 0 
 p(101) = -1

Example No. 9_02 Page No. 278

In [2]:
from numpy.linalg import solve
from numpy import mat
from sympy import Symbol

C = mat([[1, 100-100],[1, 101-100]])
p = mat([[ 3.0/7],[-4/7] ])
a = solve(C,p)
print '\n a0 = %f \n a1 = %f \n'%(a[0],a[1])

x = Symbol('x')
def horner(a,x):
    px = a[0]+ a[1]*(x - 100) 
    return px
p100 = horner(a,100)
p101 = horner(a,101)
print '\n p(100) = %f \n p(101) = %f\n'%(p100,p101)

 a0 = 0.428571 
 a1 = -1.428571 


 p(100) = 0.428571 
 p(101) = -1.000000

Example No. 9_03 Page No. 280

In [3]:
x = range(0,6)
f = [0,1, 1.4142, 1.7321, 2, 2.2361]
n = 2.5
for i in range(1,6):
    if n <= x[(i)]:
        break;
    

print '%.1f lies between points %d and %d'%(n,x[(i-1)],x[(i)])
f2_5 = f[(i-1)] + ( n - x[(i-1)] )*( f[(i)] - f[(i-1)] )/( x[(i)] - x[(i-1)] )
err1 = 1.5811 - f2_5
print 'f(2.5) = ',f2_5
print 'error1 = ',err1
print 'The correct answer is 1.5811.The difference between results is due to use of a linear model to a nonlinear use'
print 'repeating the procedure using x1 = 2 and x2 = 4'
x1 = 2
x2 = 4
f2_5 = f[(x1)] + ( 2.5 - x1 )*( f[(x2)] - f[(x1)] )/( x2 - x1 )
err2 = 1.5811 - f2_5
print 'error2 = ',err2
print 'f(2.5) = ',f2_5
print 'NOTE- The increase in error due to the increase in the interval between the interpolating data'

2.5 lies between points 2 and 3
f(2.5) =  1.57315
error1 =  0.00795
The correct answer is 1.5811.The difference between results is due to use of a linear model to a nonlinear use
repeating the procedure using x1 = 2 and x2 = 4
error2 =  0.02045
f(2.5) =  1.56065
NOTE- The increase in error due to the increase in the interval between the interpolating data

Example No. 9_04 Pg No. 282

In [3]:
from sympy import Symbol
from math import sqrt
#Lagrange Interpolation- Second order

X = [0, 1, 2 ,3 ,4 ,5]
Fx = [0, 1, 1.4142 ,1.7321 ,2, 2.2361]
X = X[2:5]
Fx = Fx[2:5]
x0 = 2.5 
x = Symbol('x')
p = 0
L=[0]
def horner(X,p,Fx,x,m):
    for i in range(1,3):
        L.append(1)
        for j in range(1,3):
            if j == i:
                continue #
            else:
                L[(i)] = L[(i)]*( x - X[(j)] )/( X[(i)] - X[(j)] ) 
        
        p = p + Fx[(i)]*L[(i)] 
    return [L[m],p]

x=2.5
L0 = horner(X,p,Fx,x,1)[0]
L1 = horner(X,p,Fx,x,2)[0]
L2 = horner(X,p,Fx,x,3)[0]
p2_5 = horner(X,p,Fx,x,3)[1]
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)

err = sqrt(2.5) - p2_5 #
print 'The error is %f'%(err)#    

For x = 2.5 we have,
 L0(2.5) = 1.500000 
 L1(2.5) = 0.250000 
 L2(2.5) = 1.000000 
 p(2.5) = 8.893756 

The error is -7.312617

Example No. 9_05 Pg No. 283

In [2]:
from sympy import Symbol
#Lagrange Interpolation

i = [ 0, 1, 2, 3 ]
X = [ 0 ,1 ,2 ,3 ]
Fx = [ 0 ,1.7183 ,6.3891, 19.0855 ]
x = Symbol('x')
n = 3 #order of lagrange polynomial 
p = 0
L=[]
for i in range(0,n+1):
    L.append(1)
    for j in range(0,n+1):
        if j == i:
            continue 
        else:
            L[i] = L[i]*( x - X[j] )/( X[i] - X[j] ) 
        
    
    p = p + Fx[i]*L[i]

print "The Lagrange basis polynomials are:"
for i in range(0,4):
        print str(L[i])


print "The interpolation polynomial is "
print str(p)

print 'The interpolation value at x = 1.5 is :', 

p1_5 = p.subs(x,1.5)
e1_5 = p1_5 + 1 #
print e1_5,'\ne**1.5 = %f'%(p1_5)#        

The Lagrange basis polynomials are:
(-x + 1)*(x - 3)*(x - 2)/6
x*(x - 3)*(x - 2)/2
-x*(x - 3)*(x - 1)/2
x*(x - 2)*(x - 1)/6
The interpolation polynomial is 
0.85915*x*(x - 3)*(x - 2) - 3.19455*x*(x - 3)*(x - 1) + 3.18091666666667*x*(x - 2)*(x - 1)
The interpolation value at x = 1.5 is : 4.36756875000000 
e**1.5 = 3.367569

Example No. 9_06 Pg No. 288

In [6]:
from numpy import zeros, prod, ones, array
from sympy.abc import x
i = [ 0, 1 ,2 ,3]
X = range(1,5)
Fx = [ 0, 0.3010, 0.4771, 0.6021] 
X = range(1,4)
Fx = Fx[0:3]
#x = poly(0,'x');
#A = Fx'
A=zeros([3,3])
A[:,0]=Fx
for i in range(2,4):
    for j in range(1,4-i+1):
        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]) ;

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])
p = A[0,0]

for i in range(2,4):
    p = p +A[0,i-1]* prod(array([x*xx for xx in ones([1,i-1])]) - array(X[0:i-1]))
print 'Polynomial is : ',str(p)
x=2.5
p=A[0,0]
for i in range(2,4):
    p = p +A[0,i-1]* prod(array([x*xx for xx in ones([1,i-1])]) - array(X[0:i-1]))
print 'p(2.5) = %.4G \n'%p

The coefficients of the polynomial are,
 a0 = 0 
 a1 = 0.301 
 a2 = -0.06245 

Polynomial is :  0.301*x - 0.06245*(1.0*x - 2)*(1.0*x - 1) - 0.301
p(2.5) = 0.4047

Example No. 9_07 Pg No. 291

In [1]:
from numpy import zeros, prod, ones, array
#Newton Divided Difference Interpolation

i = range(0,5)
X = range(1,6)
Fx = [ 0, 7 ,26 ,63 ,124]
#x = Symbol('x')
A=zeros([5,7])
A[:,0]=i
A[:,1]=X
A[:,2]=Fx

for i in range(4,8):
     for j in range(1,9-i):
        A[j-1,i-1] = ( A[j,i-2]-A[j-1,i-2] )/(X[j+i-4]-X[j-1]) 
    

p = A[0,2]
p1_5 = [p,0,0,0,0,0,0,0] 
x=1.5
for i in range(4,8):
    p = p +A[0,i-1]* prod(array([x*xx for xx in ones([1,i-3])]) - array(X[0:i-3]))
    p1_5[i-3] = p#horner(p,1.5)#

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])
print 'The function value at x = 1.5 is : ',p1_5[4]    

 p0(1.5) = 0.000000 
 p1(1.5) = 3.500000 
 p2(1.5) = 2.000000 
 p3(1.5) = 2.375000 
 p4(1.5) = 2.375000 

The function value at x = 1.5 is :  2.375

Example No. 9_08 Pg No. 297

In [5]:
from numpy import diff, zeros, prod, array, ones
from scipy.misc import factorial
#Newton-Gregory forward difference formula

X = [ 10, 20 ,30, 40, 50]
Fx = [ 0.1736, 0.3420 ,0.5000 ,0.6428, 0.7660]
#x = poly(0,'x'#

A=zeros([5,6])
A[:,0]=X
A[:,1]=Fx


for i in range(3,7):
     A[0:7-i,i-1] = diff(A[0:8-i,i-2])
        
print 'A = \n',A

x0 = X[0]
h = X[1] - X[0] #
x1 = 25
s = (x1 - x0)/h #
p = [Fx[0]] 

for j in range(0,4):
    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))

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]) 
print 'Thus sin(%d) = %.4G \n '%(x1,p[4]) 

A = 
[[  1.00000000e+01   1.73600000e-01   1.68400000e-01  -1.04000000e-02
   -4.80000000e-03   4.00000000e-04]
 [  2.00000000e+01   3.42000000e-01   1.58000000e-01  -1.52000000e-02
   -4.40000000e-03   0.00000000e+00]
 [  3.00000000e+01   5.00000000e-01   1.42800000e-01  -1.96000000e-02
    0.00000000e+00   0.00000000e+00]
 [  4.00000000e+01   6.42800000e-01   1.23200000e-01   0.00000000e+00
    0.00000000e+00   0.00000000e+00]
 [  5.00000000e+01   7.66000000e-01   0.00000000e+00   0.00000000e+00
    0.00000000e+00   0.00000000e+00]]

 p1(s) = 0.3472 
 p2(s) = 0.3472 
 p3(s) = 0.3472 
 p4(s) = 0.3472 

Thus sin(25) = 0.3472

Example No. 9_09 Pg No. 299

In [2]:
from numpy import diff, prod, array, ones, zeros
from scipy.misc import factorial
#Newton Backward difference formula

X = [ 10, 20, 30 ,40, 50]
Fx = [ 0.1736, 0.3420, 0.5000, 0.6428, 0.7660]
#x = poly(0,'x'#
#A = [X' Fx']#
A=zeros([5,6])
A[:,0]=X
A[:,1]=Fx


for i in range(2,6):
     A[i-1:5,i] = diff(A[i-2:5,i-1])

print 'A=\n',A

xn = X[4]
h = 10 #
xuk = 25#
s = (xuk - xn)/h #
print '\n s =',s
p = [Fx[4]]

for j in range(1,5):
    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) )
    
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]) 
print ' Thus sin(%d) = %.4G \n '%(xuk,p[4])     

A=
[[  1.00000000e+01   1.73600000e-01   0.00000000e+00   0.00000000e+00
    0.00000000e+00   0.00000000e+00]
 [  2.00000000e+01   3.42000000e-01   1.68400000e-01   0.00000000e+00
    0.00000000e+00   0.00000000e+00]
 [  3.00000000e+01   5.00000000e-01   1.58000000e-01  -1.04000000e-02
    0.00000000e+00   0.00000000e+00]
 [  4.00000000e+01   6.42800000e-01   1.42800000e-01  -1.52000000e-02
   -4.80000000e-03   0.00000000e+00]
 [  5.00000000e+01   7.66000000e-01   1.23200000e-01  -1.96000000e-02
   -4.40000000e-03   4.00000000e-04]]

 s = -3


 p1(s) = 0.3964 
 p2(s) = 0.2788 
 p3(s) = 0.3228 
 p4(s) = 0.3288 

 Thus sin(25) = 0.3288

Example No. 9_10 Pg No. 301

In [1]:
from sympy import symbols, degree
from sympy.polys.polyfuncs import horner

x = symbols('x')
def isitspline(f1,f2,f3,x0,x1,x2,x3):
    n1 = degree(f1)
    n2 = degree(f2)
    n3 = degree(f3)
    n = max(n1,n2,n3)
    f1_x1 = f1.subs(x,x1)
    f2_x1 = f2.subs(x,x1)
    f2_x2 = f2.subs(x,x2)
    f3_x2 = f3.subs(x,x2)
    if n ==1 and f1_x1 == f2_x1 and f2_x2 == f3_x2:
        print 'The piecewise polynomials are continuous and f(x) is a linear spline'
    elif f1_x1 == f2_x1 and f2_x2 == f3_x2:
        for i in range(1,n):
            df1 = f1.diff()
            df2 = f2.diff()
            df3 = f3.diff()
            df1_x1 = df1.subs(x,x1)
            df2_x1 = df2.subs(x,x1)
            df2_x2 = df2.subs(x,x2)
            df3_x2 = df3.subs(x,x2)
            f1 = df1
            f2 = df2
            f3 = df3
            if df1_x1 != df2_x1 or df2_x2 != df3_x2:
                print 'The %dth derivative of polynomial is not continuours'%i
                break;
            
        
        if i == n-1 and df1_x1 == df2_x1 and df2_x2 == df3_x2:
            print 'The polynomial is continuous and its derivatives from 1 to %i are continuous, f(x) is a %ith degree polynomial'%(i,i+1)
        
    else:
        print 'The polynomial is not continuous'
    
    
n = 4 
x0 = -1 
x1 = 0
x2 = 1
x3 = 2
f1 = x+1 ;
f2 = 2*x + 1 ;
f3 = 4 - x ;
print 'case 1:'
isitspline(f1,f2,f3,x0,x1,x2,x3)
n = 4
x0 = 0 
x1= 1 
x2 = 2 
x3 = 3
f1 = x**2 + 1 ;
f2 = 2*x**2 ;
f3 = 5*x - 2 ;
print 'case 2:'
isitspline(f1,f2,f3,x0,x1,x2,x3)
n = 4
x0 = 0
x1 = 1
x2 = 2
x3 = 3
f1 = x
f2 = x**2 - x + 1
f3 = 3*x - 3
print 'case 3'
isitspline(f1,f2,f3,x0,x1,x2,x3)

case 1:
The piecewise polynomials are continuous and f(x) is a linear spline
case 2:
The 1th derivative of polynomial is not continuours
case 3
The polynomial is continuous and its derivatives from 1 to 1 are continuous, f(x) is a 2th degree polynomial

Example No. 9_11 Pg No. 306

In [2]:
from numpy import array,diff, zeros, ones
from sympy import symbols
from __future__ import division
X = [ 4, 9, 16]
Fx = [ 2, 3, 4]
n = len(X)
h = diff(X)
print 'h=',h
x = symbols('x')
#A(1) = 0;
#A(n) = 0;
A=zeros(n)
#Since we do not know only a1 = A(2) we just have one equation which can be solved directly without solving tridiagonal matrix
A[1] = 6*( ( Fx[2] - Fx[1] )/h[1] - ( Fx[1] - Fx[0] )/h[0] )/( 2*( h[0] + h[1] ) )
print 'a1 = ',A[1]
xuk = 7;
for i in range(1,n):
    if xuk <= X[i]:
        break;
    

#u = x*ones([1,2]) - X[i-1:i+1]
u = array([x*xx for xx in ones([1,2])]) - array(X[i-1:i+1])
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]
print 's(x) =',s
s_7 = s.subs(x,xuk);
print 's(7) :',s_7

h= [5 7]
a1 =  -0.0142857142857
s(x) = 0.497619047619048*x - 0.0119047619047619*(1.0*x - 4)**3 + 0.00952380952380905
s(7) : 3.17142857142857

Example No. 9_12 Pg No. 313

In [1]:
from numpy import diff, diag, transpose, zeros,array, ones
from sympy import symbols
from numpy.linalg import solve
#Cubic Spline Interpolation

X = range(1,5)
Fx = [ 0.5, 0.3333, 0.25, 0.2]
n = len(X)
h = diff(X)
print 'h=',h
x = symbols('x')
A=zeros(n)
#Forming Tridiagonal Matrix
#take make diagonal below main diagonal be 1 , main diagonal is 2 and diagonal above main diagonal is 3
diag1 = h[1:n-2]
diag2 = 2*(h[0:n-2]+h[1:n-1])
diag3 = h[1:n-2]
TridiagMat = diag(diag1,-1)+diag(diag2)+diag(diag3,1)
print TridiagMat
D = diff(Fx)#
D = 6*diff(D/h)
print 'D=',D
A[1:n-1] = solve(array(TridiagMat),array(D))
print 'A=',A
xuk = 2.5;
for i in range(1,n):
    if xuk <= X[i]:
        break;
    

u = array([x*xx for xx in ones([1,2])]) - array(X[i-1:i+1])
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];
print 's(x) = ',s
s2_5 = s.subs(x,xuk)
print 's(2.5):',s2_5

h= [1 1 1]
[[4 1]
 [1 4]]
D= [ 0.5004  0.1998]
A= [ 0.       0.12012  0.01992  0.     ]
s(x) =  -0.0666*x - 0.02002*(1.0*x - 3)**2 + 0.00332*(1.0*x - 2)**3 + 0.44648
s(2.5): 0.275390000000000