python实现蒙特卡罗模拟法的实践

#模拟求近似圆周率
import random
import numpy as np
import matplotlib.pyplot as plt
#解决图标题中文乱码问题
import matplotlib as mpl
mpl.rcParams['font.sans-serif'] = ['SimHei']  # 指定默认字体
mpl.rcParams['axes.unicode_minus'] = False  # 解决保存图像是负号'-'显示为方块的问题
#进入正题
r=random.random()
num=range(0,200000,10)
mypi=np.ones((1,len(num)))
for j in range(len(num)):
    # 投点次数
    n = 10000
    # 圆的半径、圆心
    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))
    #计算pi的近似值（Monte Carlo:用统计值去近似真实值）
    mypi[0,j] = 4 * res / n
plt.plot(range(1,len(mypi[0])+1),mypi[0],'.-') 
plt.grid(ls=":",c='b',)#打开坐标网格
plt.axhline(y=np.pi,ls=":",c="yellow")#添加水平直线
# plt.axvline(x=4,ls="-",c="green")#添加垂直直线
plt.legend(['模拟', '实际'], loc='upper right', scatterpoints=1)
plt.title("近似圆周率")
plt.show()

#估算定积分
import random
import numpy as np
import matplotlib.pyplot as plt
#解决图标题中文乱码问题
import matplotlib as mpl
mpl.rcParams['font.sans-serif'] = ['SimHei']  # 指定默认字体
mpl.rcParams['axes.unicode_minus'] = False  # 解决保存图像是负号'-'显示为方块的问题
#进入正题
r=random.random()
num=range(1,10**6,500)
s = np.ones((1,len(num)))
for j in range(len(num)):
    n = 30000
    #矩形边界区域
    x_min, x_max = 0.0, 1.0
    y_min, y_max = 0.0, 1.0
    #在矩形区域内随机投点x
    x = np.random.uniform(x_min, x_max, n)#均匀分布
    y = np.random.uniform(y_min, y_max, n)
    #统计落在函数y=x^2下方的点
    res = sum(np.where(y < x**2, 1 ,0))
    #计算定积分的近似值
    s[0,j] = res / n
plt.plot(range(1,len(s[0])+1),s[0],'.-') 
plt.grid(ls=":",c='b',)#打开坐标网格
plt.axhline(y=1/3,ls=":",c="red")#添加水平直线
# plt.axvline(x=4,ls="-",c="green")#添加垂直直线
plt.legend(['模拟', '实际1/3'], loc='upper right', scatterpoints=1)
plt.title("估算定积分")
plt.show()

# 求解整数规划
import random
import numpy as np
import time
import matplotlib.pyplot as plt
#解决图标题中文乱码问题
import matplotlib as mpl
mpl.rcParams['font.sans-serif'] = ['SimHei']  # 指定默认字体
mpl.rcParams['axes.unicode_minus'] = False  # 解决保存图像是负号'-'显示为方块的问题
#进入正题
time_start=time.time() #计时开始
p0=0
for i in range(10**7):
    x=np.random.randint(0,99,(1,3))
    f=2*x[0,0]+3*x[0,0]**2+3*x[0,1]+x[0,1]**2+x[0,2]
    g=[
        x[0,0]+2*x[0,0]**2+x[0,1]+2*x[0,1]**2+x[0,2],
        x[0,0]+x[0,0]**2+x[0,1]+x[0,1]**2-x[0,2],
        2*x[0,0]+x[0,0]**2+2*x[0,1]+x[0,2],
        x[0,0]+2*x[0,1]**2
    ]
    if g[0]<=100 and g[1]<=500 and g[2]<=400 and g[3]>=10:
        if p0<f:
            x0=x
            p0=f
print('最优解：',x0)
print('最优值：',p0)
time_end=time.time() #计时结束
print('用时：',time_end-time_start)