Pi值的计算(mathematica数学实验报告)
数学实验(计算圆周率)
数学实验报告一、 实验问题人们很早就发现圆的周长和直径的比是个常数，也就是我们都熟悉的圆周率π，那么这个常数的值究竟是多少哪？我国的伟大数学家祖冲之将π精确到3.1415926-3.1415927之间，西方的数学研究者也从来没有停止过对π的计算，那么在今天计算机技术如此发达的情况下，我们如何将π的近似值精确到更多位数哪？二、问题的分析π的精确度问题困扰大家已久，所以我们不妨在MATLAB 中采用循环来实现对π的更高位数的精确。
不妨以精确到小数点后9位为例。
我们知道 31arctan 21arctan 1arctan 4+==π 因此可以利用）（π）（321-1-2n 1-2n 1n 1-n 111-2n 431arctan 21arctan 4+=⎪⎭⎫ ⎝⎛+=∑∞=三、实验内容在MATLAB 中键入以下程序clear;n=0;r=1;p=0;k=-1;a=1;b=1;while r>=1.0e-9n=n+1;k=k*(-1);a=4*a;b=9*b;pl=p+k/(2*n-1)*(2/a+3/b);r=abs(4*(pl-p));fprintf('n=%.0f,p=%.10f\n',n,4*pl);p=pl;end四、问题求解结果和结论n=1,p=3.3333333333n=2,p=3.1172839506n=3,p=3.1455761317n=4,p=3.1408505618n=5,p=3.1417411974n=6,p=3.1415615879n=7,p=3.1415993410n=8,p=3.1415911844n=9,p=3.1415929813n=10,p=3.1415925796n=11,p=3.1415926705n=12,p=3.1415926497n=13,p=3.1415926545n=14,p=3.1415926534n=15,p=3.1415926536对泰勒级数，当｜x ｜越小时级数收敛速度越快，这启示我们另外构造相关级数来逼近π。
Monte Carlo方法计算Pi
实验Monte Carlo方法计算Pi实验要求：以OpenMP实现Monte Carlo计算Pi的并行程序实验分析：通过蒙特卡罗算法计算圆周率的主导思想是：统计学（概率）一个正方形有一个内切圆，向这个正方形内随机的画点，则点落在圆内的概论为P=圆面积/正方形面积。
1. 在一个平面直角坐标系下，在点（1,1）处画一个半径为R=1的圆,以这个圆画一个外接正方形，其边长为R=1(R=1时，圆面积即Pi)。
2. 随机取一点（X，Y）使得0<=X<=2R并且0<=Y<=2R，即随机点在正方形内。
3. 判断点是否在圆内，通过公式(X-R)(X-R)+(Y-R)(Y-R)<R*R计算。
4. 设所有点的个数为N，落在圆内的点的个数为M，则P=M/N=4*R*R/Pi*R*R=4/PiPi=4*N/M当实验次数越多（N越大），所计算出的Pi也越准确。
但计算机上的随机数毕竟是伪随机数，当取值超过一定值，也会出现不随机现象，因为伪随机数是周期函数。
如果想提高精度，最好能用真正的随机数生成器(需要更深的知识)。
运行结果采用并行采用串行代码：#include"stdafx.h"#include<stdio.h>#include<omp.h>#include<stdlib.h>#include<time.h>#define NUM_THREADS 4int main(){long long max=10000000;long long i,count=0;double x,y,bulk,starttime,endtime;time_t t;starttime=clock();// 产生以当前时间开始的随机种子srand((unsigned) time(&t));omp_set_num_threads(NUM_THREADS);#pragma omp parallel for reduction(+:count) private(x,y)for(i=0;i<max;i++){x=rand(); x=x/32767;y=rand(); y=y/32767;if((x*x+y*y)<=1)count++;}bulk=4*(double(count)/max);endtime= clock();printf("pi is %f \n", bulk);printf("Running time is %f \n", endtime-starttime);return 0;}加红色代码之前为串行代码，加上红色代码之后为并行代码。
Pi值的计算(mathematica数学实验报告)
arctan x x x3 x5 (1)k1 x 2k1
35
2k 1
来计算 。 从反正切函数的泰勒级数，进行如下编程来计算 ，实验运行如下：
从实验过程可以看出，这种方法花费的时间很长。原因是当 x=1 时得到的 arctan1的
展开式收敛太慢。要使泰勒级数收敛得快，容易想到，应当使 x 的绝对值小于 1，最好
实验基本理论和方法：
1、Mathematica中常用绘图函数Plot在绘制高次函数时的方法；
2、计算圆周率 的数值积分法、泰勒级数法、蒙特卡罗法，并且利用特定的公式来
计算圆周率 。
实验内容和步骤：
（1）数值积分法计算
半径为 1 的圆称为单位圆，它的面积等于 。只要计算出单位圆的面积，就算出了 。 在坐标轴上画出以圆点为圆心，以 1 为半径的单位圆（如下图），则这个单位圆在第一 象限的部分是一个扇形，而且面积是单位圆的 1/4，于是，我们只要算出此扇形的面积， 便可以计算出 。
0
4
利用 Mathematics 编程计算上式，过程如下：
从而得到 的近似值为 3.14159265358979323846264338328，可以看出，用这种方法 计算所得到的 值是相当精确的。n 越大，计算出来的扇形面积的近似值就越接近 的 准确值。
（2）泰勒级数法计算 利用反正切函数的泰勒级数
只要计算出单位圆的面积就算出了为半径的单位圆如下图则这个单位圆在第一象限的部分是一个扇形而且面积是单位圆的14于是我们只要算出此扇形的面积便可以计算出在计算扇形面积时很容易想到使用数学分析中积分的方法第一象限中的扇形由曲线及两条坐标轴围成实际操作中我们不能准确地计算它的面积于是就通过分割的方法将其划分为许多小的梯形通过利用梯形的面积近似于扇形面积来计算利用mathematics编程计算上式过程如下
Mathematica 基本运算
Mathematica 基本运算a+mathematica数学实验(第2版)b+c 加a-b 减a b c 或a*b*c 乘a/b 除-a 负号a^b 次方Mathematica 数字的形式256 整数2.56 实数11/35 分数2+6I 复数常用的数学常数Pi 圆周率，π=3.141592654…E 尤拉常数，e=2.71828182…Degree 角度转换弧度的常数，Pi/180I 虚数，其值为√-1Infinity 无限大指定之前计算结果的方法% 前一个运算结果%% 前二个运算结果%%…%(n个%) 前n个运算结果%n 或Out[n] 前n个运算结果复数的运算指令a+bI 复数Conjugate[a+bI] 共轭复数Re[z], Im[z] 复数z的实数/虚数部分Abs[z] 复数z的大小或模数(Modulus)Arg[z] 复数z的幅角(Argument)Mathematica 输出的控制指令expr1; expr2; expr3 做数个运算，但只印出最后一个运算的结果expr1; expr2; expr3; 做数个运算，但都不印出结果expr; 做运算，但不印出结果常用数学函数Sin[x],Cos[x],Tan[x],Cot[x],Sec[x],Csc[x] 三角函数，其引数的单位为弧度Sinh[x],Cosh[x],Tanh[x],… 双曲函数ArcSin[x],ArcCos[x],ArcTan[x] 反三角函数ArcCot[x],ArcSec[x],ArcCsc[x]ArcS inh[x],ArcCosh[x],ArcTanh[x],… 反双曲函数Sqrt[x] 根号Exp[x] 指数Log[x] 自然对数Log[a,x] 以a为底的对数Abs[x] 绝对值Round[x] 最接近x的整数Floor[x] 小于或等于x的最大整数Ceiling[x] 大于或等于x的最小整数Mod[a,b] a/b所得的馀数n! 阶乘Random[] 0至1之间的随机数（最新版本已经不用这个函数，改为使用RandomReal[]）Max[a,b,c,...]，Min[a,b,c,…] a,b,c,…的极大/极小值数值设定x=a 将变数x的值设为ax=y=b 将变数x和y的值均设为bx=. 或Clear[x] 除去变数x所存的值变数使用的一些法则xy 中间没有空格，视为变数xyx y x乘上y3x 3乘上xx3 变数x3x^2y 为x^2 y次方运算子比乘法的运算子有较高的处理顺序四个处理代数指令Expand[expr] 将expr展开Factor[expr] 将expr因式分解Simplify[expr] 将expr化简成精简的式子FullSimplify[expr] Mathematica 会尝试更多的化简公式，将expr化成更精简的式子多项式/分式转换函数ExpandAll[expr] 把算式全部展开Together[expr] 将expr各项通分在并成一项Apart[expr] 把分式拆开成数项分式的和Apart[expr,var] 视var以外的变数为常数，将expr拆成数项的和Cancel[expr] 把分子和分母共同的因子消去分母/分子的运算Denominator[expr] 取出expr的分母Numerator[expr] 取出expr的分子ExpandDenominator[expr] 展开expr的分母ExpandNumerator[expr] 展开expr的分子多项式二种转换函数Collect[expr,x] 将expr表示成x的多项式，如Collect[expr,{x,y,…}] 将expr分别表示成x,y,…的多项式FactorTerms[expr] 将expr的数值因子提出，如4x+2=2(2x+1)FactorTerms[expr,x] 将expr中把所有不包含x项的因子提出FactorTerms[expr,{x,y,…}] 将expr中把所有不包含{x,y,...}项的因子提出函数和指数运算TrigExpand[expr] 将三角函数展开TrigFactor[expr] 将三角函数所组成的数学式因式分解TrigReduce[expr] 将相乘或次方的三角函数化成一次方的基本三角函数之组合ExpToTrig[expr] 将指数函数化成三角函数或双曲函数TrigToExp[expr] 将三角函数或双曲函数化成指数函数复数、次方乘积之展开ComplexExpand[expr] 假设所有的变数都是实数来对expr展开ComplexExpand[expr,{x,y,…}] 假设x,y,..等变数均为复数来对expr展开PowerExpand[expr] 将项次、系数与最高次方Coefficient[expr,form] 于expr中form的系数Exponent[expr,form] 于expr中form的最高次方Part[expr,n] 或expr[[n]] 在expr项中第n个项代换运算子expr/.x->value 将expr里所有的x均代换成valueexpr/.{x->value1,y->value2,…} 执行数个不同变数的代换expr/.{{x->value1},{x->value2},…} 将expr代入不同的x值expr//.{x->value1,y->value2,…} 重复代换到expr不再改变为止求解方程式的根Solve[lhs==rhs,x] 解方程式lhs==rhs，求xNsolve[lhs==rhs,x] 解方程式lhs==rhs的数值解Solve[{lhs1==rhs1,lhs2==rhs2,…},{x,y,…}] 解联立方程式，求x,y,…NSolve[{lhs1==rhs1,lhs2==rhs2,…},{x,y,…}] 解联立方程式的数值解FindRoot[lhs==rhs,{x,x0}] 由初始点x0求lhs==rhs的根四种括号(term) 圆括号，括号内的term先计算f[x] 方括号，内放函数的引数{x,y,z} 大括号或串列括号，内放串列的元素p[[i ]] 或Part[p,i] 双方括号，p的第i项元素p[[i,j]] 或Part[p,i,j] p的第i项第j个元素缩短输出指令expr//Short 显示一行的计算结果Short[expr,n] 显示n行的计算结果Command; 执行command，但不列出结果查询物件?Command 查询Command的语法及说明??Command 查询Command的语法和属性及选择项?Aaaa* 查询所有开头为Aaaa的物件函数定义、查询与清除f[x_]= expr 立即定义函数f[x]f[x_]:= expr 延迟定义函数f[x]f[x_,y_,…] 函数f有两个以上的引数?f 查询函数f的定义Clear[f] 或f=. 清除f的定义Remove[f] 将f自系统中清除掉含有预设值的Patterna_+b_. b的预设值为0，即若b从缺，则b以0代替x_ y_ y的预设值为1x_^y_ y的预设值为1条件式的自订函数lhs:=rhs/;condition 当condition成立时，lhs才会定义成rhsIf指令If[test,then,else] 若test为真，则回应then，否则回应elseIf[test,then,else,unknow] 同上，若test无法判定真或假时，则回应unknow 极限Limit[expr,x->c] 当x趋近c时，求expr的极限Limit[expr,x->c,Direction->1]Limit[expr,x->c,Direction->-1]微分D[f,x] 函数f对x作微分D[f,x1,x2,…] 函数f对x1,x2,…作微分D[f,{x,n}] 函数f对x微分n次D[f,x,NonConstants->{y,z,…}] 函数f对x作微分，将y,z,…视为x的函数全微分Dt[f] 全微分dfDt[f,x] 全微分Dt[f,x1,x2,…] 全微分Dt[f,x,Constants->{c1,c2,…}] 全微分，视c1,c2,…为常数不定积分Integrate[f,x] 不定积分∫f dx定积分Integrate[f,{x,xmin,xmax}] 定积分Integrate[f,{x,xmin,xmax},{y,ymin,ymax}] 定积分数列之和与积Sum[f,{i,imin,imax}] 求和Sum[f,{i,imin,imax,di}] 求数列和，引数i以di递增Sum[f,{i,imin,imax},{j,jmin,jmax}]Product[f,{i,imin,imax}] 求积Product[f,{i,imin,imax,di}] 求数列之积，引数i以di递增Product[f,{i,imin,imax},{j,jmin,jmax}]函数之泰勒展开式Series[expr,{x,x0,n}] 对expr于x0点作泰勒级数展开至(x-x0)n项Series[expr,{x,x0,m},{y,y0,n}] 对x0和y0展开关系运算子a==b 等于a>b 大于a>=b 大于等于a<b 小于a<=b 小于等于a!=b 不等于逻辑运算子!p notp||q||… orp&&q&&… andXor[p,q,…] exclusive orLogicalExpand[expr] 将逻辑表示式展开基本二维绘图指令Plot[f,{x,xmin,xmax}]画出f在xmin到xmax之间的图形Plot[{f1,f2,…},{x,xmin,xmax}]同时画出数个函数图形Plot[f,{x,xmin,xmax},option->value]指定特殊的绘图选项，画出函数f的图形Plot几种指令选项预设值说明AspectRatio 1/GoldenRatio 图形高和宽之比例，高/宽Axes True 是否把坐标轴画出AxesLabel Automatic 为坐标轴贴上标记，若设定为AxesLabel->{?ylabel?}，则为y轴之标记。
实验二怎样计算Pi
数学实验实验报告学院：数学与统计学院班级：数学与应用数学3班学号：201370010314姓名：康萍时间：2016.04.05实验二 怎样计算π一、实验目的分别用下列三种方法计算π的近似值，并比较三种方法的精确度： 数值积分法：通过使用Mathematica7.0编写梯形公式和辛普森公式的程序语言计算π。
泰勒级数法：利用反正切函数泰勒级数计算π。
蒙特卡罗（Monte Carlo ）法：通过使用Mathematica7.0编写蒙特卡罗公式的程序语言来计算π。
二、实验环境基于Windows 环境下的Mathematica7.0软件。
三、实验的基本理论和方法1、数值积分法以单位圆的圆心为原点建立直角坐标系，则单位圆在第一象限内的部分G 是一个扇形，由曲线])1,0[(12∈-=x x y 及两条坐标轴围成，它的面积4π=S 。
算出了S 的近似值，它的4倍就是π的近似值。
而扇形面积S 实际上就是定积分4112π=-⎰dx x 。
与π有关的定积分有很多，比如211x +的定积分411102π=+⎰dx x 就比21x -的定积分更容易计算，更适合于用来计算π。
一般地，要计算定积分()dx x f b a⎰，也就是计算曲线()x f y =与直线b x a x y ===,,0所围成的曲边梯形G 的面积S 。
为此，用一组平行于y 轴的直线()b x x x x x a n i x x n n i =<<<<<=-≤≤=-1210,11Λ将曲边梯形T 分成n 个小曲边梯形，总面积S 分成这些小曲边梯形的面积之和。
如果取n 很大，使每个小曲边梯形的宽度都很小，可以将它上方的边界()()i i x x x x f ≤≤-1近似的看作直线段，将每个小曲边梯形近似的看作梯形来求面积，就得到梯形公式。
如果更准确些，将每个小曲边梯形的上边界近似的看作抛物线段，就得到辛普森公式。
具体公式如下：梯形公式 设分点11,,-n x x Λ将积分区间],[b a 分成n 等份，即()n i n a b i a x i ≤≤-+=0,/。
数学实验mathematica圆周率Pi的计算方法
1.祖冲之的圆周率N 22 7,10N 355 113,103.1428571433.141592920无理数的最佳分数逼近Α Pi;Ε0 1;list ;Dop Floor q Α 0.01 ;Ε Abs p q Α q;If Ε Ε0,AppendTo list,p q ;Ε0 Ε , q,1,100000 ;list3,227,333106,355113,10399333102,10434833215,20834166317,31268999532Α Pi;Ε0 1;list ;Dop Floor q Α 0.01 ;Ε Abs p q Α ;Print "p ",p,";q ",q,";Ε ",Ε N ;If Ε Ε0,AppendTo list,p q ;Ε0 Ε , q,1,10 ;listp 3;q 1;Ε 0.141593p 6;q 2;Ε 0.283185p 9;q 3;Ε 0.424778p 12;q 4;Ε 0.566371p 15;q 5;Ε 0.707963p 18;q 6;Ε 0.849556p 22;q 7;Ε 0.00885142p 25;q 8;Ε 0.132741p 28;q 9;Ε 0.274334p 31;q 10;Ε 0.4159273,227 乐音的频率比k 2.0^ 1 12 ;music 1.0,k^2,k^4,k^5,k^7,k^9,k^11,k^121.,1.12246,1.25992,1.33484,1.49831,1.68179,1.88775,2.2数学实验圆周率Pi的计算方法.nbfreq1 1,2^ 2 12 ,2^ 4 12 ,2^ 5 12 ,2^ 7 12 ,2^ 9 12 ,2^ 11 12 ,2freq2 1,9 8,5 4,4 3,3 2,5 3,17 9,21,21 6,21 3,25 12,27 12,23 4,211 12,21,98,54,43,32,53,179,2freq freq2;m 512;Play Sin 2Pi m t freq 2 , t,0,0.8 ,PlayRange 0,1Sound SampledSoundFunction Function Play`Time3 ,Block t 0. 0.000125Play`Time3 , Sin 2Πm t freq 2 0.5 2. ,6400,8000 m 512;freq1 1,2^ 2 12 ,2^ 4 12 ,2^ 5 12 ,2^ 7 12 ,2^ 9 12 ,2^ 11 12 ,2 ;freq2 1,9 8,5 4,4 3,3 2,5 3,17 9,2 ;freq freq2;Playmusic song_ : Do Play Sin 2Pi m t freq song i ,t,0,0.8 ,PlayRange 0,1 , i,1,Length song ;music 1,2,3,4,5,6,7,8,7,6,5,4,3,2,1 ;Playmusic musicPlaysong song_ :Do x song i,1 ;w Which x 0,freq x 2,x 10,freq x 10 2,True,freq x ;y song i,2 ;Play Sin 2Pi m t w , t,0,0.4y ,PlayRange 0,1 , i,1,Length songsong2 3,4 , 5,3 , 6,1 , 1,3 , 2,1 , 6,1 ,1,1 , 5,2 , 5,3 , 11,1 , 6,1 , 5,1 , 3,1 , 5,1 , 2,8 ,2,3 , 3,1 , 7,2 , 6,2 , 5,3 , 6,1 , 1,2 , 2,2 ,3,2 , 1,2 , 6,1 , 5,1 , 6,1 , 1,1 , 5,8 ,3,3 , 5,1 , 7,2 , 2,2 , 6,1 , 1,1 , 5,3 , 3,1 ,3,1 , 5,1 , 3,2 , 5,1 , 6,1 , 7,1 , 2,1 , 6,4 ,1,2 , 1,1 , 2,1 , 5,2 , 5,1 , 3,1 , 2,2 , 3,1 , 2,1 , 1,2 ,6,1 , 5,1 , 3,4 , 1,4 , 6,1 , 1,1 , 6,1 , 5,1 ,3,1 , 5,1 , 6,1 , 1,1 , 5,4 ;Playsong song2song1 3,1 , 5,1 , 6,1 , 6,0.5 , 5,0.5 , 6,1 , 3,1 ,2,2 , 3,1 , 5,1 , 6,1 , 6,0.5 , 5,0.5 , 6,1 , 3,3 ,3,1 , 5,1 , 6,1 , 6,0.5 , 5,0.5 , 6,1 , 3,1 , 2,2 ,5,1 , 3,1 , 2,0.5 , 3,0.5 , 2,0.5 , 1,0.5 , 2,1 , 6,3 ,6,1 , 2,2 , 5,1 , 3,3 , 2,0.5 , 1,0.5 , 6,4 ,5,1 , 3,1 , 2,0.5 , 3,0.5 , 2,0.5 , 1,0.5 , 2,1 , 6,3 ,5,1 , 3,1 , 2,0.5 , 3,0.5 , 2,0.5 , 1,0.5 , 2,1 , 6,3 ;Playsong song14.单位圆的面积等于Π数学实验圆周率Pi的计算方法.nb3f x_ : 1 x2;fig Plot f x ,0 , x,0,1 ,AspectRatio 11.00.80.60.40.20.20.40.60.8 1.04数学实验圆周率Pi的计算方法.nbn 10;fig1 ;fig2 ;Do AppendTo fig1,Line 1 n i,0 , 1 n i,f 1 n i ,AppendTo fig2,Line 1 n i 1 ,f 1 n i , 1 n i,f 1 n i , Show fig,Graphics fig1 ,Graphics fig2fig3 ;fig4 ;Do AppendTo fig3,Line 1 n i,0 , 1 n i,f 1 n i 1 ,AppendTo fig4,Line 1 n i 1 ,f 1 n i 1 , 1 n i,f 1 n i 1 , i,1,n ;Show fig,Graphics fig3 ,Graphics fig4Do s1 N 4 Sum f k 10^m 10^m, k,1,10^m ;s2 N 4 Sum f k 1 10^m 10^m, k,1,10^m ;Print s1,s2, s1 s2 2 , m,1,4fig5 ;Do AppendTo fig5,Line 1 n i 1 ,f 1 n i 1 , 1 n i,f 1 n i ,i,1,n ;Show fig,Graphics fig1 ,Graphics fig5Do m 10^t;s3 N 4 f 0 f 1 2 m Sum f k m m, k,1,m 1 ,20 ;s4 N 4 f 0 f 1 m 2 Sum f k m m, k,1,m 14 Sum f k 1 2 m m, k,0,m 1 6,20 ;Print s3,s4 , t,3,4 3.1415554669110276837,3.14158751891227769063.1415914776113222011,3.14159249122019981625.数值积分法数学实验圆周率Pi的计算方法.nb5 a b f x x b a n f a f b 2 i 1n 1f x i ,x i a b a n iClear a,b,x,n ;a 0;b 1;n 100;y x_ : 4 1 x^2 ;p1 N b a n Sum y a i b a n , i,1,n 1 y a y b 2 ,50p2 N b a 6 ny a y b 2 Sum y a i b a n , i,1,n 1 4 Sum y a i 1 2 b a n , i,1,n ,50 N Pi,503.1415926535897932384626433832795028841971693993751a b f x x b a6n f a f b 4 i 1n f x i 0.5 2 i 1n 1f x i ,x i 0.5 a b a n（i 0.5 ，为 x i 1,x i 的中点Clear a,b,x,n ;a 0;b 1;n 100;y x_ : 4 1 x^2 ;p1 N b a n Sum y a i b a n , i,1,n 1 y a y b 2 ,50p2 N b a 6 ny a y b 2 Sum y a i b a n , i,1,n 1 4 Sum y a i 1 2 b a n , i,1,n ,50 N Pi,503.14159265358979323846264338327950288419716939937516.级数展开法T x_,n_ : Sum 1 ^ k 1 x^ 2k 1 2k 1 , k,1,n ;N 4 T 1,10000 ,20N Pi,203.14149265359004323853.1415926535897932384626433832795028842`20.N 4 4T 1 5,100 T 1 239,40 ,150N Pi,1503.14159265358979323853.141592653589793238462643383279502884197169399375105820974944592307816406286208998 628034825342117067982148086513282306647093844609550582231725113323563.141592653589793238462643383279502884197169399375105820974944592307816406286208998 628034825342117067982148086513282306647093844609550582231725359408137.蒙特卡罗法（随机模拟法）fig Plot 1 x2, x,0,1 ,AspectRatio 1,PlotStyle RGBColor 1,0,0 ;fig0 Graphics Line 0,1 , 1,1 , 1,0 ;Show fig,fig06数学实验圆周率Pi的计算方法.nbn 1000;fig1 ;temp 0;Do x Random ;y Random ;AppendTo fig1,Point x,y ;If x2 y2 1,temp , , i,1,n ;Show fig,fig0,Graphics fig1N temp 4 nn 10000;p ;Do m 0;Do x Random ;y Random ;If x^2 y^2 1,m m 1 , k,1,n ;AppendTo p,N 4m n , t,1,5 ;Print p ;Sum p t , t,1,5 53.1312,3.142,3.1276,3.1124,3.14483.1316。
pi的计算 实验报告
Ramanujan公式1914年，印度数学家Srinivasa Ramanujan在他的论文里发表了一系列共14条圆周率的计算公式，这是其中之一。
这个公式每计算一项可以得到8位的十进制精度。
1985年Gosper用这个公式计算到了圆周率的17,500,000位。
3、AGM(Arithmetic-Geometric Mean)算法Gauss-Legendre公式：初值：重复计算：最后计算：这个公式每迭代一次将得到双倍的十进制精度，比如要计算100万位，迭代20次就够了。
1999年9月Takahashi和Kanada用这个算法计算到了圆周率的206,158,430,000位，创出新的世界纪录。
4、Borwein四次迭代式：初值：重复计算：最后计算：这个公式由Jonathan Borwein 和Peter Borwein 于1985年发表，它四次收敛于圆周率。
5、Bailey-Borwein-Plouffe 算法014211()1681848586n n n n n n π∞==---++++∑这个公式简称BBP 公式，由David Bailey, Peter Borwein 和Simon Plouffe 于1995年共同发表。
它打破了传统的圆周率的算法，可以计算圆周率的任意第n 位，而不用计算前面的n-1位。
这为圆周率的分布式计算提供了可行性。
1997年，Fabrice Bellard 找到了一个比BBP 快40％的公式：第三部分：对于π的几种计算的研究和讨论： 1、数值积分法（I ）利用积分公式⎰-=10214dx x π计算πn=10 ans =; n=20 ans =; n=50 ans =; n=100 ans =; n=200 ans =; n=500 ans =; n=1000 ans =; n=2000 ans =;半径为1的圆称为单位圆，它的面积等于π。
只要计算出单位圆的面积，就算出了π。
mathematica 实验报告
mathematica 实验报告Mathematica 实验报告引言：Mathematica 是一款强大的数学软件，它能够帮助用户进行各种数学计算、数据分析和可视化等工作。
本实验报告将介绍我在使用 Mathematica 进行实验时的一些经验和心得。
一、实验目的本次实验的目的是通过使用 Mathematica，掌握其基本操作和功能，了解其在数学计算和数据处理方面的应用。
二、实验步骤1. 安装和启动 Mathematica首先，我在官方网站下载了 Mathematica 的安装包，并按照提示完成了安装。
然后，我启动了 Mathematica 软件，进入了主界面。
2. 基本操作在主界面中，我发现 Mathematica 提供了一个强大的交互式界面，用户可以通过键入命令和运行代码来实现各种功能。
我尝试了一些基本操作，比如进行简单的数学计算、定义变量和函数等。
3. 数据处理和分析Mathematica 提供了丰富的数据处理和分析功能，使得用户可以轻松处理和分析各种数据。
我使用了一些内置的函数和工具，对一些实验数据进行了处理和分析。
例如，我使用了 ListPlot 函数绘制了一些实验数据的散点图，并使用了Fit 函数进行了数据拟合。
4. 可视化Mathematica 还提供了强大的可视化功能，用户可以通过绘制图表和图形来展示数据和结果。
我使用了 Plot 函数绘制了一些函数的图像，并使用了 Graphics 函数绘制了一些几何图形。
5. 编程和自动化Mathematica 具有强大的编程功能，用户可以编写自己的函数和程序来实现复杂的计算和操作。
我尝试了一些简单的编程，比如编写了一个计算斐波那契数列的函数。
此外，我还了解到 Mathematica 支持自动化操作，可以通过编写脚本和批处理文件来实现自动化的计算和分析。
三、实验结果与分析通过使用 Mathematica，我成功完成了实验的各项任务，并取得了一些令人满意的结果。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
是远比 1 小。例如，因为 arctan1 arctan 1 arctan 1 ，所以我们可以计算出 arctan 1 , arctan 1
2
3
2
3
的值，从而得到 arctan1的值。这样，就使得收敛速度加快。改进后可以看出，泰勒级数
法得到的结果比数值分析法精确到小数点后更多位。
（3）蒙特卡罗法计算
0
4
利用 Mathematics 编程计算上式，过程如下：
从而得到 的近似值为 3.14159265358979323846264338328，可以看出，用这种方法 计算所得到的 值是相当精确的。n 越大，计算出来的扇形面积的近似值就越接近 的 准确值。
（2）泰勒级数法计算 利用反正切函数的泰勒级数
实验结果和结果分析：
虽然在实验过程中存在语句错写问题，但经过分析、改正均达到实验预期结果；并
在最后给出了用公式法计算π 的值，但有待于实验验证。
实验效果良好。（具体分析见实验内容与步骤）
附录：
在计算扇形面积时，很容易想到使用数学分析中积分的方法，第一象限中的扇形由
曲线 y 1 x2 (x [0,1])及两条坐标轴围成，实际操作中，我们不能准确地计算它的面
积，于是就通过分割的方法，将其划分为许多小的梯形，通过利用梯形的面积近似于扇
形面积来计算 S 1Mathematica中常用绘图函数Plot在绘制高次函数时的方法；
2、计算圆周率 的数值积分法、泰勒级数法、蒙特卡罗法，并且利用特定的公式来
计算圆周率 。
实验内容和步骤：
（1）数值积分法计算
半径为 1 的圆称为单位圆，它的面积等于 。只要计算出单位圆的面积，就算出了 。 在坐标轴上画出以圆点为圆心，以 1 为半径的单位圆（如下图），则这个单位圆在第一 象限的部分是一个扇形，而且面积是单位圆的 1/4，于是，我们只要算出此扇形的面积， 便可以计算出 。
姓名 @@@@@@@ 学院 @@@@@@@@ 班级
@@@@@@@@@ 学号 @@@@@@@@@@
实验题目
值的计算
评分
实验目的：
1、用多种方法计算圆周率 的值；
2、通过实验来体会各种方法的区别，比较各种方法的优劣；
3、尝试自己提出新的方法来计算圆周率 的值。
实验环境：
学校机房，Mathematica4.0 软件
arctan x x x3 x5 (1)k1 x 2k1
35
2k 1
来计算 。 从反正切函数的泰勒级数，进行如下编程来计算 ，实验运行如下：
从实验过程可以看出，这种方法花费的时间很长。原因是当 x=1 时得到的 arctan1的
展开式收敛太慢。要使泰勒级数收敛得快，容易想到，应当使 x 的绝对值小于 1，最好
在数值分析法中，我们利用求单位圆的 1/4 面积来得到 / 4 ，从而得到 。单位圆
的 1/4 是一个扇形，它是边长为 1 的单位正方形的一部分，单位正方形的面积 S1 1 。只
要能够求出扇形的面积 S 在正方形的面积中所占的比例 k S / S1 ，就能立即得到 S ，从 而得到 的值。下面的问题归结为如何求 k 的值，这就用到了一种利用随机数来解决此 种问题的蒙特卡罗法，其原理就是在正方形中随机的投入很多点，是所投的每个点落在 正方形中每一个位置的机会均等，看其中有多少个点落在扇形内。降落在扇形内的点的 个数 m 与所投店的总数 n 的比可以近似的作为 k 的近似值。利用 Mathetics 编程如下：
从运行结果来看，蒙特卡罗法的计算结果为 3.136，虽然精确度不太高，但运行时
间短，在很多场合下，特别是在对精确度要求不高的情况下很有用的。
（4）利用麦琴给出 4arctan1 arctan 1 ，推出π =4( 4arctan1 arctan 1 )。对比
4
5
239
5
239
以上方法，这种简单的直接用公式求的π 的方法要简单得多，所以用处更广。