西安交大 数学实验五 matlab计算圆周率pi-课件·PPT

这里是一个示例课件的简要概述，用于介绍如何使用计算机来计算圆周率：
1. 引言：介绍圆周率的概念和重要性，以及计算圆周率的方法。

2. 数学公式：讲解数学公式中圆周率的含义和计算方法，如利用圆的周长和直径的比值。

3. 近似方法：介绍一些近似计算圆周率的方法，如阿基米德法、蒙特卡洛方法和无穷级数法。

4. 计算机的角色：讨论计算机在计算圆周率中的作用，包括高精度计算、并行计算和分布式计算。

5. 编程语言和软件工具：介绍一些常用的编程语言和软件工具，如Python、MATLAB和Mathematica，用于编写计算圆周率的程序。

6. 编程实例：提供一个简单的编程实例，演示如何使用编程语言计算圆周率的近似值。

7. 程序优化：讨论如何优化计算圆周率的程序，提高计算速度和精度，如使用多线程、并行计算和算法优化。

8. 应用领域：讨论计算圆周率在实际应用中的重要性，如在科学计算、密码学和通信领域的应用。

9. 结论和总结：总结课程内容，强调计算圆周率的重要性和应用前景。

这只是一个简要的概述，您可以根据具体需要和时间安排，进一步展开和深入讲解各个部分。

在课件中可以添加图表、示意图和实例代码等辅助材料，以帮助学生更好地理解和应用计算圆周率的方法。

Matlab圆周率1. 引言圆周率（π）是数学中的一个重要常数，表示圆的周长与直径之间的比值。

圆周率的精确值是一个无限不循环的小数，被广泛应用于数学、物理和工程等领域。

在本文中，我们将使用Matlab来计算圆周率，并探讨一些与圆周率相关的数值计算方法和应用。

2. 圆周率的定义与计算方法2.1 定义圆周率可以通过圆的周长与直径的比值来定义，即π = C/d，其中C表示圆的周长，d表示圆的直径。

2.2 常用计算方法2.2.1 几何法几何法是最直观的计算圆周率的方法之一。

它基于圆的定义，将一个正方形内切于一个圆，并通过计算正方形的边长与圆的直径的比值来逼近圆周率。

几何法的思想简单易懂，但精度较低，随着划分的正方形数量增加，计算结果会逼近真实值。

2.2.2 随机法随机法是一种通过生成随机点，并统计这些点落在圆内的比例来估计圆周率的方法。

随机法利用概率统计的思想，通过大量的随机模拟来逼近圆周率。

使用Matlab可以方便地生成随机数，并进行统计计算，从而得到较为准确的圆周率估计值。

2.2.3 数列法数列法是通过一些特定公式或数列来逼近圆周率的方法。

其中最著名的数列是莱布尼兹级数和无穷乘积法。

莱布尼兹级数是一个交替级数，通过逐项相加可以逼近圆周率的值。

无穷乘积法则是使用无穷乘积的形式来逼近圆周率，其中最著名的是狄利克雷（Dirichlet）级数。

这些数列法适用于在计算机中进行快速近似计算，但对于精确值的计算效果较差。

3. Matlab中的圆周率计算函数Matlab提供了多种计算圆周率的函数，方便了用户进行圆周率的计算和应用。

下面介绍几个常用的Matlab函数。

### 3.1 pi pi是Matlab内置的常量，表示圆周率的值。

使用pi函数可以直接获得圆周率的近似值，具有较高的精度。

3.2 montecarlo_pimontecarlo_pi是一个基于随机法的圆周率计算函数。

它通过生成随机点，并统计这些点落在圆内的比例来估计圆周率。

并在1994年计算到了4044000000位．它的另一 种形式是
426880 10005 . (6n)!(545140134 13591409) n (n! )3 (3n)!(640320)3n n 0
1995 年 ， 由 David Bailey,Peter Borwein 和 Simon Plouffe 共同发表了下面的圆周率计算公式 （简称BBP公式）
x5
b x
现对每个子区间再二等分，得到2n个子区间

b a
n 1 n x i 1 x i h f ( x )dx T2 n ( f (a ) 2 f ( x i ) 2 f ( ) f (b)) , 4 2 i 1 i 1 ba h , xk a kh k 0,1,2, , n n n n 1 h x i 1 x i 1 (Tn h f ( )) (Tn h f (a ih )) 2 2 2 2 i 1 i 1
2.圆周率的数值积分计算方法 1 1 1 1 0 1 x 2 dx 4 40 1 x 2 dx
表1给出的是不同等分区间 数计算出的 圆周率的近似值.
等分区间数n 10 20 50 100 500 1000 5000
圆周率的近似值
3.14242598500110 3.14180098689309 3.14162598692300 3.14160098692312 3.14159298692312 3.14159273692313 3.14159265692313

1
0
e
x2
dx
抛物线方法 y=inline('exp(-x.^2)'); quad(y,0,1) ans = 0.746826

由Cleve Moler和John Little于1980 年代初期开发，用于解决线性代数课 程的数值计算问题。
MATLAB的优势
易于学习、使用灵活、高效的数值计 算和可视化功能、强大的工具箱支持。
发展历程
从最初的数值计算工具，逐渐发展成 为一款功能强大的科学计算软件，广 泛应用于工程、科学、经济等领域。
MATLAB工作环境与界面
MATLAB工作环境
包括命令窗口、工作空间、命令历史窗口、当 前文件夹窗口等。
界面介绍
详细讲解MATLAB界面的各个组成部分，如菜 单栏、工具栏、编辑器窗口等。
基本操作
介绍如何在MATLAB环境中创建、保存、运行脚本和函数，以及如何进行基本 的文件操作。
基本数据类型与运算
矩阵大小
使用`size`函数获取矩阵的行数 和列数。
矩阵元素访问
通过下标访问矩阵元素，如 `A(i,j)`表示访问矩阵A的第i行第j 列元素。
矩阵基本操作
包括矩阵的加、减、数乘、转置 等操作。
矩阵运算及性质
矩阵乘法 满足乘法交换律和结合律，但不满足 乘法交换律。
矩阵的逆
对于方阵，若存在一矩阵B，使得 AB=BA=I（I为单位矩阵），则称B 为A的逆矩阵。
Hale Waihona Puke 03 数据分析与可视化数据导入、导出及预处理
数据导入
介绍如何使用MATLAB导入各种格式的数据文件， 如.csv、.txt、.xlsx等。
数据导出
讲解如何将MATLAB中的数据导出为常见的数据文件格式，以 便于数据共享和交换。
数据预处理
阐述数据清洗、数据变换、数据规约等预处理技术，为后续的数 据分析和可视化奠定基础。
01
02

nl n 2 , 特别取针的长度 l d / 2 时， md m
本实验涉及的方法
3、蒙特卡罗法——蒲丰(Buffon)掷针法
x表示针的中点与最近的一条平行线间的距离, 为针与平行线 的夹角.则0 x d 2 ,0 , 确定x 平面内一个矩形 .如图：
平面内一个矩形确定的夹角为针与平行线一条平行线间的距离表示针的中点与最近的针与平行线相交本实验涉及的方法3蒙特卡罗法蒲丰buffon掷针法针与平行线相交时针在平行线矩形区域内所组成的区域面积与矩形区域面积比即为针与直线相交的概率
实验三
圆周率的求法
问
题
• 设半径为1的圆的周长与圆的直径比为A; • 设半径为2的圆的周长与圆的直径比为B. • 请问A与B哪个大?
计算历程的几个阶段
3.分析方法时期 17世纪出现了数学分析，这锐利的工具 使得许多初等数学束手无策的问题迎刃 而解。 π 的计算历史也随之进入了一 个新的阶段。
计算历程的几个阶段
1593年，韦达给出
1671年，苏格兰数学家詹姆斯.格雷戈里 公开了他发现的公式： 3 5 7
arctan x x x3 x5 x7 (1 x 1)

故
2 nl md
本实验涉及的方法
3、蒙特卡罗法——蒲丰(Buffon)掷针法
• 随机生成一个针的中点与最近一条平行直线间的 距离x,满足0≤x≤d/2; • 随机生成针与平行线的夹角a,满足0≤a≤ ； • 判断x≤L/2*sin[a],若成立，则针与直线相交， 否则不相交。
本实验涉及的方法
计算历程的几个阶段
圆周长大于内接正四边形而小于外切正四边形
计算历程的几个阶段
在我国，首先是由数学家刘徽得出较精确的 圆周率。他采用割圆术。

3
5
7
2 n 1
并利用这个公式计算到了圆周率的100位.
1914年，印度数学家Srinivasa Ramanujan 表了下面的公式：
9801 ( 4 n )! ( 1103 26390 n ) 22 4 n 4 4 n 4 ( n ! ) 99 n 0
发
在1985年，Gosper用这个公式计算到了圆周率 的17500000位．
在中国
• 祖冲之: 在刘徽研究的基础上，进一步地发展， 经过既漫长又烦琐的计算，一直算到圆内接正 24576边形，而得到一个结论： • 3.1415926 ＜ π ＜ 3.1415927 同时得到π 的两个近似分数：约率为22／7； 密率为355／113。
• 他算出的 π 的8位可靠数字，不但在当时是最精 密的圆周率，而且保持世界记录九百多年。以致 于有数学史家提议将这一结果命名为“祖率”。
1 .圆周率 的幂级数计算方法
例1
利用 arctan x 的 Maclaurin 展开式 ， 计算 的近似值 .
解
1 2 4 n 12 n 1 x x (1 ) x . 2 1 x
xx ( 1 )x arctan x x . 35 2 n 1
1989年，David 和 Gregory Chudnovsky 发表 了下面的公式
n 1 ( 1 ) ( 6 n )! 13591409 5451401 n 12 , 3 3 3 n n ( 3 n )! ( n ! ) 2 0 640320
并在1994年计算到了4044000000位．它的另一 种形式是
在中国
• 刘徽：公元263年前后，刘徽提出著名的 “割圆 术”求出了比较精确的圆周率。他发现：当圆内 接正多边形的边数不断增加后，多边形的周长会 越来越逼近圆周长，而多边形的面积也会越来越 逼近圆面积。于是，刘徽利用正多边形面积和圆 面积之间的关系，从正六边形开始，逐步把边数 加倍：正十二边形、正二十四边形，正四十八边 形……，一直到正三○七二边形，算出圆周率等 于三点一四一六，将圆周率的精度提高到小数点 后第四位。

编写M 文件,用公式 +-+-≈7
1513114π
求π的近似值,计算到第n 项,n 由用户输入 【代码】：n=input('n=');
sum=0;
for i=1:n;
sum=sum+(-1)^(i+1)/(2*i-1);
end
※ input 有两种用法（1）. 在命令窗口中输入Val=input('请输入一个整数')，这样在命
令窗口中便会显示“请输入一个整数”提示用户进行输入操作，当用户输入一个整数后，便会被赋给Val 。

（2）. 当你在命令窗口输入str=input('Please input', 's')然后
从键盘输入：1+2+3，这样str 实际得到的是：'1+2+3'而不是6。

值得注意的是：需要注意的是，如果执行本函数时，用户敲了回车而不是输入了一个数，则该函数返回一个空矩阵。

可以用MATLAB 中的isempty 函数判断输入的是否为空。

※ for i=a:b; end MATLAB 中的循环语句，i=a 开始第一次计算从for 的下一行计算到end
截止，之后开始计算a+1，a+2…….b
※ a=a+b 在MATLAB 中表示将a+b 的值赋予a 通常在循环语句中出现
※ a=b 表示将b 的值赋给a
※ a= =b 表示判断a b 的值是否相等。

MATLAB已经发展成为多学科、多种工作平台 的功能强大的大型软件;
MATLAB已经成为线性代数、自动控制理论、 概率论及数理统计、数字信号处理、时间序列 分析、动态系统仿真等高级课程的基本教学工 具，是攻读学位的大学生、硕士生、博士生必 须掌握的基本技能.
数学实验matlab
❖语言简洁紧凑，使用方便灵活。MATLAB的基本数据单元是既不需要指
主窗口除了嵌入一些子窗口外，还主要包括 菜单栏和工具栏。
1．菜单栏
在MATLAB 6.5主窗口的菜单栏，共包含File、 Edit、View、Web、Window和Help 6个菜单项。
(1) File菜单项：File菜单项实现有关文件的操作。
(2) Edit菜单项：Edit菜单项用于命令窗口的编辑操作。
(3) View菜单项：View菜单项用于设置MATLAB集成环 境的显示方式。 (4) Web菜单项：Web菜单项用于设置MATLAB的Web 操作。 (5) Window菜单项：主窗口菜单栏上的Window菜单， 只包含一个子菜单Close all，用于关闭所有打开的编辑 器窗口，包括M-file、Figure、Model和GUI窗口。
数学实验matlab
Matlab简介与入门
数学实验matlab
什么是数学实验 ？
简单讲就是利用计算机和数学软件平台. 一方面，对学习知识过程中的某些问题 进行实验探究、发现规律； 另一方面，结合已掌握的数学（微积分、 代数与几何等）知识，去探究、解决一些简单 实际问题，从而熟悉从数学建模、解法研究到 实验分析的科学研究的方法。
在MATLAB里，有很多的控制键和方向键可用于命令行 的编辑。
工作空间窗口
工作空间是MATLAB用于存储各种变量和结 果的内存空间。在该窗口中显示工作空间中所有 变量的名称、大小、字节数和变量类型说明，可 对变量进行观察、编辑、保存和删除。

实验报告（五）完成人：L.W.Yohann注：本次实验主要学习了用MATLAB求解开普勒方程和方程求根的问题，了解学习了用fzero命令、二分法、Newton迭代法、一般迭代法求解方程，以及学习了非线性方程组的求解问题，完成后，小组对第90页的上级练习题进行了程序编辑和运行。

1.绘图并观察函数零点的分布.解：在编辑窗口输入：f=inline('x-0.5*sin(x)-1');fplot(f,[0,1])grid存盘后运行得2. 利用fzero 命令求解方程.解：在编辑窗口输入：f=inline('x-0.5*sin(x)-1');c=fzero(f,[1,2])存盘后运行得c =1.49873. 用二分法求解方程.求解（1）方程x^2-2=0在（0，2）内的近似根；（2）圆x^2+y^2=2与曲线y=e^-x 的两个交点；（3）方程∫t 21+t 2x 0dt =12的近似根. （1）解：在编辑窗口输入：00.10.20.30.40.50.60.70.80.91-1-0.9-0.8-0.7-0.6-0.5-0.4-0.3f=inline('x^2-2');x1=0;x2=2;while abs(x1-x2)>10^(-5)x3=(x1+x2)/2;if f(x3)==0break;elseif f(x1)*f(x3)>0x1=x3;else f(x2)*f(x3)>0;x2=x3;endendx0=x3存盘后运行得x0 =1.4142（2）解：在编辑窗口输入：f=inline('(x^2)*exp(2*x)+1-2*exp(2*x)');x1=0;x2=2;x5=-2;x6=0;while abs(x1-x2)>10^(-5)x3=(x1+x2)/2;if f(x3)==0break;elseif f(x1)*f(x3)>0x1=x3;else f(x2)*f(x3)>0;x2=x3;endendwhile abs(x5-x6)>10^(-5)x7=(x5+x6)/2;if f(x7)==0break;elseif f(x5)*f(x7)>0x5=x7;else f(x6)*f(x7)>0;x6=x7;endendx0=x3x4=x7存盘后运行得x0 =1.3922x4 =-0.3203（3）解：在编辑窗口输入：clear;clc;syms t xf1=(t^2)/(1+t^2);f2=int(f1,t,0,x);%¼ÆËã²»¶¨»ý·Öf=inline('x - atan(x)-0.5');x1=-5;x2=5;while abs(x1-x2)>10^(-5)x3=(x1+x2)/2;if f(x3)==0break;elseif f(x1)*f(x3)>0x1=x3;else f(x2)*f(x3)>0;x2=x3;endendx0=x3存盘后运行得x0 =1.47504.用Newton迭代法求解方程求解：x=0.5sinx+1的近似根；解：在编辑窗口输入：f=inline('x-0.5*sin(x)-1');df=inline('1-0.5*cos(x)');d2f=inline('0.5*sin(x)');a=1;b=2;dlt=1.0e-5;if f(a)*d2f(a)>0x0=a;elsex0=b;endm=min(abs(df(a)),abs(df(b)));k=0;while abs(f(x0))>m*dltk=k+1;x1=x0-f(x0)/df(x0);x0=x1;fprintf('k=%d x=%.5f\n',k,x0); end存盘后运行得k=1 x=1.54858k=2 x=1.49933k=3 x=1.498705.求解非线性方程组.试求非线性方程组{2x12−x1x2−5x1+1=0x1+3lgx1−x22=0的解，初值如下：（1）x0=[1.4,−1.5]（2）x0=[3.7,2.7]解：在编辑窗口输入：function f=group5(x)f=[2*x(1)^2-x(1)*x(2)-5*x(1)+1;x(1)+3*log10(x(1))-x(2)^2];(1)：输入：[f,fval]=fsolve('group2',[1.4,-1.5]) 运行得f =1.4589 -1.3968fval =1.0e-011 *0.0759-0.6178(2)：输入：[f,fval]=fsolve('group2',[3.7,2.7])运行得f =3.4874 2.2616fval =1.0e-006 *0.0059-0.20126.解决实际问题.为了在海岛I与某城市C之间铺设一条地下光缆，每千米光缆铺设成本在水下部分使C1万元，在地下部分使C2万元，为使得该光缆的总成本最低，光缆的转折点P（海岸线上）应该取在何处？如果实际测得海岛I与城市C之间的水平距离l=30km,海岛距海岸线垂直距离h1=15km，城市距海岸线垂直距离h=10km,C1=3000万元/km,C2=1500万元/km,求P点的坐标（误差<10−3km）.解：在编辑窗口输入：f=inline('(3000*x)/(x^2 + 225)^(1/2) + (750*(2*x - 60))/((x - 30)^2 + 100)^(1/2)'); x1=5;x2=10;while abs(x1-x2)>10^(-3)x3=(x1+x2)/2;if f(x3)==0break;elseif f(x1)*f(x3)>0x1=x3;else f(x2)*f(x3)>0;x2=x3;endendfprintf('x=%.5f',x3) 存盘后运行得x=7.69104>>。

matlab圆周率计算方法及实现Matlab是一种强大的数学软件，可以用它来进行各种复杂的计算和数据分析。

其中，计算圆周率也是Matlab中常见的问题之一。

本文将介绍几种常见的计算圆周率的方法，并提供相应的Matlab代码实现。

一、Leibniz公式Leibniz公式是一种最简单的计算圆周率的方法。

该公式由德国数学家莱布尼兹在17世纪提出，其基本思想是利用无穷级数逼近圆周率。

具体公式如下：$\pi=4\sum_{n=0}^{\infty}\frac{(-1)^{n}}{2n+1}$根据该公式，我们可以编写如下代码：function pi = leibniz(n)pi = 0;for i = 0:npi = pi + (-1)^i/(2*i+1);endpi = 4*pi;end其中，n为迭代次数。

二、Monte Carlo方法Monte Carlo方法是一种随机模拟方法，可以用于估计各种复杂问题的解。

在计算圆周率时，我们可以利用Monte Carlo方法来模拟投掷点落在圆内或者圆外的概率，并通过比值来估计圆周率。

具体实现过程如下：- 在一个正方形内随机生成n个点；- 统计落在圆内的点的个数m；- 通过比值$\frac{m}{n}$来估计圆周率，即$\pi\approx4\frac{m}{n}$。

根据该方法，我们可以编写如下代码：function pi = monte_carlo(n)x = rand(1,n);y = rand(1,n);d = sqrt(x.^2 + y.^2);m = sum(d <= 1);pi = 4*m/n;end其中，n为生成点的个数。

三、Bailey-Borwein-Plouffe公式Bailey-Borwein-Plouffe公式是一种快速计算圆周率的方法。

该公式由Jonathan M. Borwein、Peter B. Borwein和Simon Plouffe在1995年提出，其基本思想是利用二进制数位上的特殊性质来逼近圆周率。

本次实验注意：《实验五MALTAB基础知识（简单）》《实验五基于Matlab的信号频谱分析（复杂）》选作一个即可实验五MALTAB基础知识（一）实验目的 (2)（二）实验设备 (2)（三）实验要求 (2)（四）实验内容 (2)1.1 MATLAB基础知识 (2)1.1.1 MATLAB程序设计语言简介 (2)1.1.2 MA TLAB界面及帮助 (2)1.2 MA TLAB基本运算 (4)1.2.1 MA TLAB内部特殊变量和常数 (4)1.2.2 变量类型 (4)1.2.3 内存变量管理 (5)1.2.4 MA TLAB常用数学函数 (5)1.2.5 MA TLAB矩阵生成 (5)1.2.6 MA TLAB矩阵运算 (8)1.2.7 MA TLAB中的矩阵分析 (10)1.3 MA TLAB程序设计 (10)1.3.1 M文件 (10)1.3.2 程序控制结构 (12)实验五MALTAB基础知识（一）实验目的●了解MA TLAB 程序设计语言的基本特点，熟悉MA TLAB软件运行环境●掌握创建、保存、打开m文件及函数的方法●掌握变量等有关概念，具备初步的将一般数学问题转化为对应的计算机模型并进行处理的能力（二）实验设备计算机，Matlab软件（三）实验要求本实验属于验证实验，请根据（四）实验内容的步骤，运行相应的指令或例子，并将仿真结果截图至文档（请自己新建一个word文档，注意，并不一定所有指令或例子的实验结果都要截图，截图数目大于等于5个即可，自己选择性截图，答案不唯一，自由发挥）请在页眉处填写班级、学号、姓名，并将实验报告命名为“实验五_学号_姓名”，并通过FTP上传至指定文件夹。

（四）实验内容1.1 MATLAB基础知识1.1.1 MATLAB程序设计语言简介MA TLAB，Matrix Laboratory的缩写，是由MathWorks公司开发的一套用于科学工程计算的可视化高性能语言，具有强大的矩阵运算能力。

乐经良
Matlab计算
创建m文件
function y=calpi4(k) m=0; for n=1:k if rand(1)^2+rand(1)^2<=1 m=m+1; end; end; y=4*m/k;
观察结果，思考为什么？
乐经良
任务4
1) 用Monte Carlo 法计算，除了加大随机数， 在随机数一定时可重复算若干次后求平均值， 看能否求得5位精确数字？ 2) 设计方案用计算机模拟Buffon实验
1
n
n
1 2!
还有
1 e n ! ( n 1) !
n
1 k e lim(1 2 ) n n n
试用上述公式或其他方法近似计算e
乐经良
谢谢各位！
2
2
乐经良
1 an 2 2
4
于是 的值
6 2n 1 Sn 1 3 2n an
(刘徽计算到96边形面积，得到 3.141)
用Matlab计算
m文件
function calpi(n) a(1)=1; for i=1:n-1 a(i+1)=sqrt(2-sqrt(4-a(i)^2)); end S=3*2^(n-1)*a(n)
乐经良
其他方法
•
1/ 的展开式 Ramanujan 公式 (4n)!(1103 26390n) (n!)23964n 9801 n1 1

2 2•Βιβλιοθήκη 算术几何平均值迭代法1 an bn a0 1, b0 , an1 , bn1 anbn , M lim an limbn n n 2 2
乐经良
计算结果 >> calpi2_1(10) ans = 3.14159257960635063255949717131 >> calpi2_1(20) ans = 3.14159265358975625659354591335 >> calpi2_1(50) ans = 3.14159265358979323846264338328

• 在二维情形用matlab计算pi（版本3）：一次性产生随机矩阵，用for 循环做选择
基于蒙特卡罗方法用Mathematica和Matlab计算圆周率pi 基于蒙特卡罗方法用Mathematica和Matlab计算圆周率pi 在二维情形用matlab计算pi（版本1）：逐个产生随机点 原理：对于3个在[0,1]区间均匀分布的随机数，它们的平方和小于1的概率是pi/6. 原理：对于3个在[0,1]区间均匀分布的随机数，它们的平方和小于1的概率是pi/6. 在二维情形用matlab计算pi（版本2）：一次性产生随机矩阵，用while循环做选择 使用Mathematica软件： 使用Mathematica软件： 算法：生成N组在[0,1]区间均匀分布的2维随机数向量(x_i, y_i)，其中i=1,2,…,N，对其中满足平方和小于1的向量进行计数，设数目为N_A ，则N_A与N的比值为pi/4，也即是说，pi=4N_A/N. 数学原理：对于两个在[0,1]区间均匀分布的随机数x，y，它们的平方和x^2+y^2<1的概率是pi/4. 基于蒙特卡罗方法用Mathematica和Matlab计算圆周率pi
• 使用Mathematica软件： 用RandomReal生成随机数
用Select选取满足内积小于或等于1的行向量 用Dot求向量的内积
基于蒙特卡罗方法用Mathematica和Matlab计算圆周率pi 在二维情形用matlab计算pi（版本3）：一次性产生随机矩阵，用for循环做选择 使用Mathematica软件： 原理：对于4个在[0,1]区间均匀分布的随机数，它们的平方和小于1的概率是pi^2/32. 数学原理：对于两个在[0,1]区间均匀分布的随机数x，y，它们的平方和x^2+y^2<1的概率是pi/4. 使用Mathematica软件： 数学原理：对于两个在[0,1]区间均匀分布的随机数x，y，它们的平方和x^2+y^2<1的概率是pi/4. 基于蒙特卡罗方法用Mathematica和Matlab计算圆周率pi 基于蒙特卡罗方法用Mathematica和Matlab计算圆周率pi 在三维情形用matlab计算pi（版本2）：一次性产生随机矩阵，用while循环做选择 在二维情形用matlab计算pi（版本3）：一次性产生随机矩阵，用for循环做选择 在三维情形用matlab计算pi（版本2）：一次性产生随机矩阵，用while循环做选择 在三维情形用matlab计算pi（版本2）：一次性产生随机矩阵，用while循环做选择

MATLAB 源码求圆周率pi 阿基米德法%% The Computation of Pi by Archimedes% Copyright 2010, Bill McKeeman, Dartmouth College%% Abstract% It is famously known that Archimedes approximated $\pi$ by computing the % perimeters of many-edged regular polygons, one polygon inside the circle % and one outside. This presentation recapitulates and explains % Archimedes' computation. The surprise to me was how many "tweaks" % Archimedes applied at various stages of an otherwise systematic % approach.%%format longhexpi = imread('diagrams/HexPi.jpg');h=image(hexpi); axis off;%% The General Plan% The circumference of a circle of radius $R$ is $2 \pi R$ which defines % $\pi$. The perimeter of the inner hexagon is $6R$ which implies $6R < 2 % \pi R$, thus providing a lower bound $3 < \pi$. Computing the perimeter % of the outer hexagon supplies an upper bound. Increasing the number of % edges tightens the bounds. Thus thought Archimedes.%% The Context% The year was earlier than 200 BCE. Archimedes, who lived in Sicily, knew % of the work in Alexandria, where Euclid had published his "Elements" some % years before.%%% Greek mathematics used a kind of decimal integer. For example % ${_\prime}\alpha\nu\kappa\delta{^\prime}$ means 1424. Mathematicians % could express rational numbers, but did not have the concept of thedigit % 0, the decimal point or the means to represent irrational numbers such as % $\pi$ and square roots. Irrational numbers were approximated by rational % upper and lower bounds.%%% Finally, algebra was unavailable, so words were used to describe % computations and algorithms. The flavor of ancient presentation is % provided by Archimedes' statement "The surface of any sphere is four % times its greatest circle." We would write $A = 4 \pi R^2$. %% The Computation of $\pi$% The result presented by Archimedes translates to%%% $3 \frac{1}{7} > \pi > 3 \frac{10}{71}$%%% which appears, using the numbers of Archimedes' time, as%%% $\gamma^{\prime}\;\zeta^{\prime\prime} > \pi >% \gamma^{\prime}\;\iota^{\prime}\;o\alpha^{\prime\prime}$%% Advice to the Reader% What follows is a detailed rendition of the steps taken by Archimedes. % While following the logic, it may help the reader to look ahead to the % Summary to get the big picture of all the computations on one page. %% The Outer n-gonouter = imread('diagrams/archimedesOuter.jpg');h=image(outer); axis off;%%% Given the edge length $e_n$ of any regular n-gon enclosing a circle, one % can compute the edge length $e_{2n}$ of the enclosing regular 2n-gon % based on the equations%%% $\begin{array}{rcl}% x_n/R &=& (e_n - e_{2n})/e_{2n} \\% (x_n+R)/R &=& e_n/e_{2n} \\% (x_n+R)/e_n &=& R/e_{2n} \\% x_n^2 &=& R^2+e_n^2% \end{array}$%%% The first equation is a consequence of similar triangles. The second and % third equations are algebraic manipulations of the first. The fourth % equation is an application of the Pythagorean Theorem.%%% The last two equations are combined into an iteration.%%% $\begin{array}{rcl}% R/e_{2n} &=& (x_n+R)/e_n \\% &=& (R + \sqrt{R^2+e_n^2})/e_n \\% &=& R/e_n + \sqrt{(R/e_n)^2+1}% \end{array}$%%% Let $a_n \stackrel{\rm def}{=} R/e_n$.% The result is a recursive formula for $a_{2n}$.%%% $a_{2n} = a_n + \sqrt{a_n^2 + 1}$%%% The value of $a_n$ more than doubles at each step, implying that each % edge $e_n$ is less than one half the length of its predecessor at each % step, which is consistent with one's intuitive understanding of the % process.%%% Archimedes carried out the computation up to the 96-gon, starting with % the 30,60,90 degree triangle bounded by sides%%% $R, x_6, e_6$%%% which has known dimensions. In MATLAB, for the unit circle, R = 1;alpha_6 = 2*pi/12; % 30 degrees, half angle of the hexagon x_6 =R/cos(alpha_6); % R*sqrt(4/3), hypotenusee_6 = x_6*sin(alpha_6); % R*sqrt(1/3), half edge of the hexagon fprintf('e_6 = %.6g*R = R*sqrt(1/3)\n', e_6);%%% Since we know the half-edge $e_6$ for $R=1$, trying the formula is % straightforward.a_6 = R/e_6;a_12 = a_6 + sqrt(a_6^2 + 1);a_24 = a_12 + sqrt(a_12^2 + 1);a_48 = a_24 + sqrt(a_24^2 + 1);a_96 = a_48 + sqrt(a_48^2 + 1);e_96 = R/a_96;%%% The values of $a_n$ are more than doubling; the values $n/a_n$ are slowly % decreasing and converging to $\pi$ as expectedfprintf('n ');fprintf('%7d ', 6, 12, 24, 48, 96);fprintf('\na_n ');fprintf('%.6g ', a_6, a_12, a_24, a_48, a_96);fprintf('\nn/a_n ');fprintf('%.6g ', 6/a_6, 12/a_12, 24/a_24, 48/a_48, 96/a_96);fprintf('... %.6g\n', pi);%%% Unfortunately for Archimedes, he did not have double precision floating % point so he, and therefore we, still have some work to do. %% Hand Calculation% It is simpler for hand calculation if the intermediate results are % improper fractions with reasonably large integer parts. This is achieved % by multiplying the iteration formula by an arbitrary integer constant % $c$, the value of which can be chosen later.%%% $\begin{array}{rcl}% c R/e_{2n} &=& c (x_n+R)/e_n \\ % &=& (c R + c\sqrt{R^2+e_n^2})/e_n \\ % &=& c R/e_n + \sqrt{(c R/e_n)^2+c^2} %\end{array}$%%% Let $a_n \stackrel{\rm def}{=} c R/e_n$ and% $b_n \stackrel{\rm def}{=} \sqrt{(c R/e_n)^2+c^2}$.% The result is new formulas for $b_n$ and $a_{2n}$.%%% $\begin{array}{rcl}% b_n &=& \sqrt{a_n^2 + c^2} \\% a_{2n} &=& a_n + b_n% \end{array}$%%% Two mutually recursive formulas are provided because (1) $b_n$ is % irrational and must be replaced by a rational bound at each step and(2) % the edge of the inner n-gon is expressed in terms of $b_n$, not $a_n$. %% MATLAB fractions% This presentation now switches to "Greek mode," that is, variable % precision integers (vpi) and fractions (fr), provided by John D'Errico % and Ben Petschel on the MathWorks File Exchange. Except for type % conversions putting values into the vpi and fr domains, their packages % are non-intrusive in the code that follows. Thank you to MathWorks, John % and Ben.%% Bounds for Irrational Numbers% Archimedes states without explanation%%% $\frac{265}{153} < \sqrt 3 < \frac{1351}{780}$%%% This statement and other similar ones are easily checked without taking % square roots in MATLAB (and by Archimedes' readers).assert(fr(265/153)^2 < 3 && 3 < fr(1351/780)^2);%% Square Roots of Rationals% Heath, in his book "The Works of Archimedes," devotes 19 pages to % speculations by famous mathematicians on how Archimedes took square % roots. At all times Archimedes had to make a tradeoff between the better % accuracy of larger denominators against the computational labor involved % dealing with larger integers. The speculation I prefer isthat % Archimedes used the known iteration for square root with a good starting % value to get a very accurate fraction, then picked among thecontinued % fraction convergents to get suitable values of lesser accuracy. %%% The following iteration for square root of $a$,%%% $x \leftarrow \frac{1}{2}(x + a/x)$%%% can be restated for rational radicand $a/b$ as%%% $\begin{array}{rcl}% x/y &\leftarrow& \frac{1}{2}(x/y+a y/b x) \\% & \leftarrow & (b x^2 + a y^2)/2 b x y% \end{array}$%%% It is not hard to pick a good-enough starting point for an improper % fraction. If the integer has a odd number of digits and leading digit % $D$, pick $S = 1,2\; \mbox{or}\; 3$ so that $S^2 <= D$. If the integer % part has an even number of digits and leading digits $DE$, pick $S = % 1,\ldots 9$ so that $S^2 <= DE$. Then form the initial value $x$ by % appending to digit $S$ a $0$ for every pair of unexamined digits. That % is, for $349450$, $DE=34, S=5, x=500, y=1$.%%% Carrying out four steps of the computation for input $3$ gives rational % approximations for the $\sqrt 3$.n = 4; % the number of steps to takea = vpi(3);b = vpi(1); % the radicand a/bx = 1; y = 1; % starting values for the rational resultfor i=1:n % all values are now vpit = b*x^2 + a*y^2;y = 2*b*x*y;g = gcd(t,y); % known as Euclid's algorithmx = t/g; % exact divisiony = y/g;fprintf('%.16f %6d/%d\n', double([fr(x)/fr(y), x, y]));endfprintf('...\n%.16f sqrt(%d/%d)\n', ...sqrt(double(a/b)), double(a), double(b));%% Getting the Continued Fraction% Taking the most accurate of the four approximations above, num = vpi(18817); den = 10864;cf = vpi([]); % cf coefficients are always integerswhile den > 0 % compute continued fractionif den == 1 % choose long form of cfif num == 1r = 1; % finish with digit 1den = 0;else % force one more digitr = num - 1;num = 1;endelse % grab a digitr = num/den; % integer divisiont = num - r*den; % new num/dennum = den;den = t;endcf(end+1) = r;end%%% gives the continued fraction for $\frac{18817}{10864}$.fprintf('%d ', double(cf));%%% Evaluating the continued fraction gives a series of convergents, % alternating greater or less than the input fraction. As can be seen in % the printout, they are good approximations to $\sqrt 3$. %% fprintf('x/y (x/y)^2+-err\n'); % table titlen = numel(cf);for i = 1:n % evaluate continued fractionif i == 1 % startup caseh(i) = cf(i);k(i) = 1;elseif i == 2 % startup caseh(i) = cf(i)*h(i-1) + 1;k(i) = cf(i)*k(i-1);elseif i > 2 % rest of cfh(i) = cf(i)*h(i-1) + h(i-2);k(i) = cf(i)*k(i-1) + k(i-2);end% express result as L+-r/K^2 (need double for fprintf)H=double(h(i)); K=double(k(i)); L=round(H^2/K^2); r=H^2-L*K^2;conv = [sprintf('%d/%d', H, K), ' '];if r == 0fprintf('%s %d\n', conv(1:16), L);elseif r < 0fprintf('%s %d-%d/%d^2\n', conv(1:16), L, abs(r), K);elsefprintf('%s %d+%d/%d^2\n', conv(1:16), L, r, K);endend%%% Perhaps it was from this list that Archimedes chose %%% $\frac{1351}{780} > \sqrt 3 > \frac{265}{153}$%%% Since there is not a unique choice, we assume that Archimedes picked % values that would lead to the elegant result for $\pi$. Mathematicians % contemporary to Archimedes mention that later Archimedes did get a more % accurate result.%%% In any case the objective here is $\pi$, not square roots, so from this % point on I shall just report the choices Archimedes made for square % roots.%% The Upper Bound%%% For the upper bound computation, Archimedes chose the initial values c = 153;a_6 = c*R/e_6; % 265.0038b_6 = c*2; % which is precisea_6 = floor(a_6); % choose rational lower boundfprintf('a_6 = %d, b_6 = %d, c = %d\n', a_6, b_6, c);%%% There are two values for $a_6$. The first result comes from the geometry % of the triangle; the second is a chosen rational approximation. In this % case the chosen $a_6$ is just a little low. Further choices for $a_n$ % will be consistently low. The consequence is that rational estimate for % $a_{96}$ is less than the fully accurate irrational value of $a_{96}$. % The value of $e_{96}$ is therefore greater than the fully accurate % half-edge length of the enclosing 96-gon. Finally $96\times e_{96}$ is % greater than the half-perimeter of the enclosing 96-gon which is itself % greater than $\pi$.%% Iterating% Following the reasoning above, each time $b_n$ is computed, it must be % replaced by a rational bound that is smaller than the computed irrational % value. The amount of the difference will be small if the values $b_n$ are % well chosen. Considering Archimedes chose them, we should expect the % best. Nevertheless, we will check. There are three values given by % Archimedes for $b_{12} \ldots b_{48}$.assert(b_6^2 <= (c*R/e_6)^2 + c^2);a_12 = a_6 + b_6; % 571%%% Archimedes chooses $\;b_{12} = 591\frac{1}{8}% < \sqrt{571^2+153^2}% = \sqrt{349450}% = 591.143\ldots$b_12 = fr(591+1/8); assert(b_12^2 < a_12^2+c^2);a_24 = a_12 + b_12; % 1162 1/8%%% Archimedes chooses $b_{24}=1172\frac{1}{8}% < \sqrt{1162\frac{1}{8}^2+153^2} % = \sqrt{1373943\frac{33}{64}} % = 1172.153\ldots$b_24 = fr(1172+1/8); assert(b_24^2 < a_24^2+c^2);a_48 = a_24 + b_24; % 2334 1/4%%% Archimedes chooses $b_{48}=2339\frac{1}{4}% < \sqrt{2334\frac{1}{4}^2+153^2} % = \sqrt{5472132\frac{1}{16}} %= 2339.259\ldots$b_48 = fr(2339+1/4); assert(b_48^2 < a_48^2+c^2);a_96 = a_48 + b_48; % 4673 1/2%%% Archimedes makes one more approximation, replacing the final (ugly) % computed upper bound by one that is less tight but more elegant. pi_hi = fr(3+1/7) % Chosen by Archimedes for publication pi_est1 = 96*c/a_96 % Computed by Archimedespi % MATLAB's best%%% Using MATLAB, we can check that the values are ordered as expected. assert(pi_hi > pi_est1 && pi_est1 > pi);%% The Lower Bound% There is another diagram (in Heath), and another set of geometric% equations for the lower bounds which leads to the same recursive % formulas. One should not be surprised. In both cases the n-gons are the % same; it is the scale that differs.%%% This time the chosen rational values of $a_n$ must be greater than the % true values, causing the computed half-perimeter $96\times e_{96}$ to be % less than the actual half-perimeter of the inner 96-gon. %% % The initial values area_6 = 1351; b_6 = 1560; c = 780; assert(a_6^2+c^2 > b_6^2); %% Iterating% The values chosen for $b_n$ therefore must also be greater thanthe % irrational true values. Repeating the iteration, and in this case also % removing some common factors:a_12 = a_6 + b_6; % 2911%%% Archimedes chooses $b_{12}=3013\frac{3}{4}% > \sqrt{2911^2+780^2}% = \sqrt{9082321}% = 3013.689\ldots$b_12 = fr(3013+3/4); assert(b_12^2 > a_12^2+c^2);a_24 = a_12 + b_12; % 5924 3/4%%% Archimedes removes the common factor $\frac{13}{4}$ from $a_{24}$ and % $c$.a_24 = a_24*(4/13); % 1823c = c*(4/13); % 240%%% Archimedes chooses $b_{24}=1838\frac{9}{11} % > \sqrt{1823^2+240^2} % = \sqrt{3380929} % = 1838.730\ldots$ b_24 = fr(1838+9/11);assert(b_24^2 > a_24^2+c^2); a_48 = a_24 + b_24; % 3661 9/11 %%% Archimedes removes the common factor $\frac{40}{11}$ a_48 =a_48*(11/40); % 1007c = c*(11/40); % 66%%% Archimedes chooses $b_{48}=1009\frac{1}{6} % > \sqrt{1007^2+66^2} % = \sqrt{1018405} % = 1009.161\ldots$ b_48 = fr(1009+1/6);assert(b_48^2 > a_48^2+c^2); a_96 = a_48 + b_48; % 2016 1/6 %%% Archimedes chooses $b_{96}=2017\frac{1}{4} % >\sqrt{2016\frac{1}{6}^2+66^2}% = \sqrt{4069284\frac{1}{36}} % = 2017.247\ldots$ b_96 =fr(2017+1/4); assert(b_96^2 > a_96^2+c^2);pi_lo = fr(3+10/71); % Chosen by Archimedes for publication pi_est2 = 96*c/b_96; % Computed by Archimedes %%% The five values aredisp(pi_hi);disp(pi_est1);disp(pi);disp(pi_est2);disp(pi_lo);%%% The difference between the published bounds is disp(pi_hi - pi_lo) %%% Check that the five values are in decreasing order assert(pi_lo < pi_est2 && pi_est2 < pi && pi < pi_est1 && pi_est1 < pi_hi);%% QED% Therefore, $3\frac{1}{7} > \pi > 3\frac{10}{71}$.%%% Archimedes might have written: $OE\Delta$%% Summary% Here are all of the intermediate results in a table. Most of the action % is in the choosing of $b_n$. Recall that the chosen values of $b_n$ for % the outer case must be less than the computed values and conversely % greater for the inner case.%%% $\begin{array}{lllcllrl}% & \rm outer &&&& \rm inner && \rm remark \\ % a & b & c &\mathbf{n} & a & b & c \\ % 265 & 306 & 153 & \mathbf{6} & 1351 & 1560 & 780 & \rm start \\ % 571 & \sqrt{571^2+153^2} & 153 & \mathbf{12} % & 2911 & \sqrt{2911^2+780^2}\makebox[2pt]{} & 780 \\ % &591.1430\ldots &&&& 3013.6889\ldots && \rm compute \\ % & 591\frac{1}{8} &&&& 3013\frac{3}{4} && \rm choose \\ % &&& \mathbf{24}% & 5924\frac{3}{4}\makebox[5pt]{} && 780% & \mbox{factor out}\;\frac{13}{4} \\ % 1162\frac{1}{8} &\sqrt{1162\frac{1}{8}^2+153^2} & 153 & % & 1823 & \sqrt{1823^2+240^2} & 240 \\ % & 1172.1534\ldots &&&& 1838.7303\ldots && \rm compute \\ % & 1172\frac{1}{8} &&&& 1838\frac{9}{11} && \rm choose % \end{array}$%% $\begin{array}{lllcllrl}% && & \mathbf{48}% & 3661\frac{9}{11} && 240 & \mbox{factor out}\;\frac{40}{11} \\ % 2334\frac{1}{4} & \sqrt{2334\frac{1}{4}^2+153^2} & 153 % && 1007 &\sqrt{1007^2+66^2} & 66 \\ % & 2339.2589\ldots &&&& 1009.1605\ldots && \rm compute \\ % & 2339\frac{1}{4} &&&& 1009\frac{1}{6} && \rm choose \\ % 4673\frac{1}{2} & (b\;\mbox{not needed}) & 153 & \mathbf{96} % &2016\frac{1}{6} & \sqrt{2016\frac{1}{6}^2+66^2} & 66 \\ % &&&&&2017.2467\ldots && \rm compute \\ % &&&&& 2017\frac{1}{4} && \rm choose % \end{array}$%% $\begin{array}{cc}% \rm outer\;result & \rm inner\;result \\ % 96\times 153/4673\frac{1}{2}=3\frac{1335}{9347}=3.1428265982\ldots & % 96\times 66 / 2017\frac{1}{4}=3 \frac{1137}{8069}=3.1409096542\ldots% \end{array}$%% References% My principal references are "The Works of Archimedes" (1897), "The Method % of Archimedes Recently Discovered by Heiberg" (1912) and "A History of % Greek Mathematics" (1921) all by Sir Thomas Heath. The 1921 book is % still in print; Dartmouth Library dug wonderful dusty old tomes from % storage for the rest. I also watched an enjoyable lectureon Archimedes % in "Great Thinkers, Great Theorems" by Professor William Dunham (Great % Courses series) and read some entries in Wikipedia.。