求π的方法
求π的方法
11级数教一班 高志磊 [1**********]06
(一)古典方法
用圆内接正多边形和圆外切正多边形来逼近
6边形 12边形
24边形
源代码:
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 圆
yuanzhoulv1 =
3.[1**********]000
yuanzhoulv2 =
3.[1**********]775
i =
2
k =
12
yuanzhoulv1 =
3.[1**********]025
yuanzhoulv2 =
3.[1**********]347
i =
3
k =
24
yuanzhoulv1 =
3.[1**********]124
yuanzhoulv2 =
3.[1**********]750
i =
4
k =
48
yuanzhoulv1 =
3.[1**********]687
yuanzhoulv2 =
3.[1**********]143
i =
5
k =
96
yuanzhoulv1 =
3.[1**********]051
yuanzhoulv2 =
3.[1**********]537
(二)分析方法
用分析方法来求π的近似值,其中应用的主要工具是收敛的无穷乘积 和无穷级数。
1.利用walis 方法证明:
2
2244661335572k2k22k12k1k1
源代码:
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.[1**********]923
输入n=10
ans =
3.[1**********]350
2.利用泰勒展开求π
2k1x3x5
kxarctaxnx(1)352k1 k0
当x=1时,
1111(1)k
352k1 k0 4
即
1114*(1)4*(1)k
352k1 k0
源代码:
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.[1**********]760
k=10
ans =
3.[1**********]559
k=40
ans =
3.[1**********]322
k=100
ans =
3.[1**********]099
k=200
ans =
3.[1**********]296
k=400
ans =
3.[1**********]876
分析:觉得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
m=m+1;
end
end
m
disp('所求的圆周率的值PI=')
4*m/n
format long
运行结果:
输入想要重复的试验次数n=100
n =
100
m =
77
所求的圆周率的值PI=
ans =
3.[1**********]000
输入想要重复的试验次数n=200
n =
200
m =
159
所求的圆周率的值PI=
ans =
3.[1**********]000
输入想要重复的试验次数n=300
n =
300
m =
231
所求的圆周率的值PI=
ans =
3.[1**********]000
结果分析:随着试验次数的增大,pi的值变小,比
(四)数值积分解法
利用公式
41
0x2dx,利用复合梯形公式。
源代码:
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)); 3.14小。
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.[1**********]498
n=20
Using complex of the trapezoid method.
ans =
3.[1**********]701
n=30
Using complex of the trapezoid method.
ans =
3.[1**********]987
n=1000
Using complex of the trapezoid method.
ans =
3.[1**********]196
n=2000
Using complex of the trapezoid method.
ans =
3.[1**********]810
>>
结果分析:
此种方法收敛速度让然很慢,当n=1000,n=2000时,pi的精确度才到小数点后3位。