通过推广矩形方法进行 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)

与 类型错误：“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)