h = [0,0.8,0.4,0.2]#
I = [0,0.1728,1.0688,1.4848]#
E = [0,89.5,34.9,9.5]#
I1 = 4 * I[(2)] / 3 - I[(1)] / 3#
t = 1.640533#
et1 = t - I1 #
Et1 = et1 * 100/t#
print "Error of the improved integral for segment 1 and 2 = ",Et1,"%"
I2 = 4 * I[(3)] / 3 - I[(2)] / 3#
et2 = t - I2 #
Et2 = et2 * 100/t#
print "Error of the improved integral for segment 4 and 2 = ",Et2,"%"
I1 = 1.367467#
I2 = 1.623467#
I = 16 * I2 /15 - I1 / 15#
print "Obtained integral which is the correct answer till the seventh decimal",I
#f(x) = 0.2 + 25*x - 200*x**2 + 675*x**3 - 900*x**4 + 400*x**5
# for using two point gauss legendre formulae, the intervals have to be changed to -1 and 1
#therefore, x = 0.4 + 0.4 * xd
#thus the integral is transferred to
#(0.2 + 25*(0.4+0.4*x) - 200*(0.4 + 0.4*x)**2 + 675*(0.4 + 0.4*x)**3 - 900*(0.4 + 0.4*x)**4 + 400*(0.4 + 0.4*x)**5)*0.4
#for three point gauss legendre formulae
x1 = -(1/3) ** 0.5#
x2 = (1/3) ** 0.5#
I1 = (0.2 + 25*(0.4+0.4*x1) - 200*(0.4 + 0.4*x1)**2 + 675*(0.4 + 0.4*x1)**3 - 900*(0.4 + 0.4*x1)**4 + 400*(0.4 + 0.4*x1)**5)*0.4#
I2 = (0.2 + 25*(0.4+0.4*x2) - 200*(0.4 + 0.4*x2)**2 + 675*(0.4 + 0.4*x2)**3 - 900*(0.4 + 0.4*x2)**4 + 400*(0.4 + 0.4*x2)**5)*0.4#
I = I1 + I2#
print "Integral obtained using gauss legendre formulae =",I
t = 1.640533#
e = (t - I)*100/t#
print "error = ",e,"%"
# f(x) = 0.2 + 25*x - 200*x**2 + 675*x**3 - 900*x**4 + 400*x**5
# for using three point gauss legendre formulae, the intervals have to be changed to -1 and 1
#therefore, x = 0.4 + 0.4 * xd
#thus the integral is transferred to
#(0.2 + 25*(0.4+0.4*x) - 200*(0.4 + 0.4*x)**2 + 675*(0.4 + 0.4*x)**3 - 900*(0.4 + 0.4*x)**4 + 400*(0.4 + 0.4*x)**5)*0.4
#for three point gauss legendre formulae
x1 = -0.7745967#
x2 = 0#
x3 = 0.7745967#
c0 = 0.5555556#
c1 = 0.8888889#
c2 = 0.5555556#
I1 = (0.2 + 25*(0.4+0.4*x1) - 200*(0.4 + 0.4*x1)**2 + 675*(0.4 + 0.4*x1)**3 - 900*(0.4 + 0.4*x1)**4 + 400*(0.4 + 0.4*x1)**5)*0.4#
I2 = (0.2 + 25*(0.4+0.4*x2) - 200*(0.4 + 0.4*x2)**2 + 675*(0.4 + 0.4*x2)**3 - 900*(0.4 + 0.4*x2)**4 + 400*(0.4 + 0.4*x2)**5)*0.4#
I3 = (0.2 + 25*(0.4+0.4*x3) - 200*(0.4 + 0.4*x3)**2 + 675*(0.4 + 0.4*x3)**3 - 900*(0.4 + 0.4*x3)**4 + 400*(0.4 + 0.4*x3)**5)*0.4#
I = c0 * I1 + c1 * I2 + c2 * I3#
print "integral obtained using three point gauss legendre formulae = ",I
from math import exp
#f(t) = g*m*(int(0,10,(1-exp(-c*t/m))))/c
#for using gauss quadrature method, limits are changed to -1 to 1 by replcing x = 5 + 5*xd
#the new integral obtained is
#(1 - exp(-c*(5 + 5*x)/m ))*5
g = 9.8#
c = 12.5#
m = 68.1#
#for two point method
x1 = -(1/3)**0.5#
x2 = (1/3)**0.5#
I1 = g*m*(1 - exp(-c*(5 + 5*x1)/m ))*5 / c#
I2 = g*m*(1 - exp(-c*(5 + 5*x2)/m ))*5 / c#
I = I1 + I2#
print "integral by two point method = ",I
x1 = -0.7745967#
x2 = 0#
x3 = 0.7745967#
c0 = 0.5555556#
c1 = 0.8888889#
c2 = 0.5555556#
I1 = g*m*(1 - exp(-c*(5 + 5*x1)/m ))*5 / c#
I2 = g*m*(1 - exp(-c*(5 + 5*x2)/m ))*5 / c#
I3 = g*m*(1 - exp(-c*(5 + 5*x3)/m ))*5 / c#
I = c0*I1 + c1 * I2 + c2 * I3#
print "integral by three point method =",I
x1 = -0.861136312#
x2 = -0.339981044#
x3 = 0.339981044#
x4 = 0.861136312#
c1 = 0.3478548#
c2 = 0.6521452#
c3 = 0.6521452#
c4 = 0.3478548#
I1 = g*m*(1 - exp(-c*(5 + 5*x1)/m ))*5 / c#
I2 = g*m*(1 - exp(-c*(5 + 5*x2)/m ))*5 / c#
I3 = g*m*(1 - exp(-c*(5 + 5*x3)/m ))*5 / c#
I4 = g*m*(1 - exp(-c*(5 + 5*x4)/m ))*5 / c#
I = c1*I1 + c2 * I2 + c3 * I3 + c4 * I4#
print "integral by four point method =",I
x1 = -0.906179846#
x2 = -0.538469310#
x3 = 0#
x4 = 0.538469310#
x5 = 0.906179846
c1 = 0.2369269#
c2 = 0.4786287#
c3 = 0.5688889#
c4 = 0.4786287#
c5 = 0.2369269#
I1 = g*m*(1 - exp(-c*(5 + 5*x1)/m ))*5 / c#
I2 = g*m*(1 - exp(-c*(5 + 5*x2)/m ))*5 / c#
I3 = g*m*(1 - exp(-c*(5 + 5*x3)/m ))*5 / c#
I4 = g*m*(1 - exp(-c*(5 + 5*x4)/m ))*5 / c#
I5 = g*m*(1 - exp(-c*(5 + 5*x5)/m ))*5 / c#
I = c1*I1 + c2 * I2 + c3 * I3 + c4 * I4 + c5 * I5#
print "integral by five point method =",I
from math import exp,pi
#N(x) = (int(-infinity,-2,exp(-(x**2)/2)) + int(-2,1,exp(-(x**2)/2)))/(2*pi)**0.5
#first integral can be solved as
#int(-infinity,-2,exp(-(x**2)/2)) = int(-0.5,0,exp(-1/(2*t**2))/t**2)
h = 1/8#
#int(-0.5,0,exp(-1/(2*t**2))/t**2) = h*(f(x-7/16) + f(x-5/16) + f(x-3/16) + f(x-1/16))
t1 = -7/16#
t2 = -5/16#
t3 = -3/16#
t4 = -1/16#
m1 = exp(-1/(2*t1**2))/t1**2#
m2 = exp(-1/(2*t2**2))/t2**2#
m3 = exp(-1/(2*t3**2))/t3**2#
m4 = exp(-1/(2*t4**2))/t4**2#
I1 = h*(m1 + m2 + m3 + m4)#
print "the value of first integral = ",I1
#simpsons 1/3rd sule is applied for the second integral
h1 = 0.5#
x1 = -2#
x2 = -1.5#
x3 = -1#
x4 = -0.5#
x5 = 0#
x6 = 0.5#
x7 = 1#
n1 = exp(-(x1**2)/2)#
n2 = exp(-(x2**2)/2)#
n3 = exp(-(x3**2)/2)#
n4 = exp(-(x4**2)/2)#
n5 = exp(-(x5**2)/2)#
n6 = exp(-(x6**2)/2)#
n7 = exp(-(x7**2)/2)#
I2 =(1-(-2)) * (n1 + 4 *(n2 + n4 + n6) + 2*(n3 + n5) + n7)/(18)#
print "The value of second integral = ",I2
f = (I1 + I2)/(2 * pi)**0.5#
print "Therefore the final result can be computed as :",f
N = 0.8413#
e = (N - f) * 100 / N#
print "error = ",e,"%"