求π的方法
- 格式：doc
- 大小：61.24 KB
- 文档页数：10
求π的方法
11级数教一班 高志磊 2011020231006
(一)古典方法
用圆内接正多边形和圆外切正多边形来逼近
源代码：
format long
disp('利用圆的正内接外切多边形的周长来逼近圆的周长');
sym n;
n=input('n= ');
for i=1:n
k=3.*2^i;%正多边形的边数
njb=2*sind (180./k);
zc1 = k.*njb;
yuanzhoulv1=zc1./2;
wqb=2*tand(180./k);
zc2=k.*wqb;
yuanzhoulv2=zc2./2;
i,k,yuanzhoulv1,yuanzhoulv2
end
运行结果：%当输入n=5时
利用圆的正内接外切多边形的周长来逼近圆的周长
n= 5
i =
1
k =
6 6边形 12边形 24边形
圆
yuanzhoulv1 =
3.00000000000000 yuanzhoulv2 =
3.46410161513775
i =
2
k =
12
yuanzhoulv1 =
3.10582854123025 yuanzhoulv2 =
3.21539030917347
i =
3
k =
24
yuanzhoulv1 =
3.13262861328124
yuanzhoulv2 =
3.15965994209750
i =
4
k =
48
yuanzhoulv1 =
3.13935020304687 yuanzhoulv2 =
3.14608621513143
i =
5
k =
96
yuanzhoulv1 =
3.14103195089051 yuanzhoulv2 =
3.14271459964537
（二）分析方法
用分析方法来求π的近似值，其中应用的主要工具是收敛的无穷乘积 和无穷级数。
1.利用walis 方法证明：
∏∞=⎪⎭⎫ ⎝⎛+⋅-⋅=⎪⎭⎫ ⎝⎛⋅⋅⎪⎭⎫ ⎝⎛⋅⋅⎪⎭⎫ ⎝⎛⋅⋅=112212227656543432122k k k k k
π
源代码：
function PI=wallis( )
k=input('n=');
PI=2;
for i=1:k
PI=PI*(2*i/(2*i-1)*2*i/(2*i+1));
end
运行结果：
输入n=20
ans =
3.10351696153923
输入n=10
ans =
3.06770380664350
2.利用泰勒展开求π
12)1(53arctan 1205
3+-=-+-=+∞=∑k x x x x x k k k 当x=1时，
121)1(5131140+-=-+-=∑∞=k k k
π
即
121)1(*4)51311(*40+-=-+-=∑∞=k k k
π
源代码：
function PI=taylor()
k=input ('k=');
sn=1;
for i=1:k
sn=sn+((-1)^i)*(1/(2*i+1));
end
PI=4*sn;
运行结果：
k=20
ans =
3.18918478227760
k=10
ans =
3.23231580940559
k=40
ans =
3.16597927284322
k=100
ans =
3.15149340107099
k=200
ans =
3.14656774718296
k=400
ans =
3.14408641529876
分析：觉得taylor 级数展开收敛性不太好，收敛速度很慢。
（三）概率方法
取一个二维数组（x,y），取一个充分大的正整数n，重复n次，每次独立地从（0，1）中随机地取一对数x和y ，分别检验x^2+y^2≤1是否成立。
设n 次试验中等式成立的共有m次，令π≈4m/n。
源代码：
n=input('输入想要重复的试验次数n=')
m=0;
for i=1:n
A=rand(1,2);
x=A(1,1);
y=A(1,2);
if (x^2+y^2<1)
m=m+1;
end
end
m
disp('所求的圆周率的值PI=')
4*m/n
format long
运行结果：
输入想要重复的试验次数n=100
n =
100
m =
77
所求的圆周率的值PI=
ans =
3.08000000000000
输入想要重复的试验次数n=200 n =
200
m =
159
所求的圆周率的值PI=
ans =
3.18000000000000
输入想要重复的试验次数n=300 n =
300
m =
231
所求的圆周率的值PI=
ans =
3.08000000000000
结果分析：随着试验次数的增大，pi 的值变小，比
3.14小。
(四）数值积分解法
利用公式
⎰-=1
0214dx x π，利用复合梯形公式。
源代码：
function bo = Trapezoid_M()
n=input('n=');
disp('Using complex of the trapezoid method.');
tn_tmp = 0;
a = 0;
b = 1;
h = (b - a) / n;
for k = 1:(n - 1)
xk = a + k * h;
tn_tmp = tn_tmp + f(xk);
end
tn = h / 2 * (f(a) + 2 * tn_tmp + f(b));
tn_tmp = 0;
for k = 1:n
tn_tmp = tn_tmp + f(a + k * h - h / 2);
end
t2n = (tn + h * tn_tmp) / 2;
bo = 4 * t2n;
function y = f(x)
y = (1 - x^2)^(1 / 2);
（其中Trapezoid_M( )函数中的n是等分区间的份数，因为是复化梯形，实际上是将区间等分成2n份）
运行结果如下：
n=10
Using complex of the trapezoid method
ans =
3.12846487975498
n=20
Using complex of the trapezoid method.
ans =
3.13694773421701
n=30
Using complex of the trapezoid method.
ans =
3.13906363461987
n=1000
Using complex of the trapezoid method.
ans =
3.14157950591196
n=2000
Using complex of the trapezoid method.
ans =
3.14158800514810
>>
结果分析：
此种方法收敛速度让然很慢，当n=1000，n=2000时，pi的精确度才到小数点后3位。
