在 Python 中使用 Monte Carlo 积分计算体积和质心

计算科学 Python 蒙特卡洛
2021-11-27 18:47:27

特别是，我想专注于找到体积，因为我需要它来开始解决质心问题V

这个同质体（环面部分）定义为3D

x2+(y2+z2−5)2<4
其中，

y<1,z>2

我将形状限制在一个矩形区域，其中3D

−2<x<2−7<y<12<z<7

体积定义如下V=∭dxdydz

下面是我尝试对积分进行编码以找到体积。

我是蒙特卡洛编码的新手，所以任何反馈都将不胜感激。

import numpy as np
from scipy import random 

#Limits of Integration
x_min = -2
x_max = 2 

y_min = -7
y_max = 1 

z_min = 2
z_max = 7 

#Number of Iterations 
N = 1000 

Xrand = random.uniform(x_min,x_max,N)
Yrand = random.uniform(y_min,y_max,N)
Zrand = random.uniform(z_min,z_max,N)

for i in range(len(Xrand)):
    Xrand[i] = random.uniform(x_min,x_max)
for i in range(len(Yrand)):
    Yrand[i] = random.uniform(y_min,y_max)    
for i in range(len(Zrand)):
    Zrand[i] = random.uniform(z_min,z_max)   

def func(x,y,z):
    return x**2 + (np.sqrt(y**2 + z**2)-5)**2

integral = 0.0

for i in range(N):
    integral += func(Xrand[i],Yrand[i],Zrand[i])
    
Volume =(1/N)*float(N)*integral 
print(Volume)
2个回答

这是一个更正和稍微改进的代码，在这里设置用于计算完整的环面体积以验证结果。

import numpy as np
from scipy import random 

def func(x,y,z):
    return x**2 + (np.sqrt(y**2 + z**2)-5)**2

def MCvolume(N=1000):
    #Limits of Integration
    x_min = -2
    x_max = 2 

    y_min = -7
    y_max = 7 #1 

    z_min = -7 #2
    z_max = 7 

    Xrand = random.uniform(x_min,x_max,N)
    Yrand = random.uniform(y_min,y_max,N)
    Zrand = random.uniform(z_min,z_max,N)

    for i in range(len(Xrand)):
            Xrand[i] = random.uniform(x_min,x_max)
    for i in range(len(Yrand)):
            Yrand[i] = random.uniform(y_min,y_max)    
    for i in range(len(Zrand)):
            Zrand[i] = random.uniform(z_min,z_max)   

    integral = 0.0

    for i in range(N):
        if (func(Xrand[i],Yrand[i],Zrand[i]) < 4.0):
            integral += 1.0
    
    VolumeBox=(x_max-x_min)*(y_max-y_min)*(z_max-z_min)
    Volume = VolumeBox*integral/float(N) 

    return Volume


#-stats from M trials, in each one using N random points
Mtry=10
Nrand=1000
vols=np.zeros(Mtry)

for k in range(0,Mtry):
    vols[k]=MCvolume(N=Nrand)

print("For Mtry, Nrand =", Mtry, ", ", Nrand, "; mean=", np.mean(vols), ", standard deviation=", np.std(vols))
print("Exact torus volume = ", 2*np.pi*5.0*np.pi*2.0**2)

您缺少构成蒙特卡洛积分的组成部分。

整合该区域（在您的情况下为体积）的步骤是在包括所需体积的整个域中生成随机样本（您已经这样做了）。现在，你想总结你想要整合的区域内的所有点，并丢弃那些不在其中的点。如我所见，您目前不检查您的随机样本是否在环面内！两个检查它​​们是否位于其中，您只需测试您所说的环面定义是否适用于样本。

其次，您必须考虑要分配给构成积分的那些点击的权重。您正在通过您的域分布总共 1000 个随机点，这些点由 x_max、x_min 等定义。然后权重应该是：1/Samples * Volume_sampled。当您在代码中修复这两个问题时，您应该得到正确的体积积分。

另外，我建议你大幅增加样本量，MC 方法收敛缓慢。

其它你可能感兴趣的问题
上一篇Navier--Stokes 方程的不同形式 下一篇如何检测由于浮点加法和乘法的舍入而导致的精度损失？