import numpy as np
from scipy import integrate

##  (x1)^2(x2)^2(x3)^2 在 [0,1] 的积分

a1, b1 = 1, 2
a2, b2 = 3, 4
a3, b3 = 1, 5


# 三重积分计算
def f(x1, x2, x3):
    return x1 ** 2 * x2 ** 2 * x3 ** 2


I_exact, Error = integrate.tplquad(f, a1, b1, a2, b2, a3, b3)
np.random.seed(1)

# 平均值法
N = 10000
x1_sample = a1 + (b1 - a1) * np.random.rand(N)
x2_sample = a2 + (b2 - a2) * np.random.rand(N)
x3_sample = a3 + (b3 - a3) * np.random.rand(N)

h_x = f(x1_sample, x2_sample, x3_sample)
I_approx_stat = (b3 - a3) * (b2 - a2) * (b1 - a1) / N * np.sum(h_x)

# 数值积分
M = 200
h1 = (b1 - a1) / (M - 1)
h2 = (b2 - a2) / (M - 1)
h3 = (b3 - a3) / (M - 1)
x1 = np.linspace(a1, b1, M)
x2 = np.linspace(a2, b2, M)
x3 = np.linspace(a3, b3, M)
x1_mesh, x2_mesh, x3_mesh = np.meshgrid(x1, x2, x3)
I_approx_rec = np.sum(f(x1_mesh, x2_mesh, x3_mesh) * h1 * h2 * h3)

print('多重积分值：', I_exact)
print('\n平均值法结果：', I_approx_stat)
print('\n数值积分结果：', I_approx_rec)