Numerical Integration

from matplotlib import pyplot as plt

plt.figure(figsize=(10,10))
file = open("input.txt","r")

n = file.readline().rstrip()
n = int(n)
X=[]
Y=[]
X_1=[]

for line in file:
    line  = line.rstrip().split()
    X.append(float(line[0]))
    Y.append(float(line[1]))

a=0
b=0
c=0
def Trapezoid(a,b):
    return 0.5*(X[b]-X[a])*(Y[a]+Y[b])

def Simpson13(a,b,c):
    return (X[c]-X[a])*(Y[a]+4*Y[b]+Y[c])/6.0

def Simpson38(a,b,c,d):
    return (X[d]-X[a])*(Y[a]+3*Y[b]+3*Y[c]+Y[d])/8.0

summ = 0
idx = 0
while idx<n-1:
    if((idx+3<=n-1) and abs((X[idx+3]-X[idx+2])-(X[idx+2]-X[idx+1]))<0.00001 and abs((X[idx+2]-X[idx+1])-(X[idx+1]-X[idx]))<0.00001):
        summ +=Simpson38(idx,idx+1,idx+2,idx+3)
        X_1=[]
        Y_1=[]
        X_1.extend([X[idx],X[idx+1],X[idx+2],X[idx+3]])
        Y_1.extend([Y[idx],Y[idx+1],Y[idx+2],Y[idx+3]])
        a = a+1
        idx = idx+3
        plt.plot(X_1,Y_1,color="black")
        plt.fill_between(X_1,Y_1,color="cyan",label="Simpson's 3/8 rule")


    elif((idx+2<=n-1) and abs((X[idx+2]-X[idx+1])-(X[idx+1]-X[idx]))<0.00001):
        summ +=Simpson13(idx,idx+1,idx+2)
        X_2=[]
        Y_2=[]
        X_2.extend([X[idx],X[idx+1],X[idx+2]])
        Y_2.extend([Y[idx],Y[idx+1],Y[idx+2]])
        b = b+1
        idx = idx+2
        plt.plot(X_2,Y_2,color="black")
        plt.fill_between(X_2,Y_2,color="pink",label="Simpson's 1/3 rule")
    elif(idx+1<=n-1):
        summ +=Trapezoid(idx,idx+1)
        X_3=[]
        Y_3=[]
        X_3.extend([X[idx],X[idx+1]])
        Y_3.extend([Y[idx],Y[idx+1]])
        c = c+1
        idx = idx+1
        plt.plot(X_3,Y_3,color="black")
        plt.fill_between(X_3,Y_3,color="yellow",label="Trapezoid rule")

print(X)
print(Y)
print("Trapezoid: ",c," intervals")
print("1/3 rule: ",b*2," intervals")
print("3/8 rule: ",a*3," intervals")
print("Integral Value: ",summ)
plt.xlabel(" x co-ordinates")
plt.ylabel(" f(x)")
plt.legend()
plt.grid("true")
plt.show()