再说圆的面积-蒙特卡洛（蒙特卡洛方法求圆周率的matlab程序）

在微积分-圆的面积和周长(1)介绍微积分方法求解圆的面积，本文使用蒙特卡洛方法求解圆面积。

取（0，1）* (0,1)区间，也就是单位圆第一象限的端点区间对应的正方形区间。下面是计算机给出的结果：很显然第一象限的面积是1/4单位圆面积，即Π/4

在这里插入图片描述

本文使用蒙特卡洛方法对圆面积进行求解。 最小标椎随机生成器算法：

xi = axi-1(mod m)
ui = xi/m
其中m=2^31-1,   a=7^5=16807

这里使用最小标椎生成器算法生成随机数。代码中（ux,uy） = (ui, ui+1)。

# #!/usr/bin/python
import matplotlib.pyplot as plt
import numpy as np
xi = 3
count=0
list_x = []
list_y = []
for i in range(1000001):
    xi = 16807*xi%( 2**31-1)
    ux = xi/(2**31-1)
    # print(ux)
    xi = 16807*xi%( 2**31-1)
    uy = xi/(2**31-1)
    if (ux**2+uy**2 <= 1):
        list_x.append(ux)
        list_y.append(uy)
        count = count + 1
print(count/1000000)

x = np.array(list_x)
y = np.array(list_y)

plt.scatter(x, y)
plt.show()

在这里插入图片描述

运行结果：0.785573

直接使用python的random进行求解：

#!/usr/bin/python
# Write Python 3 code in this online editor and run it.
import random
print(random.uniform(0, 1))
count = 0
for i in range(1000001):
    ux = random.uniform(0, 1)
    uy = random.uniform(0, 1)
    if (ux**2+uy**2 <= 1):
        count = count + 1
print(count/1000000)

运行结果：0.786203

从上面的结果来看，第一种方法的效果更好一点。

蒙特卡洛方法看起来高大上的东西，本质上就是概率问题。

上面的代码是蒙特卡洛二型方法，这个代码是一型方法，利用平均值求面积。

#!/usr/bin/python
# Write Python 3 code in this online editor and run it.
xi = 3
result=0
for i in range(1000001):
        xi = 16807*xi%( 2**31-1)
        ux = xi/(2**31-1)
        # print(ux)
        result += (1-ux**2)**0.5
print(result/1000000)

运行结果：0.7852800345401167

更多内容，欢迎关注我的微信公众号：半夏之夜的无情剑客