Implementing composite Gauss quadrature in Python

I want to implement the composite Gaussian quadrature in Python to evaluate the integral ∫₀¹ e^x2 dx. Evaluting this using Python’s quad command, I get ∫₀¹ e^x2 dx ≈ 1.46

Below is my attempt at implementing this in Python. What I expect is that as n gets larger, the closer the quadrature gets to the ‘real’ integral. However, as I vary n, the results gets smaller and smaller. What mistake am I making?

def gauss(f,a,b,n):
    h = float(b-a)/n
    [x,w] = np.polynomial.legendre.leggauss(n)
    result =  0
    for i in range(n):
        result += w[i] * f(x[i])
    result *= h
    return result


    for i in range(1,10):  
        print(gauss(exp,0,1,i))
    
    2.0
    1.3956124250860895
    0.9711551112557443
    0.731135234765899
    0.5850529284514102
    0.4875503086966867
    0.41790049038666144
    0.36566293624426005
    0.32503372130693325

Answer

Here is the working code based on your attempt:

import numpy as np

def gauss(f,a,b,n):
    half = float(b-a)/2.
    mid = (a+b)/2.
    [x,w] = np.polynomial.legendre.leggauss(n)
    result =  0.
    for i in range(n):
        result += w[i] * f(half*x[i] + mid)
    result *= half
    return result

def fun(x):
    return np.exp(x**2)

for i in range(1,10):  
    print(gauss(fun,0,1,i))

The problem was that the integral has to be rescaled to the domain of the Legendre polynomials you are using. Also, you have not used the quadrature rule correctly (it does not use the stepsize h that you have in your code.) The correct rule is:

∫_a^b f(x) dx ≈ (b-a)/2 ∑_i=1ⁿw_i f((b-a)/2ξ_i + (a+b)/2).

With this the code above outputs:

1.2840254166877414
1.4541678892391303
1.462409711477322
1.462646815656646
1.4626516680186825
1.4626517449041974
1.4626517458962964
1.4626517459070796
1.462651745907181