补充实验

补充实验 圆周率的近似计算

vpa(pi,101) ans = 3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170680

１.圆周率π的计算历程 • 所谓“圆周率”是指一个圆的周长与其直径的比值。古今中外，许多人致力于圆周率的研究与计算。为了计算出圆周率的越来越好的近似值，一代代的数学家为这个神秘的数贡献了无数的时间与心血。 • 回顾历史，人类对 π 的认识过程，反映了数学和计算技术发展情形的一个侧面。 π 的研究，在一定程度上反映这个地区或时代的数学水平。德国数学家康托说：“历史上一个国家所算得的圆周率的准确程度，可以作为衡量这个国家当时数学发展水平的指标。” • 直到19世纪初，求圆周率的值应该说是数学中的头号难题。为求得圆周率的值，人类走过了漫长而曲折的道路。

实验时期 • 基于对一个圆的周长和直径的实际测量而得出的。 • 在古代世界，实际上长期使用 π ＝3这个数值。 • 最早见于文字记载的有基督教《圣经》中的章节，其上取圆周率为3。这一段描述的事大约发生在公元前950年前后。

几何法时期 真正使圆周率计算建立在科学的基础上，首先应归功于阿基米德。他是科学地研究这一常数的第一个人，是他首先提出了一种能够借助数学过程而不是通过测量的、能够把 π 的值精确到任意精度的方法。由此，开创了圆周率计算的第二阶段。　圆周长大于内接正多边形周长而小于外切正多边形周长．　据说阿基米德用到了正96边形才算出他的值域。

在中国 • 刘徽：公元263年前后，刘徽提出著名的 “割圆术”求出了比较精确的圆周率。他发现：当圆内接正多边形的边数不断增加后，多边形的周长会越来越逼近圆周长，而多边形的面积也会越来越逼近圆面积。于是，刘徽利用正多边形面积和圆面积之间的关系，从正六边形开始，逐步把边数加倍：正十二边形、正二十四边形，正四十八边形……，一直到正三○七二边形，算出圆周率等于三点一四一六，将圆周率的精度提高到小数点后第四位。

在中国 • 祖冲之:　在刘徽研究的基础上，进一步地发展，经过既漫长又烦琐的计算，一直算到圆内接正24576边形，而得到一个结论： • 3.1415926 ＜ π ＜ 3.1415927 同时得到π 的两个近似分数：约率为22／7；　　密率为355／113。 • 他算出的 π 的8位可靠数字，不但在当时是最精密的圆周率，而且保持世界记录九百多年。以致于有数学史家提议将这一结果命名为“祖率”。

分析法时期 • 这一时期人们开始摆脱求多边形周长的繁难计算，利用无穷级数或无穷连乘积来算 π 。 • 1593年，韦达（VieTa）给出这一不寻常的公式是 π 的最早分析表达式。甚至在今天，这个公式的优美也会令我们赞叹不已。它表明仅仅借助数字2，通过一系列的加、乘、除和开平方就可算出 π 值。

接着有多种表达式出现。如沃利斯1650年给出：接着有多种表达式出现。如沃利斯1650年给出：

1706年，英国天文学教授John Machin 利用 发现了下面的公式并利用这个公式计算到了圆周率的100位.

1914年，印度数学家Srinivasa Ramanujan　发表了下面的公式：在1985年，Gosper用这个公式计算到了圆周率的17500000位．

1989年，David 和 Gregory Chudnovsky 发表了下面的公式并在1994年计算到了4044000000位．它的另一种形式是

1995年，由David Bailey,Peter Borwein 和 Simon Plouffe 共同发表了下面的圆周率计算公式（简称BBP公式）该公式的最大优点在于：经后来人将该公式变形后打破了传统的计算方法，可以直接计算圆周率的任意第n位数，而不是先计算前面的n-1位数．

1997年，Fabrice Bellard发表了一个比BBP算法更快的公式从而，大大降低了圆周率近似值的计算量.

（1） 解

编写下面的程序: n=10; %选择展开式的次数 s=0; digits(22); %定义计算过程中的精度 for k=1:n s=s+4*(-1)^(k+1)/(2*k-1); end vpa(s,20) %定义显示精度为20位

（2） （3） Machin公式

练习

n=50; %定义等分积分区间数，可以更改 i=0:1/n:1; s=0; for k=1:length(i)-1 s=s+(1/(1+((i(k)+i(k+1))/2)^2))*1/n; end 4*s

表１

数值积分简介 在高等数学中有一类积不出的积分，如（概率积分）　　（椭圆积分）利用数值积分法

MATLAB命令 • 梯形方法 ——trapz(x,y) • 抛物线方法——quad(f,a,b) 如求积分的近似值梯形方法抛物线方法输入：x=0:0.1:1; y= exp(-x.^2); trapz(x,y) 输出：ans = 0.746211 y=inline('exp(-x.^2)'); quad(y,0,1) ans = 0.746826

梯形方法 条件：公式推导： ——trapz(x,y) 当区间划分为n等分时

复化梯形方法 当区间等距划分为n个子区间时现对每个子区间再二等分，得到2n个子区间

举例利用复化梯形算法求Pi的近似值. 取运用复化梯形算法

输入初值： • clear;a=0;b=1; • f=inline('4/(1+x*x)'); • t1=(b-a)/2*(f(a)+f(b)); • er=1;n=1; • whileer>1.0e-6 • h=(b-a)/n; • s=0; • fori=1:n • s=s+f(a+i*h-h/2); • end • t2=(t1+h*s)/2; • er=abs(t2-t1); • fprintf('t=%.6f,r=%.6f\n',t2,er); • n=2*n; t1=t2; • end No Yes 输出结果：STOP

抛物线方法 条件：分析：梯形法是对每个子区间用梯形面积近似曲线下面积累加而成；而抛物线法是对每个子区间考虑其中点，用三点决定的抛物线下面积来近似.

当区间划分为n(n>1)等分时 ——辛浦生(Simpson)方法

clear; a=0;b=1; f=inline('4/(1+x*x)'); t1=(b-a)/2*(f(a)+f(b)); er=1;s1=0; n=1; • whileer>1.0e-6 • h=(b-a)/n; • s=0; • fori=1:n • s=s+f(a+i*h-h/2); • end • t2=(t1+h*s)/2; • s=t2+(t2-t1)/3; • er=abs(s-s1); • fprintf('s=%.6f,r=%.6f\n',s,er); • s=s1; • n=2*n; • t1=t2; • end

s=3.133333,r=3.133333 s=3.141569,r=3.141569 s=3.141593,r=3.141593 s=3.141593,r=3.141593 s=3.141593,r=3.141593 s=3.141593,r=3.141593 s=3.141593,r=3.141593 s=3.141593,r=3.141593 s=3.141593,r=3.141593 s=3.141593,r=3.141593 s=3.141593,r=3.141593 s=3.141593,r=3.141593 s=3.141593,r=3.141593 s=3.141593,r=3.141593 s=3.141593,r=3.141593 s=3.141593,r=3.141593

cs=0 n=500 %随机取点数 for i=1:n a=rand(1,2); if a(1)^2+a(2)^2<=1 cs=cs+1 end end 4*cs/n １０１

3.18400000000000 3.10400000000000 3.13866666666667 3.12080000000000 3.14376000000000 　　从计算结果看：这种数据模拟算法收敛的速度很慢．虽然在实验次数较少时，算得的近似值距离真实值误差较大，但这种数据模拟方法简单易行．在精度要求不是很高的情况下，这种取随机数进行数据模拟的方法还是有一定的实用价值，通过简单编程我们可以模拟出许多数学现象．

补充实验

补充实验

Presentation Transcript