《实验2π的计算》由会员分享，可在线阅读，更多相关《实验2π的计算（26页珍藏版）》请在金锄头文库上搜索。
1、. 22|21 2188 .224 4242)21(11 ()21(| . 22. 11121111221212221122NNnaSaOEABSSnaaaaaEFAFAEaSanSannnOAEnnnnnn，故又注意到，相应的边数，故，注意到，边长和面积，则步单位圆内接多边形的第分别是，设值也就越接近圆周率近单位圆的面积，其数边形的面积就越来越接那么内接正多正多边形的边数开始，逐步成倍地增加从单位圆的内接正方形实验实验 2 2： 的计算的计算学院：机械与动力工程学院 姓名：唐子彦 学号：5150209101371 1、实验目的实验目的通过求 的近似值，了解历史上计算 值的一些方法，包括刘徽割
2、圆术、级数展开、数值积分和 Monte Carlo 法等. 复习微积分中相关知识，比较它们的差异，了解计算方法对提高计算效率的意义.2 2、实验内容实验内容（1 1）利用鲁道夫割圆术计算利用鲁道夫割圆术计算 值值【问题】德国人鲁道夫一生计算圆周率，他同样是用圆的内接多边形逼近圆周. 不过，他是从正方形开始成倍增加边数. 试推导出他计算所采用的递推公式，然后求 的近似值到 10 位和 20 位有效数字.【解】图 1.1.1 ,21 , 21 421 ,211 12 nnanS nnana nnnnn且， 且，的递推公式：由此可得鲁道夫割圆术NN为通过“鲁道夫法”计算 的近似值，现编写 M 文件如
3、下：为求出指定有效数字位数的 值，还需编写 M 文件如下：运行结果如下：图 1.2【答】由图 1.2 不难发现， “鲁道夫法”计算 值效率较高，n=16 时即可得到 的 10 位有效数字，n=32 时即可得到 的 20 位有效数字. 然而，由于中途过程涉及开方等运算，因此计算较为复杂繁琐.function PI=m2_7(n) %该函数利用鲁道夫割圆术计算 的近似值，n 为计算次数a(1)=sym(sqrt(2); %为尽可能消除累计误差，此处采用符号运算 for i=1:n-1a(i+1)=sym(sqrt(2-sqrt(4-a(i)2);endS=sym(2n*a(n);PI=vpa(S,
4、30);function m2_8(m) %该函数用于计算指定有效数字位数的 的近似值，m 为有效数字位数 n=1;PI=vpa(pi,100);while abs(PI-m2_7(n)=10(-m+1)n=n+1;enddisp(迭代次数 n=,num2str(n)vpa(m2_7(n),m) %仅显示 的近似值的指定位有效数字，下同（2 2）利用幂级数展开式计算利用幂级数展开式计算 值值【问题】简单公式，Machin 公式，31arctan21arctan4 2391arctan51arctan44以及公式. 试验证上述三个公式（分别记为公式81arctan51arctan21arctan
5、41、2、3） ，并利用反正切函数的幂级数展开式求 值，比较上述三公式的计算效率. 此外，再找出一种利用幂级数展开式求 的方法并验证之.【解】.1424404tan131 21131 21tantan1tantan)tan(.31tan,21tan31arctan,21arctan) 1 (得证，公式，故且，因为，则记 .244203404tan12391 11912012391 119120tan4tan1tan4tan)4tan(1191204tan125)51(1512tan1tan22tan.2391tan,51tan2391arctan,51arctan)2(22得证，公式，故且，因
6、为，则再记 yxyxyxyxyxxxxxyxyx.34.431 tan1tan1)4(tan3181 51181 51tantan1tantan)tan(.81tan,51tan,21tan81arctan)3(得证，公式即，故且，因为，则再记zxzxzxzxzxzxz 为利用上述三个公式求出 值，现编写以下三个 M 文件：第一个 M 文件，对应于公式 1：第二个 M 文件，对应于公式 2：第三个 M 文件，对应于公式 3：为比较它们计算指定有效数字位数的 值的效率，还需编写 M 文件如下（详见第五页）：function PI=m2_10(n)%该函数利用简单公式：PI/4=arctan(1/
7、2)+arctan(1/3)的幂级数展开式计算 S=0;for i=0:n-1S=S+(-1)i/(2*i+1)*(1/(2(2*i+1)+1/(3(2*i+1);endPI=vpa(4*S,50);function PI=m2_11(n) %该函数利用 Machin 公式：PI/4=4*arctan(1/5)-arctan(1/239)的幂级数展开式计算 S=0;for i=0:n-1S=S+(-1)i/(2*i+1)*(4*(1/5)(2*i+1)-(1/239)(2*i+1);endPI=vpa(4*S,50);function PI=m2_12(n) %该函数利用公式：PI/4=arc
8、tan(1/2)+arctan(1/5)+arctan(1/8)的幂级数展开式计算 S=0;for i=0:n-1S=S+(-1)i/(2*i+1)*(1/2)(2*i+1)+(1/5)(2*i+1)+(1/8)(2*i+1);endPI=vpa(4*S,50);function m2_13(m)%该函数用于计算对指定有效数字位数的 ，利用不同幂级数展开式所需的项数 n1=1;n2=1;n3=1;PI=vpa(pi,100);while abs(PI-m2_10(n1)=10-(m-1)n1=n1+1;endwhile abs(PI-m2_11(n2)=10-(m-1)n2=n2+1;endw
9、hile abs(PI-m2_12(n3)=10-(m-1)n3=n3+1;enddisp(有效数字位数 m=,num2str(m)disp(公式一所需的项数 n1=,num2str(n1)disp(公式二所需的项数 n2=,num2str(n2)disp(公式三所需的项数 n3=,num2str(n3) disp( )运行结果如下（每五个为一组）：图 2.1 图 2.2 图 2.3 图 2.4 从上述四幅图中可以看出，公式 1、3 的计算效率基本相同，而公式 2 的计算效率高于其他两个. 然而，从有效数字位数 m=15 开始，所得结果中公式 1、2、3 所需项数不再增加，这可能是 MATLA
10、B 本身的原因. 个人猜测：当 MATLAB 计算到一个与 极为相近的数时，可能将其自动补全为 ，而没有继续计算. 为试图解决该问题，可改用 C+进行编程，所需的 CPP 文件如下：/利用 arctan 的 Taylor 级数展开式计算 的近似值#include #include #define PI 3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117068const long double ONE=1;using namespace std;lon
11、g double fun_1(int n);long double fun_2(int n);long double fun_3(int n); /下接第 7 页int main() /上接第 6 页int n1=0,n2=0,n3=0;int m;coutm;for(int i=1;i=epsilon)n1+;while(abs(PI-fun_2(n2)=epsilon)n2+;while(abs(PI-fun_3(n3)=epsilon)n3+;cout17 的情形更是如此（参见图 2.8）.return 4*S; /上接第 7 页long double fun_2(int n) /公式
12、2 利用 /4=4arctan(1/5)-arctan(1/239)long double S=0;for(int i=0;i=10(-m+1)S=S+a/(2*n+1)*b)*(1/2)(2*n+1);n=n+1;a=a*(2*n-1);b=b*(2*n);enddisp(为达到,num2str(m),位有效数字，所需项数为 n=,num2str(n) vpa(6*S,m)输入图 2.9 所示命令行，运行结果如下：图 2.10【答】对比以上四公式的计算结果不难发现，公式 1、3 计算效率大致相同且较低，公式 2 的计算效率最高，公式 4 的计算效率介于它们之间比公式 1、3 略高。(3)(3
13、)数值积分计算数值积分计算 值值【问题】利用数值积分计算 ，分别用“梯形法”和“Simpson 法”精确到 10位有效数字，再用“Simpson 法”精确到 15 位有效数字.【解】在本题中， “梯形法”和“Simpson 法”的计算原理如下：为利用数值积分计算 值，现编写如下两个 M 文件：function PI=m2_17(n)%该函数利用“梯形法”进行数值计算，从而求出 的近似值，n 为划分份数 x=0:1/n:1;y=1./(1+x.2);S=2*sum(y)-1-0.5;PI=vpa(2*S/n,50);function PI=m2_18(n) %该函数利用“Simpson 法”进行
14、数值计算，从而求出 的近似值，n 为划分份数 x=0:1/(2*n):1;y=1./(1+x.2);x1=1/n:1/n:(n-1)/n;y1=1./(1+x1.2);S1=sum(y1);S2=sum(y)-S1-1.5;S=1/(6*n)*(1.5+2*S1+4*S2);PI=vpa(4*S,50);.42)(32.42)(61.112 , 2 , 1 , 0 ,22 1 , 0)2().23(2. )23(21)(211)(21.11, 2 , 1 , 0 ,0,1) 1 (.1 , 0,11.41d 41111222011112220201000012210def2mkmkkkmmkmkkkmkkknkknknkknknkkkkkkyyyymyyyymAxymkmkxmSimpsonynynyyynnyyAxynknkxnxxyAAxx故计算，可以得到：曲线用抛物线代替进行将相邻两个区间的函数，则等分，令将区间法故算，可以得到：曲线用直线段代替来计将每个小区间上的函数，则等分，令将区间梯形法记，则由于为计算指定有效数字位数的 值，还需编写以下两个 M 文件：(1)梯形法：function m2_19(m) %该函数利用“梯形法”计算指定有
