integration 3 dimensions by generalizing the method of rectangle

Question

def RectanglesPointMilieu(f,a,b,n):
    interval = 1.* (b-a) / n
    sumfct = 0
    for i in np.arange(a,b,interval):
            sumfct += f(i + interval/2.)
    return interval * sumfct

how do i change it into 3 dimensions calculation?

def rec3(f,X1,X2,X3,a1,a2,a3,b1,b2,b3,N1,N2,N3):
    interval1 = 1.* (b1-a1) / N1
    interval2 = 1.* (b2-a2) / N2
    interval3 = 1.* (b3-a3) / N3
    Sum1 = 0
    Sum2 = 0
    Sum3 = 0
    for i in np.arange(a1,b1,interval1):
        Sum1 += f((X1[i]+ interval1))
    for i in np.arange(a2,b2,interval2):
        Sum2 += f((X2[i]+ interval2))
    for i in np.arange(a3,b3,interval3):
        Sum3 += f((X3[i]+ interval3))

    return interval * float(Sum)

i did something like this but i'm just so lost and confused...i don't know how to continue...[I'm new at this]

def rec3(f,X1,X2,X3,a1,a2,a3,b1,b2,b3,N1,N2,N3):
    interval1 = 1.* (b1-a1) / N1
    interval2 = 1.* (b2-a2) / N2
    interval3 = 1.* (b3-a3) / N3
    Sum = 0
    for i in np.arange(a1,b1,interval1):
        for j in np.arange(a2,b2,interval2):
            for k in np.arange(a3,b3,interval3):
                Sum += f((X1[i]+ interval1),(X2[j]+ interval2),(X3[k]+ interval3))
    return interval1 * interval2 * interval3 * float(Sum)

with TypeError: 'numpy.ndarray' object is not callable

#def TEST_Q2():
    # Créer des tableaux
N = 1E5
X1 = rd.normal(0,1,N)
X2 = rd.normal(0,1,N)
X3 = rd.normal(0,1,N)
XX1 = X1 * np.sqrt(2)
XX2 = X2 * np.sqrt(2)
XX3 = X3 * 2.


def fct(a,b,c): # [x1**2 * x2**4 * exp(-x1**2 -x2**2 -2x3**2)]
    return ((a**2. * b**4.)/32.) * np.exp(-a**2./2.-b**2./2.-c**2./2.)




# Calculer l’intégrale
F = fct(XX1,XX2,XX3)
print rec3(F,XX1,XX2,XX3,0,0,0,1,1,1,N,N,N)


# Rlt = 11.8416033988

i use this as my Test

Answer 1

This is called multidimensional numerical integration. In order to perform such on an n-dimensional function, you will need to use n nested loops. In your case where n=3, the naive approach is to simply do something like:

def rec3(f,a1,a2,a3,b1,b2,b3,N1,N2,N3):
    interval1 = 1.* (b1-a1) / N1
    interval2 = 1.* (b2-a2) / N2
    interval3 = 1.* (b3-a3) / N3
    Sum = 0

    for x in np.arange(a1,b1,interval1):
        for y in np.arange(a2,b2,interval2):
            for z in np.arange(a3,b3,interval3):
                Sum += f(x+interval1/2., y+interval2/2., z+interval3/2.)

    return abs(interval1*interval2*interval3)*Sum

I don't know python very well, so there may be a syntax error in there, so use with care.

There are other approaches which may yield better results. You can try a few of the links in the search I suggested above to see if there is something more appropriate for your particular need.

UPDATE: Here is how it would be used:

N = 1E3

def fct(a,b,c):
    return ((a**2. * b**4.)/32.) * np.exp(-a**2./2.-b**2./2.-c**2./2.)

print rec3(fct,0,0,0,1,1,1,N,N,N)