Simpson's rule 3/8 for n intervals in Python

Question

im trying to write a program that gives the integral approximation of e(x^2) between 0 and 1 based on this integral formula: Formula

i've done this code so far but it keeps giving the wrong answer (Other methods gives 1.46 as an answer, this one gives 1.006).

I think that maybe there is a problem with the two for cycles that does the Riemman sum, or that there is a problem in the way i've wrote the formula. I also tried to re-write the formula in other ways but i had no success

Any kind of help is appreciated.

import math
import numpy as np

def f(x):
    y = np.exp(x**2)
    return y

a = float(input("¿Cual es el limite inferior? 
"))
b = float(input("¿Cual es el limite superior? 
"))
n = int(input("¿Cual es el numero de intervalos? "))

x = np.zeros([n+1])
y = np.zeros([n])
z = np.zeros([n])

h = (b-a)/n

print (h)

x[0] = a
x[n] = b


suma1 = 0
suma2 = 0

for i in np.arange(1,n):
    x[i] = x[i-1] + h
    suma1 = suma1 + f(x[i])

alfa = (x[i]-x[i-1])/3

for i in np.arange(0,n):
    y[i] = (x[i-1]+ alfa)
    suma2 = suma2 + f(y[i])
    z[i] = y[i] + alfa


int3 = ((b-a)/(8*n)) * (f(x[0])+f(x[n]) + (3*(suma2+f(z[i]))) + (2*(suma1)))

print (int3)

user_2538939 · Answer

I'm not a math major but I remember helping a friend with this rule for something about waterplane area for ships.

Here's an implementation based on Wikipedia's description of the Simpson's 3/8 rule:

# The input parameters
a, b, n = 0, 1, 10

# Divide the interval into 3*n sub-intervals
# and hence 3*n+1 endpoints
x = np.linspace(a,b,3*n+1)
y = f(x)

# The weight for each points
w = [1,3,3,1]

result = 0
for i in range(0, 3*n, 3):
    # Calculate the area, 4 points at a time
    result += (x[i+3] - x[i]) / 8 * (y[i:i+4] * w).sum()

# result = 1.4626525814387632

user_13743233 · Answer

You can do it using numpy.vectorize (Based on this wikipedia post):

a, b, n = 0, 1, 10**6
h = (b-a) / n
x = np.linspace(0,n,n+1)*h + a

fv = np.vectorize(f)

(
3*h/8 * (
    f(x[0]) +
    3 * fv(x[np.mod(np.arange(len(x)), 3) != 0]).sum() + #skip every 3rd index
    2 * fv(x[::3]).sum() + #get every 3rd index
    f(x[-1])
    )
)
#Output: 1.462654874404461

If you use numpy's built-in functions (which I think is always possible), performance will improve considerably:

a, b, n = 0, 1, 10**6
x = np.exp(np.square(np.linspace(0,n,n+1)*h + a))

(
3*h/8 * (
    x[0] +
    3 * x[np.mod(np.arange(len(x)), 3) != 0].sum()+ 
    2 * x[::3].sum() + 
    x[-1]
    )
)
#Output: 1.462654874404461

Simpson's rule 3/8 for n intervals in Python

2 Answers