Matlab计算圆周率

n 1
练习
利用下面的等式计算 的值，并与例 1比较.
2 1 (a ) 2, 8 n1 ( 2n 1) 2 ( 1)n1 (b) , 2 12 n1 n 3 ( 1)n1 (c ) 3, 32 n1 ( 2n 1)
( 1) sin(2n 1) (d ) 3 . 8 n 1 ( 2n 1)
举例 利用复化梯形算法求Pi的近似值.
y
取
4
运用复化梯形算法
ba T1 ( f (a ) f (b)) 2
1
0
1
x
n 4 1 h f (a ih )) 0 1 x 2 dx T2n 2 (Tn h 2 i 1
输入初值：
n 1 h T2 n (Tn h f (a ih )) 2 2 i 1
在中国
• 祖冲之: 在刘徽研究的基础上，进一步地发展， 经过既漫长又烦琐的计算，一直算到圆内接正 24576边形，而得到一个结论： • 3.1415926 ＜ π ＜ 3.1415927 同时得到π 的两个近似分数：约率为22／7； 密率为355／113。
• 他算出的 π 的8位可靠数字，不但在当时是最精 密的圆周率，而且保持世界记录九百多年。以致 于有数学史家提议将这一结果命名为“祖率”。
在中国
• 刘徽：公元263年前后，刘徽提出著名的 “割圆 术”求出了比较精确的圆周率。他发现：当圆内 接正多边形的边数不断增加后，多边形的周长会 越来越逼近圆周长，而多边形的面积也会越来越 逼近圆面积。于是，刘徽利用正多边形面积和圆 面积之间的关系，从正六边形开始，逐步把边数 加倍：正十二边形、正二十四边形，正四十八边 形……，一直到正三○七二边形，算出圆周率等 于三点一四一六，将圆周率的精度提高到小数点 后第四位。
x2 x3
x4
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
分析法时期
• 这一时期人们开始摆脱求多边形周长的繁难 计算，利用无穷级数或无穷连乘积来算 π 。 • 1593年，韦达给出
2 2 2 2 2 2 2 2 2 2
这一不寻常的公式是 π 的最早分析表达式。甚至 在今天，这个公式的优美也会令我们赞叹不已。它 表明仅仅借助数字2，通过一系列的加、乘、除和 开平方就可算出 π 值。
接着有多种表达式出现。如沃利斯1650 年给出：

22446688 2 133 4 55 7 7
1706年，英国天文学教授John Machin 利用
x x x n1 x arctan x x 1 3 5 7 2n 1
发现了下面的公式
1 1 16 arctan 4 arctan , 5 239
No Yes
clear;a=0;b=1; f=inline('4/(1+x*x)'); t1=(b-a)/2*(f(a)+f(b)); er=1;n=1; while er>1.0e-6 h=(b-a)/n; s=0; for i=1:n s=s+f(a+i*h-h/2); end t2=(t1+h*s)/2; er=abs(t2-t1);
3
5
7
2 n1
并利用这个公式计算到了圆周率的100位.
1914年，印度数学家Srinivasa Ramanujan 表了下面的公式：
发
9801 (4n)! (1103 26390n) 2 2 4n 4 4n 4 ( n ! ) 99 n 0
在1985年，Gosper用这个公式计算到了圆周率 的17500000位．
1.圆周率的幂级数计算方法
例1
利用arctan x的Maclaurin展开式， 计算的近似值.
1 2 4 1 x x 2 1 x
3 5
解
(1) n 1 x 2 n
n 1
.
x x ( 1) x arctan x x 3 5 2n 1
2 n1
1989年，David 和 Gregory Chudnovsky 发表 了下面的公式
1 ( 1)n (6n)! 13591409 545140134n 12 , 3 3 3 n n 0 ( 3n)!( n! ) 640320 2
并在1994年计算到了4044000000位．它的另一 种形式是
数学实验五
圆周率π的近似计算
主 讲 :魏 平
１.圆周率π的计算历程
• 所谓“圆周率”是指一个圆的周长与其直径的比值。 古今中外，许多人致力于圆周率的研究与计算。为 了计算出圆周率的越来越好的近似值，一代代的数 学家为这个神秘的数贡献了无数的时间与心血。。 • 回顾历史，人类对 π 的认识过程，反映了数学和计 算技术发展情形的一个侧面。 π 的研究，在一定程 度上反映这个地区或时代的数学水平。德国数学家 康托说：“历史上一个国家所算得的圆周率的准确 程度，可以作为衡量这个国家当时数学发展水平的 指标。” • 直到19世纪初，求圆周率的值应该说是数学中的头 号难题。为求得圆周率的值，人类走过了漫长而曲 折的道路。
当区间划分为n(n>1)等分时
o
a
x1
x2 x3
Fra Baidu bibliotekx4
x5
b x

b a
n 1 n x k 1 x k h f ( x )dx S n ( f ( x0 ) 2 f ( xk ) f ( xn ) 4 f ( )) , 6 2 k 1 k 1 ba h , xk a kh k 0,1,2, , n n
将积分区间 [0， 1]分成n等份，在每一个小区间 上， 选取中点为 i , 编写下面的程序计算 的近似值.
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
实验时期
• 基于对一个圆的周长和直径的实际测量而 得出的。 • 在古代世界，实际上长期使用 π ＝3这个数 值。 • 最早见于文字记载的有基督教《圣经》中 的章节，其上取圆周率为3。这一段描述的 事大约发生在公元前950年前后。
几何法时期

真正使圆周率计算建立在 科学的基础上，首先应归 功于阿基米德。他是科学 地研究这一常数的第一个 人，是他首先提出了一种 能够借助数学过程而不是 通过测量的、能够把 π 的 圆周长大于内接正多边 值精确到任意精度的方法。 形周长而小于外切正多边 由此，开创了圆周率计算 形周长． 的第二阶段。 据说阿基米德用到了正 96边形才算出他的值域。
.
1 1 ( 1) 4 arctan 1 4(1 ). 3 5 2n 1
n 1
1 1 ( 1) 4 arctan 1 4(1 ). 3 5 2n 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位
o
a
x1s2
s3
s1
x2 x3
b x
y0
y1
y2
S 3
x3 x2 x x3 ( f ( x2 ) 4 f ( 2 ) f ( x3 )) 6 2
h/ 2
h/ 2
h S ( y0 4 y1 y2 ) 6
y
S 3 x x3 h ( f ( x2 ) 4 f ( 2 ) f ( x3 )) 6 2
426880 10005 . (6n)!(545140134 n 13591409) 3 3n ( n ! ) ( 3 n )! ( 640320 ) n 0
1995 年 ， 由 David Bailey,Peter Borwein 和 Simon Plouffe 共同发表了下面的圆周率计算公式 （简称BBP公式）
——trapz(x,y)

b a
ba f ( x )dx T1 ( f (a ) f (b)) 2
复化梯形方法
y
a
x1
当区间等距划分为n个子区间时

b a
n 1 h f ( x )dx Tn ( f (a ) 2 f ( xi ) f (bo )) , 2 i 1 ba h , xk a kh k 0,1,2, , n n
1 4 2 1 1 n . 8n 4 8n 5 8n 6 n 0 16 8n 1

该公式的最大优点在于：经后来人将该公式变 形后打破了传统的计算方法，可以直接计算圆周率 的任意第n位数，而不是先计算前面的n-1位数．
1997 年， Fabrice Bellard 发表了一个比 BBP 算 法更快的公式

1
0
e
x2
dx
抛物线方法 y=inline('exp(-x.^2)'); quad(y,0,1) ans = 0.746826
梯形方法
条件： f ( x)在[a, b]上连续或有有限个第一 类间断点 . 公式推导： b f ( x )dx s1 s2 sn
a
y
a
x1s2
1 1 32 1 256 64 n 0 1024 4n 1 4n 3 10n 1
n
64 4 4 1 . 10n 3 10n 5 10n 7 10n 9
从而，大大降低了圆周率近似值的计算量.
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

2
1 1 x
4
0
dx

（概率积分）
（椭圆积分）
抛物线法 利用数值积分法 梯形法
MATLAB命令 •梯形方法 ——trapz(x,y) •抛物线方法——quad(f,a,b) 如求积分的近似值 梯形方法 输入： x=0:0.1:1; y= exp(-x.^2); trapz(x,y) 输出： ans = 0.746211
fprintf('t=%.6f,r=%.6f\n',t2,er);
输出结果：STOP
n=2*n; t1=t2;
end
条件： 分析：
抛物线方法
y
s4
f ( x )在[a, b]上连续或有有限个第一 类间断点 .
梯形法是对每个子区间用梯形 面积近似曲线下面积累加而成； 而抛物线法是对每个子区间考虑 其中点，用三点决定的抛物线下 面积来近似.
s3
s4
f ( x i 1 ) f ( x i ) ( ( xi xi 1 )) o 2 i 1
n
s1
x2 x3
b x
当区间划分为n等分时
b
n 1 h f ( xi ) f (b)) , a f ( x )dx Tn 2 ( f (a ) 2 i 1 ba 其中 h , xk a kh k 1,2,, n 1 n