Chapter03:Interpolation

Ex3.4:pg-86

In [7]:
#example 3.4
#interpolation
#page 86
import math
from __future__ import division
x=[1, 3, 5, 7]
y=[24, 120, 336, 720]
d1=[0,0,0]
d2=[0,0,0]
d3=[0,0,0]
h=2     #interval between values of x
c=0
for i in range(0,3):
    d1[c]=y[i+1]-y[i]
    c=c+1
c=0
for i in range(0,2):
    d2[c]=d1[i+1]-d1[i]
    c=c+1
c=0
for i in range(0,1):
    d3[c]=d2[i+1]-d2[i]
    c=c+1
d=[0,d1[0],d2[0],d3[0]]
x0=8   #value at 8
pp=1
y_x=y[0]
p=(x0-1)/2
for i in range(1,4):
    pp=1
    for j in range(0,i):
        pp=pp*(p-(j))  
    y_x=y_x+(pp*d[i])/math.factorial(i)
print "value of function at %f is :%f" %(x0,y_x)

value of function at 8.000000 is :990.000000

Ex3.6:pg-87

In [26]:
#example 3.6
#interpolation
#page 87
import math
x=[15, 20, 25, 30, 35, 40]
y=[0.2588190, 0.3420201, 0.4226183, 0.5, 0.5735764, 0.6427876]
d1=[0,0,0,0,0]
d2=[0,0,0,0]
d3=[0,0,0]
d4=[0,0]
d5=[0]
h=5  #interval between values of x
c=0
for i in range(0,5):
    d1[c]=y[i+1]-y[i]
    c=c+1
c=0
for i in range(0,4):
    d2[c]=d1[i+1]-d1[i]
    c=c+1
c=0
for i in range(0,3):
    d3[c]=d2[i+1]-d2[i]
    c=c+1
c=0
for i in range(0,2):
    d4[c]=d3[i+1]-d3[i]
    c=c+1
c=0
for i in range(0,1):
    d5[c]=d4[i+1]-d4[i]
    c=c+1
c=0
d=[0,d1[0], d2[0], d3[0], d4[0], d5[0]]
x0=38    #value at 38 degree
pp=1
y_x=y[0]
p=(x0-x[0])/h
for i in range(1,6):
    pp=1
    for j in range(0,i):
        pp=pp*(p-(j))   
    y_x=y_x+((pp*d[i])/math.factorial(i));
print "value of function at %i is :%f" %(x0,y_x)

value of function at 38 is :0.615661

Ex3.7:pg-89

In [1]:
#example 3.7
#interpolation
#page 89
x=[0, 1, 2, 4]
y=[1, 3, 9, 81]
#equation is y(5)-4*y(4)+6*y(2)-4*y(2)+y(1)
y3=(y[3]+6*y[2]-4*y[1]+y[0])/4
print "the value of missing term of table is :%d" %(y3)

the value of missing term of table is :31

Ex3.8:pg-89

In [33]:
#example 3.8
#interpolation
#page 89
import math
x=[0.10, 0.15, 0.20, 0.25, 0.30]
y=[0.1003, 0.1511, 0.2027, 0.2553, 0.3093]
d1=[0,0,0,0,0]
d2=[0,0,0,0,0]
d3=[0,0,0,0,0]
d4=[0,0,0,0,0]
h=0.05 #interval between values of x
c=0
for i in range(0,4):
    d1[c]=y[i+1]-y[i]
    c=c+1
c=0
for i in range(0,3):
    d2[c]=d1[i+1]-d1[i]
    c=c+1
c=0
for i in range(0,2):
    d3[c]=d2[i+1]-d2[i]
    c=c+1
c=0
for i in range(0,4):
    d4[c]=d3[i+1]-d3[i]
    c=c+1
d=[0,d1[0], d2[0], d3[0], d4[0]]
x0=0.12   #value at 0.12;
pp=1
y_x=y[0]
p=(x0-x[0])/h
for i in range(1,5):
    pp=1;
    for j in range(0,i):
        pp=pp*(p-(j))    
    y_x=y_x+(pp*d[i])/math.factorial(i)
print "value of function at %f is :%0.4g\n \n" %(x0,y_x)
x0=0.26  #value at 0.26;
pp=1
y_x=y[0]
p=(x0-x[0])/h
for i in range(1,5):
    pp=1
    for j in range(0,i):
        pp=pp*(p-(j))    
    y_x=y_x+(pp*d[i])/math.factorial(i);
print "value of function at %f is :%0.4g\n \n" %(x0,y_x)
x0=0.40   #value at 0.40;
pp=1
y_x=y[0]
p=(x0-x[0])/h
for i in range(1,5):
    pp=1
    for j in range(0,i):
         pp=pp*(p-(j))    
    y_x=y_x+(pp*d[i])/math.factorial(i)
print "value of function at %f is :%0.4g\n \n" %(x0,y_x)
x0=0.50  #value at 0.50;
pp=1
y_x=y[0]
p=(x0-x[0])/h
for i in range(1,5):
    pp=1
    for j in range(0,i):
         pp=pp*(p-(j))    
    y_x=y_x+(pp*d[i])/math.factorial(i)
print "value of function at %f is :%0.5g\n \n" %(x0,y_x)

value of function at 0.120000 is :0.1205
 

value of function at 0.260000 is :0.266
 

value of function at 0.400000 is :0.4241
 

value of function at 0.500000 is :0.5543

Ex3.9:pg-93

In [35]:
#example 3.9
#Gauss' forward formula
#page 93
import math
x=[1.0, 1.05, 1.10, 1.15, 1.20, 1.25, 1.30];
y=[2.7183, 2.8577, 3.0042, 3.1582, 3.3201, 3.4903, 3.66693]
d1=[0,0,0,0,0,0]
d2=[0,0,0,0,0]
d3=[0,0,0,0]
d4=[0,0,0]
d5=[0,0]
d6=[0]
h=0.05  #interval between values of x
c=0
for i in range(0,6):
    d1[c]=y[i+1]-y[i]
    c=c+1
c=0
for i in range(0,5):
    d2[c]=d1[i+1]-d1[i]
    c=c+1
c=0
for i in range(0,4):
    d3[c]=d2[i+1]-d2[i]
    c=c+1
c=0
for i in range(0,3):
    d4[c]=d3[i+1]-d3[i]
    c=c+1
c=0
for i in range(0,2):
    d5[c]=d4[i+1]-d4[i]
    c=c+1
c=0
for i in range(0,1):
    d6[c]=d5[i+1]-d5[i]
    c=c+1
d=[0,d1[3], d2[2], d3[2], d4[1], d5[0], d6[0]]
x0=1.17 #value at 1.17;
pp=1
y_x=y[3]
p=(x0-x[3])/h
for i in range(1,6):
    pp=1;
    for j in range(0,i):
        pp=pp*(p-(j))    
    y_x=y_x+(pp*d[i])/math.factorial(i)
print "value of function at %f is :%0.4g\n \n" %(x0,y_x)

value of function at 1.170000 is :3.222

Ex3.10:pg-97

In [37]:
#practical interpolation
#example 3.10
#page 97
import math
x=[0.61, 0.62, 0.63, 0.64, 0.65, 0.66, 0.67]
y=[1.840431, 1.858928,1.877610, 1.896481, 1.915541, 1.934792, 1.954237]
d1=[0,0,0,0,0,0]
d2=[0,0,0,0,0]
d3=[0,0,0,0]
d4=[0,0,0]
h=0.01 #interval between values of x
c=0
for i in range(0,6):
    d1[c]=y[i+1]-y[i]
    c=c+1
c=0
for i in range(0,5):
    d2[c]=d1[i+1]-d1[i]
    c=c+1
c=0
for i in range(0,4):
    d3[c]=d2[i+1]-d2[i];
    c=c+1
c=0
for i in range(0,3):
    d4[c]=d3[i+1]-d3[i];
    c=c+1
d=[d1[0], d2[0], d3[0], d4[0]]
x0=0.644
p=(x0-x[3])/h;
y_x=y[3]
y_x=y_x+p*(d1[2]+d1[3])/2+p**2*(d2[1])/2  #stirling formula
print "the value at %f by stirling formula is : %f\n\n" %(x0,y_x)
y_x=y[3]
y_x=y_x+p*d1[3]+p*(p-1)*(d2[2]+d2[3])/2
print " the value at %f by bessels formula is : %f\n\n" %(x0,y_x)
y_x=y[3]
q=1-p
y_x=q*y[3]+q*(q**2-1)*d2[2]/2+p*y[4]+p*(q**2-1)*d2[4]/2
print "the value at %f by everrets formula is : %f\n\n" %(x0,y_x)

the value at 0.644000 by stirling formula is : 1.904082


 the value at 0.644000 by bessels formula is : 1.904059


the value at 0.644000 by everrets formula is : 1.904044

Ex3.11:pg-99

In [39]:
#practical interpolation
#example 3.11
#page 99
x=[0.61, 0.62, 0.63, 0.64, 0.65, 0.66, 0.67]
y=[1.840431, 1.858928, 1.877610, 1.896481, 1.915541, 1.934792, 1.954237]
d1=[0,0,0,0,0,0]
d2=[0,0,0,0,0]
d3=[0,0,0,0]
d4=[0,0,0]
h=0.01   #interval between values of x
c=0
for i in range(0,6):
    d1[c]=y[i+1]-y[i]
    c=c+1
c=0
for i in range(0,5):
    d2[c]=d1[i+1]-d1[i]
    c=c+1
c=0
for i in range(0,4):
    d3[c]=d2[i+1]-d2[i]
    c=c+1
c=0
for i in range(0,3):
    d4[c]=d3[i+1]-d3[i]
    c=c+1
d=[d1[0], d2[0], d3[0], d4[0]]
x0=0.638
p=(x0-x[3])/h
y_x=y[3]
y_x=y_x+p*(d1[2]+d1[3])/2+p**2*(d2[1])/2   #stirling formula
print "value at %f by stirling formula is : %f\n\n" %(x0,y_x)
y_x=y[2]
p=(x0-x[2])/h
y_x=y_x+p*d1[2]+p*(p-1)*(d2[1])/2
print "the value at %f by bessels formula is : %f\n\n" %(x0,y_x)

value at 0.638000 by stirling formula is : 1.892692


the value at 0.638000 by bessels formula is : 1.892692

Ex3.12:pg-99

In [41]:
#practical interpolation
#example 3.12
#page 99
x=[1.72, 1.73, 1.74, 1.75, 1.76, 1.77, 1.78]
y=[0.1790661479, 0.1772844100, 0.1755204006, 0.1737739435, 0.1720448638, 0.1703329888, 0.1686381473]
d1=[0,0,0,0,0,0]
d2=[0,0,0,0,0]
d3=[0,0,0,0]
d4=[0,0,0]
h=0.01  #interval between values of x
c=0
for i in range(0,6):
    d1[c]=y[i+1]-y[i]
    c=c+1
c=0
for i in range(0,5):
    d2[c]=d1[i+1]-d1[i]
    c=c+1
c=0
for i in range(0,4):
    d3[c]=d2[i+1]-d2[i]
    c=c+1
c=0
for i in range(0,3):
    d4[c]=d3[i+1]-d3[i]
    c=c+1
x0=1.7475
y_x=y[2]
p=(x0-x[2])/h
y_x=y_x+p*d1[2]+p*(p-1)*((d2[1]+d2[2])/2)/2
print "the value at %f by bessels formula is : %0.10f\n\n" %(x0,y_x)
y_x=y[3]
q=1-p
y_x=q*y[2]+q*(q**2-1)*d2[1]/6+p*y[3]+p*(p**2-1)*d2[1]/6
print "the value at %f by everrets formula is : %0.10f\n\n" %(x0,y_x)

the value at 1.747500 by bessels formula is : 0.1742089204


the value at 1.747500 by everrets formula is : 0.1742089122

Ex3.13:pg-104

In [43]:
#example 3.13
#lagrange's interpolation formula
#page 104
x=[300, 304, 305, 307]
y=[2.4771, 2.4829, 2.4843, 2.4871]
x0=301
log_301=(-3*-4*-6*2.4771)/(-4*-5*-7)+(-4*-6*2.4829)/(4*-1*-3)+(-3*-6*2.4843)/(5*-2)+(-3*-4*2.4871)/(7*3*2)
print "valie of log x at 301 is =%f" %(log_301)

 

valie of log x at 301 is =2.478597

Ex3.14:pg-105

In [44]:
#example 3.14
#lagrange's interpolation formula
#page 105
y=[4, 12, 19]
x=[1, 3, 4];
y_x=7
Y_X=(-5*-12)/(-8*-15)+(3*3*-12)/(8*-7)+(3*-5*4)/(15*7)
print "values is %f" %(Y_X)
 

values is 1.857143

Ex3.15:pg-105

In [2]:
#example 3.15
#lagrange's interpolation formula
#page 105
import math
x=[2, 2.5, 3.0]
y=[0.69315, 0.91629, 1.09861]
def l0(x):
    return (x-2.5)*(x-3.0)/(-0.5)*(-1.0)
def l1(x):
    return ((x-2.0)*(x-3.0))/((0.5)*(-0.5))
def l2(x):
    return ((x-2.0)*(x-2.5))/((1.0)*(0.5))
f_x=l0(2.7)*y[0]+l1(2.7)*y[1]+l2(2.7)*y[2];
print "the calculated value is %f:" %(f_x)
print "\n\n the error occured in the value is %0.9f" %(abs(f_x-log(2.7)))

the calculated value is 0.994116:


 the error occured in the value is 0.000864627

Ex3.16:pg-106

In [4]:
#example 3.16
#lagrange's interpolation formula
#page 106
import math
x=[0, math.pi/4,math.pi/2]
y=[0, 0.70711, 1.0];
x0=math.pi/6
sin_x0=0
for i in range(0,3):
    p=y[i]
    for j in range(0,3):
        if j!=i:
            p=p*((x0-x[j])/( x[i]-x[j]))
    sin_x0=sin_x0+p
print "sin_x0=%f" %(sin_x0)
 

sin_x0=0.517431

Ex3.18:pg-107

In [13]:
#error in lagrange's interpolation formula
#example 3.18
#page 107
import math
x=[2, 2.5, 3.0]
y=[0.69315, 0.91629, 1.09861]
def l0(x):
    return (x-2.5)*(x-3.0)/(-0.5)*(-1.0)
def l1(x):
    return ((x-2.0)*(x-3.0))/((0.5)*(-0.5))
def l2(x):
    return ((x-2.0)*(x-2.5))/((1.0)*(0.5))
f_x=l0(2.7)*y[0]+l1(2.7)*y[1]+l2(2.7)*y[2]
print "the calculated value is %f:" %(f_x)
err=math.fabs(f_x-math.log10(2.7))
def R_n(x):
    return (((x-2)*(x-2.5)*(x-3))/6)
est_err=abs(R_n(2.7)*(2/8))
if est_err<err:
    print "\n\n the error agrees with the actual error"

the calculated value is 0.994116:


 the error agrees with the actual error

Ex3.19:pg-107

In [14]:
#error in lagrenge's interpolation
#example 3.19
#page 107
import math
x=[0, math.pi/4 ,math.pi/2]
y=[0, 0.70711, 1.0]
def l0(x):
    return ((x-0)*(x-math.pi/2))/((math.pi/4)*(-1*math.pi/4))
def l1(x):
    return ((x-0)*(x-math.pi/4))/((math.pi/2)*(math.pi/4))
f_x=l0(math.pi/6)*y[1]+l1(math.pi/6)*y[2]
err=abs(f_x-math.sin(math.pi/6))
def f(x):
    return ((x-0)*(x-math.pi/4)*(x-math.pi/2))/6
if abs(f(math.pi/6))>err:
    print "\n\n the error agrees with the actual error"


 the error agrees with the actual error

Ex3.21:pg-110

In [35]:
#hermite's interpolation formula
#exammple 3.21
#page 110
from __future__ import division
import math
x=[2.0, 2.5, 3.0]
y=[0.69315, 0.91629, 1.09861]
y1=[0,0,0]
def f(x):
    return math.log(x)
h=0.0001
for i in range(0,3):
    y1[i]=(f(x[i]+h)-f(x[i]))/h
def l0(x):
    return (x-2.5)*(x-3.0)/(-0.5)*(-1.0)
def l1(x):
    return ((x-2.0)*(x-3.0))/((0.5)*(-0.5))
def l2(x):
    return ((x-2.0)*(x-2.5))/((1.0)*(0.5))
dl0=(l0(x[0]+h)-l0(x[0]))/h
dl1=(l1(x[1]+h)-l1(x[1]))/h
dl2=(l2(x[2]+h)-l2(x[2]))/h
x0=2.7
u0=(1-2*(x0-x[0])*dl0)*(l0(x0))**2
u1=(1-2*(x0-x[1])*dl1)*(l1(x0))**2
u2=(1-2*(x0-x[2])*dl2)*(l2(x0))**2
v0=(x0-x[0])*l0(x0)**2
v1=(x0-x[1])*l1(x0)**2
v2=(x0-x[2])*l2(x0)**2
H=u0*y[0]+u1*y[1]+u2*y[2]+v0*y1[0]+v1*y1[1]+v2*y1[2]
print "the approximate value of ln(%0.2f) is %f:" %(x0,H)

the approximate value of ln(2.70) is 0.993362:

Ex3.22:pg-114

In [37]:
#newton's general interpolation formula
#example 3.22
#page 114
x=[300, 304, 305, 307]
y=[2.4771, 2.4829, 2.4843, 2.4871]
d1=[0,0,0]
d2=[0,0]
for i in range(0,3):
    d1[i]=(y[i+1]-y[i])/(x[i+1]-x[i])
for i in range(0,2):
    d2[i]=(d1[i+1]-d1[i])/(x[i+2]-x[i])
x0=301
log301=y[0]+(x0-x[0])*d1[0]+(x0-x[1])*d2[0]
print "valure of log(%d) is :%0.4f" %(x0,log301)

valure of log(301) is :2.4786

Ex3.23:pg-114

In [18]:
#example 3.23
#newton's divided formula
#page 114
from sympy import simplify
x=[-1, 0, 3, 6, 7]
y=[3, -6, 39, 822, 1611]
d1=[0,0,0,0,0]
d2=[0,0,0,0,0]
d3=[0,0,0,0,0]
d4=[0,0,0,0,0]
X=0
for i in range(0,4):
    d1[i]=(y[i+1]-y[i])/(x[i+1]-x[i])
for i in range(0,3):
    d2[i]=(d1[i+1]-d1[i])/(x[i+2]-x[i])  
for i in range(0,2):
    d3[i]=(d2[i+1]-d2[i])/(x[i+3]-x[i])  
for i in range(0,1):
    d4[i]=(d3[i+1]-d3[i])/(x[i+4]-x[i]) 
X=Symbol('X')
f_x=simplify(y[0]+(X-x[0])*d1[0]+(X-x[1])*(X-x[0])*d2[0]+(X-x[0])*(X-x[1])*(X-x[2])*d3[0]+(X-x[0])*(X-x[1])*(X-x[2])*(X-x[3])*d4[0])
print f_x

1.0*X**4 - 3.0*X**3 + 5.0*X**2 - 6.0

Ex3.24:pg-116

In [2]:
#interpolation by iteration
#example 3.24
#page 116
import numpy
x=[300, 304, 305, 307]
y=[2.4771, 2.4829, 2.4843, 2.4871]
x0=301
d1=[0,0,0]
d2=[0,0]
d3=[0]
for i in range(0,3):
    a=y[i]
    b=x[i]-x0
    c=y[i+1]
    e=x[i+1]-x0
    d=matrix([[a,b],[c,e]])
    d11=numpy.linalg.det(d)
    d1[i]=d11/(x[i+1]-x[i])
for i in range(0,2):
    a=d1[i]
    b=x[i+1]-x0
    c=d1[i+1]
    e=x[i+2]-x0
    d=matrix([[a,b],[c,e]])
    d22=numpy.linalg.det(d)
    f=(x[i+2]-x[i+1])
    d2[i]=d22/f
for i in range(0,1):
    a=d2[i]
    b=x[i+2]-x0
    c=d2[i+1]
    e=x[i+3]-x0
    d=matrix([[a,b],[c,e]])
    d33=numpy.linalg.det(d)
    d3[i]=d33/(x[i+3]-x[i+2])
print "the value of log(%d) is : %f" %(x0,d3[0])

the value of log(301) is : 2.476900

Ex3.25:pg-118

In [3]:
#inverse intrpolation
#example 3.25
#page  118
from __future__ import division
x=[2, 3, 4, 5]
y=[8, 27, 64, 125]
d1=[0,0,0]
d2=[0,0]
d3=[0]
for i in range(0,3):
    d1[i]=y[i+1]-y[i]
for i in range(0,2):
    d2[i]=d1[i+1]-d1[i]
for i in range(0,1):
    d3[i]=d2[i+1]-d2[i]
yu=10   #square rooot of 10
y0=y[0]
d=[d1[0], d2[0] ,d3[0]]
u1=(yu-y0)/d1[0]
u2=((yu-y0-u1*(u1-1)*d2[0]/2)/d1[0])
u3=(yu-y0-u2*(u2-1)*d2[0]/2-u2*(u2-1)*(u2-2)*d3[0]/6)/d1[0]
u4=(yu-y0-u3*(u3-1)*d2[0]/2-u3*(u3-1)*(u3-2)*d3[0]/6)/d1[0]
u5=(yu-y0-u4*(u4-1)*d2[0]/2-u4*(u4-1)*(u4-2)*d3[0]/6)/d1[0]
print "%f \n %f \n %f \n %f \n %f \n " %(u1,u2,u3,u4,u5)
print "the approximate square root of %d is: %0.3f" %(yu,x[0]+u5)

0.105263 
 0.149876 
 0.153210 
 0.154107 
 0.154347 
 
the approximate square root of 10 is: 2.154

Ex3.26:pg-119

In [4]:
#double interpolation 
#example 3.26
#page 119
y=[0, 1, 2, 3, 4]
z=[0,0,0,0,0]
x=[[0, 1, 4, 9, 16],[2, 3, 6, 11, 18],[6, 7, 10, 15, 22],[12, 13, 16, 21, 28],[18, 19, 22, 27, 34]]
print "X="
print x
#for x=2.5
for i in range(0,5):
    z[i]=(x[i][2]+x[i][3])/2
#y=1.5
Z=(z[1]+z[2])/2
print "the interpolated value when x=2.5 and y=1.5 is : %f" %(Z)

X=
[[0, 1, 4, 9, 16], [2, 3, 6, 11, 18], [6, 7, 10, 15, 22], [12, 13, 16, 21, 28], [18, 19, 22, 27, 34]]
the interpolated value when x=2.5 and y=1.5 is : 10.500000
In [ ]: