интеграция 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

Это называется multidimensional numerical integration. Чтобы выполнить такую ​​функцию на 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 

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

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

UPDATE: Вот как он будет использоваться:

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)