博主专业是电气,电气科学实验中,有时用到蒙特卡罗模拟,今天讲解下。
先附上个人电气博文目录传送门。见下。点击即可
学好电气全靠它,个人电气博文目录(持续更新中…)
蒙特卡罗模拟
蒙特卡罗法是一种以抽样和随机数的产生为基础的随机性方法,因此也称为随机抽样法、计算机随机模拟法等。蒙特卡罗方法的基本原理是通过数字模拟试验,得到所要求解的出现某种事件的概率,作为问题的近似解。很久之前人们就已经开始使用蒙特卡罗方法来解决问题了,早在17世纪,发生事件的概率是用发生事件的频率来决定的。在20世纪计算机的使用,使得这样的实验可以大量快速的进行模拟。近年来,随着蒙特卡罗方法不仅被用于数学、物理、化学等领域,还应用于核科学与技术、机械工程、天文学和仪器科学与技术等。该法产生了不同的分支,有直接蒙特卡罗法,动力学蒙特卡罗法等 。
本质就是通过模拟事件发生,然后通过计算某一特定事件发生的频率来计算得到某值(=概率)
本文通过两个简单算例感受下蒙特卡罗模拟。
算例1:计算π
import numpy as np
import matplotlib.pyplot as plt
from pylab import *
import matplotlib; matplotlib.use('TkAgg')
mpl.rcParams['font.sans-serif'] = ['SimHei']
mpl.rcParams['axes.unicode_minus'] = False
from matplotlib.patches import Circle#导入圆
n = 100000
# 投点次数
r = 1.0 # 半径
a, b = (0., 0.) # 圆心
# 圆的信息
x_min, x_max = a-r, a+r
y_min, y_max = b-r, b+r
# 正方形区域边界
x = np.random.uniform(x_min, x_max, n) # 均匀分布
y = np.random.uniform(y_min, y_max, n)
d = np.sqrt((x-a)**2 + (y-b)**2)# 计算点到圆心的距离
res = sum(np.where(d < r, 1, 0)) #np.where(d < r, 1, 0) 如果d小于r,返回1.否则返回0# 统计落在圆内的点的数目
pi = 4 * res / n #res / n 得到是 =圆面积/正方形面积 。4是本例中正方形面积。
#公式 π*r*r=正方形面积(4)*(落在圆内的个数/落在正方形中的个数)
print('蒙特卡罗模拟得到的π值: ', pi) #得到的π值
#画图
fig = plt.figure(figsize = (6,6))
axes = fig.add_subplot(1,1,1)
plt.scatter(x,y,marker='*',linewidth = 0.1)
plt.axis('equal')
circle = Circle(xy = (a,b),radius = r, alpha = 0.5 ,color = 'red')
axes.add_patch(circle)
plt.title('蒙特卡罗模拟得到的π值:%f ' %pi)
plt.show()
算例2:计算积分 y = sin(x)
import numpy as np
import matplotlib.pyplot as plt
from pylab import *
import matplotlib; matplotlib.use('TkAgg')
mpl.rcParams['font.sans-serif'] = ['SimHei']
mpl.rcParams['axes.unicode_minus'] = False
n = 10000#投点次数
#矩形边界
x_min, x_max = 0.0, 1.0
y_min, y_max = 0.0, 1.0
x = np.random.uniform(x_min, x_max, n) # 均匀分布
y = np.random.uniform(y_min, y_max, n)
def f(x):
return np.sin(x)
# 创建函数 y = sin(x)
res = sum(np.where(y < f(x), 1, 0)) #统计点在sin(x)下面的个数
integral = res / n
print('蒙特卡罗模拟得到的sin(x)积分值:%f: ', integral)
plt.scatter(x,y)#散点图
xi = np.linspace(0,1,100)
yi = np.sin(xi)
plt.plot(xi,yi,'--k')
plt.fill_between(xi, yi, 0, color ='red',alpha=0.5,label='area')
plt.grid()
plt.title('蒙特卡罗模拟得到的sin(x)积分值:%f: '%integral)
plt.show()
作者:电气-余登武