#linear splines
#example 5.1
#page 182
from numpy import matrix
X=matrix([[1],[2], [3]])
y=matrix([[-8],[-1],[18]])
m1=(y[1][0]-y[0][0])/(X[1][0]-X[0][0])
m2=(y[2][0]-y[1][0])/(X[2][0]-X[1][0])
def s1(x):
return y[0][0]+(x-X[0][0])*m1
def s2(x):
return y[1][0]+(x-X[1][0])*m2
print "the value of function at 2.5 is %0.2f: " %(s2(2.5))
#cubic splines
#example 5.3
#page 188
from numpy import matrix
import math
X=matrix([[1],[2],[3]])
y=matrix([[-8],[-1],[18]])
M1=0
M2=8
M3=0
h=1
#deff('y=s1(x)','y=3*(x-1)^3-8*(2-x)-4*(x-1)')
def s1(x):
return 3*math.pow(x-1,3)-8*(2-x)-4*(x-1)
#deff('y=s2(x)','y=3*(3-x)^3+22*x-48');
def s2(x):
return 3*math.pow(3-x,3)+22*x-48
h=0.0001
n=2.0
D=(s2(n+h)-s2(n))/h;
print " y(2.5)=%f" %(s2(2.5))
print "y1(2.0)=%f" %(D)
#cubic spline
#example 5.4
#page 189
from numpy import matrix
import math
x=matrix([[0],[math.pi/2],[math.pi]])
y=matrix([[0],[1],[0]])
h=x[1][0]-x[0][0]
M0=0
M2=0
M1=((6*(y[0][0]-2*y[1][0]+y[2][0])/math.pow(h,2))-M0-M2)/4
X=math.pi/6.0
s1=(math.pow(x[1][0]-X,3)*(M0/6)+math.pow(X-x[0][0],3)*M1/6+(y[0][0]-math.pow(h,2)*M0/6)*(x[1][0]-X)+(y[1][0]-math.pow(h,2)*M1/6)*(X-x[0][0]))/h;
x=matrix([[0],[math.pi/4], [math.pi/2], [3*math.pi/4], [math.pi]])
y=matrix([[0], [1.414], [1] ,[1.414]])
M0=0
M4=0
A=matrix([[4, 1, 0],[1, 4, 1],[0 ,1 ,4]]) #calculating value of M1 M2 M3 by matrix method
C=matrix([[-4.029],[-5.699],[-4.029]])
B=A.I*C
print "M0=%f\t M1=%f\t M2=%f\t M3=%f\t M4=%f\t\n\n" %(M0,B[0][0],B[1][0],B[2][0],M4)
h=math.pi/4;
X=math.pi/6;
s1=(-0.12408*math.pow(X,3)+0.7836*X)/h;
print "the value of sin(pi/6) is:%f" %(s1)
#cubic spline
#example 5.5
#page 191
import math
from numpy import matrix
x=[1,2,3]
y=[6,18,42]
m0=40
s1=0
m1=(3*(y[2]-y[0])-m0)/4
X=0
s1=m0*((x[1]-X)**2)*(X-x[0])-m1*((X-x[0])**2)*(x[1]-X)+y[0]*((x[1]-X)**2)*(2*(X-x[0])+1)+y[1]*((X-x[0])**2)*(2*(x[1]-X)+1)
print "s1= %f+261*x-160X^2+33X^3" %(s1)
#surface fitting by cubic spline
#example 5.7
#page 195
import math
from numpy import matrix
def L0(y):
return math.pow(y,3)/4-5*y/4+1
def L1(y):
return (math.pow(y,3)/2)*-1+3*y/2
def L2(y):
return math.pow(y,3)/4-y/4
A=[ [1,2,9], [2,3,10], [9,10,17] ]
x=0.5
y=0.5
S=0.0
S=S+L0(x)*(L0(x)*A[0][0]+L1(x)*A[0][1]+L2(x)*A[0][2])
S=S+L1(x)*(L0(x)*A[1][0]+L1(x)*A[1][1]+L2(x)*A[1][2])
S=S+L2(x)*(L0(x)*A[2][0]+L1(x)*A[2][1]+L2(x)*A[2][2])
print "approximated value of z(0.5 0.5)=%f\n\n" %(S)
print " error in the approximated value : %f" %((abs(1.25-S)/1.25)*100)
#cubic B-splines
#example 5.8
#page 200
import math
k=[0.0, 1, 2, 3, 4]
pi=[0.0, 0, 4, -6, 24]
x=1
S=0
for i in range(2,5):
S=S+math.pow(k[i]-x,3)/(pi[i])
print "the cubic splines for x=1 is %f\n\n" %(S)
S=0
x=2
for i in range(2,5):
S=S+math.pow(k[i]-x,3)/(pi[i])
print "the cubic splines for x=2 is %f\n\n" %(S)
#cubic B-spline
#example 5.9
#page 201
k=[0, 1, 2, 3, 4];
x=1 #for x=1
s11=0
s13=0
s14=0
s24=0
s12=1/(k[2]-k[1])
s22=((x-k[0])*s11+(k[2]-x)*s12)/2.0 #k[2]-k[0]=2
s23=((x-k[1])*s11+(k[3]-x)*s13)/(k[3]-k[1])
s33=((x-k[0])*s22+(k[3]-x)*s23)/(k[3]-k[0])
s34=((x-k[1])*s23+(k[4]-x)*s24)/(k[4]-k[1])
s44=((x-k[0])*s33+(k[4]-x)*s34)/(k[4]-k[0])
print "s11=%f\t s22=%f\t s23=%f\t s33=%f\t s34=%f\t s44=%f\n\n" %(s11,s22,s23,s33,s34,s44)
x=2 #for x=2
s11=0
s12=0
s14=0
s22=0
s13=1/(k[3]-k[2])
s23=((x-k[1])*s12+(k[3]-x)*s13)/2.0 # k[3]-k[1]=2
s24=((x-k[2])*s13+(k[4]-x)*s14)/(k[2]-k[0])
s33=((x-k[0])*s22+(k[3]-x)*s23)/(k[3]-k[0])
s34=((x-k[1])*s23+(k[4]-x)*s24)/(k[4]-k[1])
s44=((x-k[0])*s33+(k[4]-x)*s34)/(k[4]-k[0])
print "s13=%f\t s23=%f\t s24=%f\t s33=%f\t s34=%f\t s44=%f\n\n" %(s13,s23,s24,s33,s34,s44)