用蒙特卡洛方法计算π,计算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 算法设计

  1. 基本思想
    在这里插入图片描述
    在这里插入图片描述在这里插入图片描述
    在这里插入图片描述
    在这里插入图片描述

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)   # 控制撒点速度
版权声明:本文为博主原创文章,遵循 CC 4.0 BY-SA 版权协议,转载请附上原文出处链接和本声明。
本文链接:https://blog.csdn.net/aprilzj123/article/details/122099682

智能推荐

如何自定义NavigationBar的高度_rk framework navigationbarview 高度-程序员宅基地

文章浏览阅读1.2k次。UINavigationBar的高度在苹果官方的SDK中是固定的44个点,但是实际项目中我们却有可能遇到这样的情况,如下图:这样的一个UINavigationBar的高度达到了84个点,这就需要我们自定义系统自带的UINavigationBar的高度,但是系统并没有直截了当的方法来调整这个NavigationBar的Height,于是我进行了以下的尝试。在需要进行调整的那个ViewCo_rk framework navigationbarview 高度

python毕业设计作品基于django框架 教室实验室预约系统毕设成品(8)毕业设计论文模板_django项目毕设答辩ppt-程序员宅基地

文章浏览阅读359次。python毕业设计作品基于django框架 教室实验室预约系统毕设成品(8)毕业设计论文模板_django项目毕设答辩ppt

android wifi源码分析,android wifi打开过程源码解析及Wifi打开失败原因分析-程序员宅基地

文章浏览阅读464次。在android中wifi打开的状态从DISABLED-->ENABLING-->ENABLED1 WifiSettings.java--wifi界面wifi开关wifi开关定义在SettingsActivity.java中,传入WifiEnabler.java,并在WifiEnabler.java中响应。private SwitchBar mSwitchBar;public Swit..._android wifi打开过程源码解析及wifi打开失败原因分析

Java编程:删除 List 元素的三种正确方法_java synchronized list 元素删除-程序员宅基地

文章浏览阅读942次。删除 List 中的元素会产生两个问题:删除元素后 List 的元素数量会发生变化;对 List 进行删除操作可能会产生并发问题;_java synchronized list 元素删除

D语言游戏编程(11):D语言基础之模板和混入(mixin)技术_d语言 模板实例化-程序员宅基地

文章浏览阅读3.7k次。 D语言通过模板,很好的支持泛型编程。与C++的模板相比较,各有优略。总体上说,D语言的模板在很多方面还是很方便的。 D语言还支持模板的混入(mixin),简单的讲就是把模板实例化之后,将模板中的代码插入到当前的位置。这是一个非常方便的工具! 具体的,请看下面的演示代码。import std.stdio;void main()...{ tryTemplate();_d语言 模板实例化

目标 linux 服务器提权,史上最全Linux提权后获取敏感信息方法 (zhuan)-程序员宅基地

文章浏览阅读364次。(Linux)的提权是怎么一回事:收集 – 枚举,枚举和一些更多的枚举。过程 – 通过数据排序,分析和确定优先次序。搜索 – 知道搜索什么和在哪里可以找到漏洞代码。适应 – 自定义的漏洞,所以它适合。每个系统的工作并不是每一个漏洞“都固定不变”。尝试 – 做好准备,试验和错误。系统类型系统是什么版本?cat /etc/issuecat /etc主机上有哪些工作计划?crontab -lls -al..._linux服务器被提权如何解决

随便推点

C# 四个字节十六进制数和单精度浮点数之间的相互转化_bitconverter.tosingle-程序员宅基地

文章浏览阅读3.9k次。C# 四个字节十六进制数和单精度浮点数之间的相互转化即是所谓的IEEE754标准,这也是大多数硬件存储浮点数的标准。单精度浮点数占4个字节,表示范围为:在负数的时候是从 -3.402823E38 到 -1.401298E-45,而在正数的时候是从 1.401298E-45 到 3.402823E38 。一、在C#中的转换函数为:1,由四个字节的十六机制数组转浮点数: byte[] b..._bitconverter.tosingle

iOS[swift]防止单点手势连续快速触发_swift 修改手势响应时间-程序员宅基地

文章浏览阅读874次。按钮被连点,单点手势被连续触发有时候我们的APP会出现各种糟糕的现象(例如: 绑定事件是弹出个页面,连点之后连续弹出多个相同的页面,或者绑定某个指令,连点之后连续下发多个指令等等)等等...防止按钮的连点我已经在之前的文章中总结过了,本文总结的是防止自定义绑定的单点手势的连点:废话不说,老规矩上代码:import UIKitclass TapGestureManager:UITap..._swift 修改手势响应时间

使用 Nginx+SpringBoot+Redis 实现负载均衡以及session共享_spring boot 多节点连接redis 实现负载均衡-程序员宅基地

文章浏览阅读1.7k次。一、所需的环境Redis 存储Session信息两个Spring Boot 应用,均将session信息托管到上面的RedisNginx 配置负载均衡策略 反向代理到SpringBoot应用二、Nginx 配置在使用Nginx实现负载均衡主要使用了upstream 参数。我们在定义了一个名字为huizi的服务器池。我们本地分别启动了两个不同端口的应用。并分别指定不同应用的分..._spring boot 多节点连接redis 实现负载均衡

JLINK仿真器与ST-LINK仿真器的安装与配置.pdf-程序员宅基地

文章浏览阅读154次。JLINK仿真器与ST-LINK仿真器的安装与配置.pdf工欲善其事,。。。。。。stm32的开发环境搭建 观看地址说到仿真器,首先要了解一下JTAG。JTAG协议JTAG(Joint Test Action Group,联合测试行动小组)是一种国际标准测试协议(IEEE 1149.1兼容),主要用于芯片内部测试。现在多数的高级器件都支持JTAG协议,如ARM、DS..._jink仿真器输入和输出交叉

RDD与DataFrame与Dataset之间的关系及转换关系_dataframe、dataset、rdd之间的转换-程序员宅基地

文章浏览阅读403次。RDD与DataFrame与Dataset之间的转换关系:_dataframe、dataset、rdd之间的转换

关系型数据库和非关系型数据库的区别与联系_谈谈关系型数据库与非关系型数据库的联系和区别-程序员宅基地

文章浏览阅读2.2k次。数据库一、概念数据库是以一定方式储存在一起、能与多个用户共享、具有尽可能小的冗余度、与应用程序彼此独立的数据集合。数据库管理系统是一种操纵和管理数据库的大型软件,用于建立、使用和维护数据库,简称DBMS。它对数据库进行统一的管理和控制,以保证数据库的安全性和完整性。二、分类关系型数据库NOSQL三、NoSQL与关系型数据库的区别存储方式传统的关系型数..._谈谈关系型数据库与非关系型数据库的联系和区别

推荐文章

热门文章

相关标签