用蒙特卡洛方法计算π，计算e,计算γ

最新推荐文章于 2023-04-25 11:40:39 发布

aprilzj123

最新推荐文章于 2023-04-25 11:40:39 发布

阅读量2.7k

点赞数 4

分类专栏：小白算法解析文章标签：几何学算法

本文链接：https://blog.csdn.net/aprilzj123/article/details/122099682

版权

小白算法解析专栏收录该内容

1 篇文章 0 订阅

订阅专栏

用蒙特卡洛方法计算π，e和欧拉常数γ，并可视化

算法设计与分析的作业，还挺有意思的，给大家分享下

1.用蒙特卡洛方法计算π

import numpy as np
import math 
import matplotlib.pyplot as plt

list=[]
N=[]

for n in range(9999990,10000000): #为了看清楚概率的变化，可以做一个循环，n自动取整数
   x=np.random.uniform(-0.5,0.5,n) #产生的随机点的坐标定位。产生随机点的功能在numpy里
   y=np.random.uniform(-0.5,0.5,n) #产生的随机点的坐标定位。产生随机点的功能在numpy里
   
    #统计落在圆形里面的点
   pos=sum(np.where(x**2+y**2<0.25,1,0))#通过求和来得到落在圆里面的点的总数，条件写在where里面，如果满足条件返回1，不满足返回0
   pis=pos/n/0.25#计算π
   res=abs(pis-math.pi) #可以通过计算误差来看一下真实值和计算值的差距，在math库中有π，所以要引入math库
   list.append(res)
   N.append(n)
   print(pis)

print(list)
print(N)
plt.plot(N,list)

2.计算e的代码（这个来源于同济子豪兄，后续我的作业计算欧拉常数γ的代码就是参考的这位仁兄的代码）

# 蒙特卡洛法计算自然常数e——python编程及可视化

> 张子豪  同济大学

蒙特卡洛方法是一种用野蛮粗暴的蛮力对抗精致数学的一种计算思维，能够将复杂数学问题转化为简单粗暴的重复步骤，在工程上有很多应用。我还用蒙特卡洛方法计算了圆周率，请看我另一篇博客。

![蒙特卡洛方法计算自然常数e](https://upload-images.jianshu.io/upload_images/13714448-0e73a7a806a876ac.png?imageMogr2/auto-orient/strip%7CimageView2/2/w/1240)

![蒙特卡洛方法计算自然常数e](https://upload-images.jianshu.io/upload_images/13714448-800677c161c98f44.gif?imageMogr2/auto-orient/strip)



# 原理

![曲面四边形面积即为积分之后的值](https://upload-images.jianshu.io/upload_images/13714448-89504a4f4ce5f7a5.png?imageMogr2/auto-orient/strip%7CimageView2/2/w/1240)原理：用蒙特卡洛方法随机在左图矩形方格中撒点，统计y=1/x内外点的个数，
由几何概型，估算出曲线下曲面四边形的面积。
同时，由定积分可算出这部分面积为ln2，即e**(估算出的面积)== 2，即可求出e。

撒点越多，e的计算值也越来越趋近2.71828的真实值。

# 源代码

```python
# 张子豪 2019-3-14
import random
import matplotlib.pyplot as plt
import numpy as np


DARTS = 1024*1024 # 总共撒点的个数
counts = 0 # 落在曲线下方的点数
e = 0 # e的计算值
xs = [0,0]
ys = [0,0]

# 开始画左边的图：撒点估计曲线下方的面积
plt.subplot(121)
x = np.arange(0.5,2.5,0.001)
plt.ylim(0,1.25) # y轴坐标范围
plt.xlabel('x') # x轴标签
plt.ylabel('y') # y轴标签
plt.plot(x,1/x) # 绘制反比例函数曲线
plt.legend(loc=1) # 在右上角增加图例
plt.legend(['y = 1 / x']) # 图例的内容
plt.plot([1,1,2,2],[0,1,1,0],'r',linewidth=0.2) # 绘制撒点范围框

for i in range(DARTS):
    x = random.uniform(1,2)
    y = random.uniform(0,1)
    if y < 1/x: # 点落在曲线下方
        counts += 1
        plt.subplot(121)
        plt.plot(x,y,'g.')
    else: # 点落在曲线上方
        plt.subplot(121)
        plt.plot(x,y,'r.')
    if counts>0:
        e = pow(2,i/counts)

    # 开始画右边的图：e的计算值随投掷次数的关系
    plt.subplot(122)
    xs[0] = xs[1] # 上一个e值与下一个e值，通过xs与ys列表中的两个元素进行两点连线
    xs[1] = i
    ys[0] = ys[1]
    ys[1] = e
    plt.ylim(0,4.5) # y轴坐标范围
    plt.xlabel('Number of try') # x轴标签
    plt.ylabel('Estimation of e') # y轴标签
    plt.yticks(np.arange(0,4.5,0.5)) # y轴刻度线
    plt.title('e:{:.10f}\ncount:{}'.format(e,i)) # 图的标题动态更新
    plt.axhline(np.e,linewidth=0.05,color='r') # 绘制2.71828参考线
    plt.plot(xs,ys,'b--',linewidth=0.3) # 绘制e的计算值随撒点次数变化的曲线
    plt.ion() # 保持图像处于交互更新状态 
    plt.pause(0.2) # 控制撒点速度

2.1 可视化

蒙特卡洛方法计算自然常数e

3.我的原创——用蒙特卡洛算法计算欧拉常数γ

3.1 欧拉常数Euler-Mascheroni constant

欧拉常数的产生：
欧拉常数最先由瑞士数学家莱昂哈德·欧拉(Leonhard Euler)在1735年发表的文章 De Progressionibus harmonicus observationes 中定义。欧拉曾经使用C作为它的符号，并计算出了它的前6位小数。1761年他又将该值计算到了16位小数。1790年，意大利数学家马歇罗尼（Lorenzo Mascheroni)引入了γ作为这个常数的符号，并将该常数计算到小数点后32位。但后来的计算显示他在第20位的时候出现了错误。欧拉数以世界著名数学家欧拉名字命名。

欧拉常数的定义： 调和级数与自然对数的差值的极限。

在这里插入图片描述
其近似值为0.577 215 664 9…

3.2用蒙特卡洛算法计算定积分

蒙特卡洛方法的另一个重要应用就是求定积分。
在这里插入图片描述

当我们在[a,b]之间随机取一点x时，它对应的函数值就是f(x)。接下来我们就可以用来粗略估计曲线下方的面积，也就是我们需要求的积分值，当然这种估计（或近似）是非常粗略的。

在这里插入图片描述

3.3 算法设计

基本思想

4.代码部分

import random
import matplotlib.pyplot as plt
import numpy as np


DARTS = 1024*1024  # 总共撒点的个数
n = int(input('请输入一个较大的整数'))   # 要确保输入的整数足够大
countln = 0   # 落在对数函数下方的点数
countharm = 0  # 落在调和级数部分的点数
count1 = 0  # 落在单位正方形的点数
γ = 0   # γ的计算值
xs = [0, 0]
ys = [0, 0]


# 开始画左边的图：撒点估计面积
plt.subplot(121)
x = np.arange(1/n, 1, 1/n)
plt.ylim(0, n)  # y轴坐标范围
plt.xlabel('x')  # x轴标签
plt.ylabel('y')  # y轴标签
plt.plot(x, 1/x)  # 绘制反比例函数曲线
plt.legend(loc=1)  # 在右上角增加图例
plt.legend(['y = 1 / x'])  # 图例的内容
plt.plot([0, 0, 1, 1], [0, n, n, 0], 'r', linewidth=0.2)  # 绘制撒点范围框

for i in range(DARTS):
    x = random.uniform(0, 1)
    y = random.uniform(0, n)
for j in range(n-1):
    if j / n <= x <= (j + 1) / n and y <= n / (j + 1):  # 点落在调和级数部分
        countharm += 1
    plt.subplot(121)
    plt.plot(x, y, 'b.')
    if y < 1/x and x >= 1/n:  # 点落在对数曲线下方
        countln += 1
        plt.subplot(121)
        plt.plot(x, y, 'g.')
    if 0 <= x <= 1 and 0 <= y <= 1:
        count1 += 1
if countln > 0 and countharm > 0 and count1 > 0:
    γ = (countharm-countln)/count1


# 开始画右边的图：γ的计算值随投掷次数的关系
plt.subplot(122)
xs[0] = xs[1]   # 上一个γ值与下一个γ值，通过xs与ys列表中的两个元素进行两点连线
xs[1] = i
ys[0] = ys[1]
ys[1] = γ
plt.ylim(0, 1.0)   # y轴坐标范围
plt.xlabel('Number of try')   # x轴标签
plt.ylabel('Estimation of γ')   # y轴标签
plt.yticks(np.arange(0, 1.0, 0.5))   # y轴刻度线
plt.title('γ:{:.10f}\ncount:{}'.format(γ, i))   # 图的标题动态更新
plt.axhline(np.euler_gamma, linewidth=0.05, color='r')   # 绘制0.57721 参考线
plt.plot(xs, ys, 'b--', linewidth=0.3)   # 绘制γ的计算值随撒点次数变化的曲线
plt.ion()   # 保持图像处于交互更新状态
plt.pause(0.2)   # 控制撒点速度