Chapter-19 : Fourier Approximation¶
Ex:19.1 Pg: 529¶
In [5]:
from math import sin,cos,atan
def f(t):
y=1.7+cos(4.189*t+1.0472)
return y
deltat=0.15#
t1=0#
t2=1.35#
omega=4.189#
Del=(t2-t1)/9#
t=[]
for i in range(1,11):
t.append(t1+Del*(i-1))
sumy=0#
suma=0#
sumb=0#
y=[];a=[];b=[]
for i in range(1,11):
y.append(f(t[i-1]))
a.append(y[(i-1)]*cos(omega*t[(i-1)]))
b.append(y[i-1]*sin(omega*t[i-1]))
sumy=sumy+y[i-1]
suma=suma+a[i-1]
sumb=sumb+b[i-1]
A0=sumy/10#
A1=2*suma/10#
B1=2*sumb/10#
print "The least square fit is y=A0+A1*cos(w0*t)+A2*sin(w0*t), where"
print "A0=",A0
print "A1=",A1
print "B1=",B1
theta=atan(-B1/A1)#
C1=(A1**2 + B1**2)**0.5#
print "Alternatively, the least square fit can be expressed as"
print "y=A0+C1*cos(w0*t + theta), where"
print "A0=",A0
print "Theta=",theta
print "C1=",C1
print "Or"
print "y=A0+C1*sin(w0*t + theta + pi/2), where"
print "A0=",A0
print "Theta=",theta
print "C1=",C1
Ex:19.2 Pg: 532¶
In [6]:
a0=0#
#f(t)=-1 for -T/2 to -T/4
#f(t)=1 for -T/4 to T/4
#f(t)=-1 for T/4 to T/2
#ak=2/T* (integration of f(t)*cos(w0*t) from -T/2 to T/2)
#ak=2/T*((integration of f(t)*cos(w0*t) from -T/2 to -T/4) + (integration of f(t)*cos(w0*t) from -T/4 to T/4) + (integration of f(t)*cos(w0*t) from T/4 to T/2))
#Therefore,
#ak=4/(k*pi) for k=1,5,9,.....
#ak=-4/(k*pi) for k=3,7,11,.....
#ak=0 for k=even integers
#similarly we find the b's.
#all the b's=0
print "The fourier approximtion is:"
print "4/(pi)*cos(w)*t) - 4/(3*pi)*cos(3*(w)*t) + 4/(5*pi)*cos(5*(w)*t) - 4/(7*pi)*cos(7*(w)*t) + ....."
Ex:19.3 Pg: 550¶
In [16]:
from numpy import arange,log
%matplotlib inline
from matplotlib.pyplot import plot, title,xlabel,ylabel,show
x=arange(0.5,5.6,0.5)
y=[]
for i in range(1,12):
y.append(0.9846*log(x[i-1])+1.0004)
plot(x,y)
title("y vs x")
xlabel("x")
ylabel("y")
show()
Ex:19.6 Pg: 555¶
In [32]:
from numpy import nditer,mat,divide
x=[0.05, 0.12, 0.15, 0.3, 0.45 ,0.7 ,0.84 ,1.05]
y=[0.957, 0.851, 0.832, 0.72 ,0.583, 0.378, 0.295, 0.156]
sx=sum(x)#
sxx=0
for xx in x:
sxx+=xx*xx
sx3=0
for xx in x:
sx3+=xx*xx*xx
sx4=0
for xx in x:
sx4+=xx*xx*xx*xx
sx5=0
for xx in x:
sx5+=xx*xx*xx*xx*xx
sx6=0
for xx in x:
sx6+=xx*xx*xx*xx*xx*xx
n=8#
sy=sum(y)#
sxy=0
for xx,yy in nditer([x,y]):
sxy+=xx*yy
sx2y=0
for xx,yy in nditer([x,y]):
sx2y+=xx*xx*yy
sx3y=0
for xx,yy in nditer([x,y]):
sx3y+=xx*xx*xx*yy
m=mat([[n, sx, sxx ,sx3],[sx, sxx, sx3, sx4],[sxx, sx3, sx4, sx5],[sx3, sx4, sx5, sx6]])
p=mat([[sy],[sxy],[sx2y],[sx3y]])
a=divide(m,p)
print "The cubic polynomial is y=a0 + a1*x + a2*x**2 + a3*x**3, where a0, a1, a2 and a3 are"
print a