数学建模 π的计算

实验报告 实验目的 1. 掌握用MATLAB软件求π值的，并对结果作初步分析； 2.掌握用MATLAB软件实现蒙特卡洛法的过程。 实验内容 运用数值积分法、泰勒级数法、模拟蒲丰实验和随机整数互素法 计算 π 值，取 n = （k=3, 4, …,7），所得结果和程序 运行时间，并作简要分析。 实验过程 一、 数值积分法 单位圆的面积等于以单位圆的圆心为原点建立直角坐标系， 单位圆在第一象限内的部分 G 是一个扇形，其面积 S = π/4 只需计算出 S 的近似值，它的4倍就是 π 的近似值 1、 梯形公式法 将扇形 G 分为 n 个宽度相等的部分 Gk（1 ≤ k ≤ n） Gk：曲边梯形 左右边界平行，上方边界为曲线 n→∞ Gk→梯 形 Gk 梯形面积： 等价于 π 取 n =10k （k=3, 4, …,7），按上述方法，通过计算扇形面积计 算 π 的近似值。 程序 Clear tic n=1e3(3, 4, …,7); x=linspace(0,1,n+1); y=(1-x.^2).^0.5; pai=(0.5*y(1)+0.5*y(n+1)+sum(y(2:n)))*4/n toc 2、辛普森（Simpson）公式法 G G 仍用分点 xi = a + i(b-a)/n (1≤ i ≤ n-1) 将区间[a, b]分成 n 等份， 直线 x = xi (1≤ i ≤ n-1) 将曲边梯形分成 n 个小曲边梯形，再作 每个小区间[xi-1, xi]的中点 将第 i 个小曲边梯形的上边界 y=f(x) (xi-1≤x≤xi) 近似地看作经 过这三点的抛物线段，则可求得： 其中 于 是 得 到 S ; 即 程序： clear tic n=1e3(3,4,5,6,7); x=linspace(0,1,n+1); x2=x(1:n)+1/(2*n); y=4./(1+x.^2); y2=4./(1+x2.^2); pai=((y(1)+y(n+1))+2*sum(y(2:n))+4*sum(y2))/(6*n) toc 二、 泰勒级数法 考虑反正切函数的泰勒级数 取 x =1，n =1e3(3,4,5,6,7) 程序： clear tic k=1e3(3,4,5,6,7); x=1; n=1:k; pai=sum((-1).^(n-1).*x.^(2.*n-1)./(2.*n-1)) toc 结果显示：花费时间长！准确度差！！ 原因：x =1 时得到的 arctan1 的展开式收敛太慢！ Maqin 公式法 提高收敛速度： x << 1 随着指数的增加，x 的幂快速接近于 0， 泰勒级数就会快速收敛 取 2k-1 = 63，其误差已经非常小，表明收敛速度非常快 有 再考虑收敛更快的： /4=arctan 1/2+arctan 1/3 程序： clear tic k=1e3(3,4,5,6,7); n=1:k; x=1/2; at_1=sum((-1).^(n-1).*x.^(2.*n-1)./(2.*n-1)); x=1/3; at_2=sum((-1).^(n-1).*x.^(2.*n-1)./(2.*n-1)); pai_1=4*at_1+4*at_2 toc /4=4arctan 1/5-arctan 1/239 程序： clear tic k=1e3(3,4,5,6,7); n=1:k; x=1/5; at_5=sum((-1).^(n-1).*x.^(2.*n-1)./(2.*n-1)); x=1/239; at_239=sum((-1).^(n-1).*x.^(2.*n-1)./(2.*n-1)); pai=16*at_5-4*at_239 toc 三、蒙特卡罗（Monte Carlo）法 在数值积分中，我们利用求 1/4 单位圆的面积来得到 π/4，进而 得到 π 1/4 单位圆是一个扇形 G，它是边长为 1 的单位正方形 G1 的一部分。G1 的面积 S1=1，只要能求出扇形 G 的面积 S 在正方形 G1 的面积 S1 中所占 的比例即可求得 π 值。 4 问题转化为求扇形面积在正方形面积中所占比例 k 在正方形中随机地投很多点，使所投的每 个点落在正方形中每一个位置的机会均等， 看其中有多少个点落在扇形内。 记投点总数为 n，落在扇形内 的点的个数为 m，则 k 在计算机上实现随机投点，需要能产生[0, 1]上均匀分布的随机数。 在 Matlab 中，使用 rand 函数实现这一功能。 rand 产生 1 个随机数 rand(n) 产生 n×n 的随机方阵 rand(m,n) 产生 m×n 的随机矩阵 需要产生两个[0, 1]上均匀分布的随机数 x，y，表示一个点的坐 标，这个点落在单位正方形内每个位置的机会均等，当 x2 + y2 ≤ 1 时，该点落在扇形内。 程序： clear tic n=1e7; x=rand(1,n); y=rand(1,n); c=x.^2+y.^2; m=sum(c<=1); k=m/n; pai=4*k toc 四、蒲丰（Buffon）掷针试验 1777 年法国数学家蒲丰提出了一种用于计算 π 的随机掷针试验， 步骤如下： （1）取一张白纸，在上面画许多间距为 d 的等距平行线； （2） 取一根长度为 l ( l < d ) 的均匀直针，随机地向画有平行线的纸 上掷去，一共投掷 n 次（n 是一个很大的整数），观察针和直线 相交的次数 m； （3）针和直线相交的概率 p = 2l / (πd )，取m/n 为 p 的近似值， 则 π 针和直线相交的概率为 p=2l/( d) 的证明过程 设 M 表示针的中点，x 表示针落下后 M 到最近一条平 行线的距离，表示针与直线的夹角。显然有 0≤x≤d/2 Ox 平面上表示一个矩形 Ω。针和某一条平 行线相交的充要条件是： Ox 平面上表示一个区域 A，点 M 落在区域 A 内的概率为 程序： clear tic n=1e7; d=2; l=1; x=abs(d*rand(1,n)-1); y=0.5*pi*rand(1,n); m=sum(x<=0.5*sin(y)); pai=n/m toc 五、随机整数互素法 蒙特卡洛法随机整数互素的概率： 取一个大的整数 N，在 1 到 N 之间随机地取一对整数 a，b，找 出它们的最大公约数(a, b)，当(a, b)=1 时称 a，b 互素。做 n 次 这样的实验，记录其中(a, b)=1 的情况出现的次数 m，算出 p = m/n 的值。注：随机整数互素的概率： 程序： Clear tic n=1e3(3,4,5,6,7); N=1e7; a=round(1+N*rand(1,n)); b=round(1+N*rand(1,n)); m=sum(gcd(a,b)==1); p=m/n pai=(6*n/m)^0.5 toc 实验结论 实验结果统计表 π=3.141592653589793 3.141555466911023 3.141591477611321 3.14159261640199 3.141592653552502 3.141592653552502 0.000345 s 0.001115 s 0.008649 s 0.113040 s 0.925687 s 3.141592653589790 3.141592653589802 3.141592653589803 3.141592653589796 3.141592653589856 0.000345 s 0.001375 s 0.010334 s 0.128711 s 0.931969 s 3.140592653839794 3.141492653590035 3.141582653589780 3.141591653589692 3.141592553589780 0.000409 s 0.001900 s 0.012589 s 0.113065 s 1.236871 s 3.141592653589792 3.141592653589792 3.141592653589792 3.141592653589792 3.141592653589792 0.001705 s 0.008424 s 0.075838 s 0.669452 s 5.293582 s 3.141592653589794 3.141592653589794 3.141592653589794 3.141592653589794 3.141592653589794 0.002017 s 0.009191 s 0.078560 s 0.694382 s 5.282906 s 3.152000000000000 3.142000000000000 3.149040000000000 3.141124000000000 3.142203600000000 0.000293 s 0.001201 s 0.011782 s 0.110932 s 1.068707 s 3.174603174603174 3.089280197713933 3.123048094940662 3.135533432625226 3.142013018616741 0.000337 s 0.001807 s 0.016089 s 0.150625 s 1.423036 s 3.216337604513385 3.144508978712735 3.144042691862183 3.140239850316813 3.141148071086062 0.039859 s 0.367654 s 3.461486 s 33.917678 s 239.278874 s 蒙特卡罗（M onte C arlo）法 模拟蒲松实验法 随机整数互素法 方法 n 数 值 积 分 法 泰 勒 级 数 法 4arctan 1/5-arctan 1/239 arctan 1/2+arctan 1/3 梯形公式 辛普森公式 普通算法 Maq in 公 式 结果分析： 由以上计算结果知：利用概率统计原理的三种方法(模拟蒲松实 验法 、蒙特卡罗（Monte Carlo）法、随机整数互素法)的计算 结果耗时长，运算次数多，精确度差，但该原理适用性强；而数 值积分法和普通的泰勒级数法的耗时短，运算次数较小，精确度 相对较高；泰勒级数法中的 Maqin 公式法耗时很短，运算次数小， 精确度高。 注：泰勒级数法中的 可将 算到 200 位。