Интеграция 3 измерения путем обобщения метода прямоугольника

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

как я могу изменить его на 3 измерения расчета?

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)

я сделал что-то подобное, но я просто потерян и растерян... я не знаю, как продолжить...[Я новичок в этом]

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)

с TypeError: объект 'numpy.ndarray' не вызывается

#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

я использую это как мой тест

Answer 1

Это называется многомерным численным интегрированием. Для того, чтобы выполнить такое на n-функция, вам нужно будет использовать n вложенные циклы. В вашем случае, где n=3Наивный подход - просто сделать что-то вроде:

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

Я не очень хорошо знаю Python, поэтому там может быть синтаксическая ошибка, так что будьте осторожны.

Есть и другие подходы, которые могут дать лучшие результаты. Вы можете попробовать несколько ссылок в поиске, который я предложил выше, чтобы увидеть, есть ли что-то более подходящее для ваших конкретных потребностей.

ОБНОВЛЕНИЕ: Вот как это будет использоваться:

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)