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
()1681
848586n n n n n n π∞
==---++++∑
这个公式简称BBP 公式，由David Bailey, Peter Borwein 和Simon Plouffe 于1995年共同发
表。它打破了传统的圆周率的算法，可以计算圆周率的任意第n 位，而不用计算前面的n-1位。这为圆周率的分布式计算提供了可行性。1997年，Fabrice Bellard 找到了一个比BBP 快40％的公式：
第三部分：
对于π的几种计算的研究和讨论： 1、
数值积分法（I ）利用积分公式?-=1
0214dx 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的圆称为单位圆，它的面积等于π。只要计算出单位圆的面积，就算出了π。在坐标轴上画出以圆点为圆心，以1为半径的单位圆（如下图），则这个单位圆在第一象限的部分是一个扇形，而且面积是单位圆的1/4，于是，我们只要算出此扇形的面积，便可以计算出π。但计算量较大。 数值积分法（II ）利用公式?+=102
11
4dx 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 =; 设分点x 1,x 2,…x n-1将积分区间[0,1]分成n 等分。
所有的曲边梯形的宽度都是h=1/n 。记yi=f(xi).则第i 个曲边梯形的面积A 近似地等于梯形面积，即：
A=(y(i-1)+yi)h/2。
将所有这些梯形面积加起来就得到：
A ≈2/n[2(y 1+y 2+…y n-1)+y 0+y n ]
3、
利用复化梯形算法求Pi 的近似值.
1
41
x 2
x
=1/2 (T n +h
i 1
n
f a
i h
h 2
)
4、
泰勒级数法计算π：
利用反正切函数的泰勒级数ΛΛ+--+-+-=--1
2)1(53arctan 1
2153k x x x x x k k 来计算π (I)x=1时
Λ+-+-
≈7
1
513114
π
π=4* k 1
n
1k 1
2k
1
n = 10 ans= ;
n = 20 ans = ; n = 50 ans = ; n = 100 ans = ; n = 200 ans = ; n = 500 ans = ; n = 1000 ans = ; n = 2000 ans = ;
x=1时得到的arctan1的展开式收敛太慢,逼近速度太慢，运算庞大，对速度造成了很大影响。 5、
泰勒级数法计算π： 利用反正切函数的泰勒级数
ΛΛ+--+-+-=--1
2)1(53arctan 1
2153k x x x x x k k 来计算π (II)进一步精细化 11
arctan1arctan arctan 23
=+
n = 10 ans= ;
n = 20 ans = ; n = 50 ans =; n = 100 ans = ; n = 200 ans = ; n = 500 ans = ; n = 1000 ans = ; n = 2000 ans =;
当x=1时得到的arctan1的展开式收敛太慢。要使泰勒级数收敛得快，容易想到，
应当使x 的绝对值小于1，最好是远比1小。例如，因为11
arctan1arctan arctan 23
=+，
所以我们可以计算出11
arctan ,arctan 23
的值，从而得到arctan1的值。这样，就使得收敛
速度加快，逼近的速度大大增加.
利用麦琴给出
239
1arctan
51arctan 44-=π
，推出π=4(2391
arctan 51arctan 4-)
n = 10 ans= ;
n = 20 ans = ; n = 50 ans =; n = 100 ans = ; n = 200 ans = ; n = 500 ans = ; n = 1000 ans = ; n = 2000 ans =;
对泰勒级数，随着｜x ｜的减小，级数的收敛速度明显加快，这启示我们另外构造相关级数来逼近π。
7、 蒙特卡罗法计算π
单位圆的1/4是一个扇形，它是边长为1的单位正方形的一部分，单位正方形的面积
11S =。只要能够求出扇形的面积S 在正方形的面积中所占的比例1/k S S =，就能立即
得到S ，从而得到π的值。下面的问题归结为如何求k 的值，这就用到了一种利用随机数来解决此种问题的蒙特卡罗法，其原理就是在正方形中随机的投入很多点，是所投的每个点落在正方形中每一个位置的机会均等，看其中有多少个点落在扇形内。降落在扇形内的点的个数m 与所投店的总数n 的比可以近似的作为k 的近似值。
n = 100 ans= ;
n = 200 ans = ;
n = 500 ans =
n = 1000 ans = ;
n = 2000 ans = ;
n = 5000 ans =;
n = 10000 ans = ;
n = 20000 ans =;
这种数据模拟算法收敛的速度很慢，从运行结果来看，蒙特卡罗法的计算结果为，虽然精确度不太高，但运行时间短，在很多场合下，特别是在对精确度要求不高的情况下很有用的。
除以上几种方法之，还有下列一些的求其他π的方法：
*利用高斯公式
π=481
arctan
18+321
arctan
57
—201
arctan
239
*、π=2+1/3*(2+2/5*(2+3/7*(2+ ……（2+k/(2k+1)*(2+….)）)……..)))….(当k=2799时可精确到800位)
*、π/6=1/2+1/2*1/(3*2^3)+((1*3)/(2*4))*(1/(5*2^5))+……
*、e^(π*i)+1=0 (欧拉公式，也称世界上最杰出的公式) *、1+(1/2)^2+(1/3)^2+(1/4)^2+….(1/n)^2=π^2/6
*、1+(1/2)^4+(1/3)^4+(1/4)^4+….(1/n)^4=π^4/90
*、1+(1/2)^6+(1/3)^6+(1/4)^6+….(1/n)^6=π^6/945
*、1+(1/2)^8+(1/3)^8+(1/4)^8+….(1/n)^8=π^8/9450
*、1+(1/2)^10+(1/3)^10+(1/4)^10+….(1/n)^10=π^10/93555
第四部分：
求 的小数点后第位数字的几种方法：1、
反正切函数的泰勒级数
Λ
Λ+
-
-
+
-
+
-
=
-
-
1
2
)1
(
5
3
arctan
1
2
1
5
3
k
x
x
x
x
x
k
k
得到
1
2
1
)1
(
5
1
3
1
1
40+
-
=
-
+
-
=∑∞
=
k
k
k
Λ
π
；
推出：=4 *1
2
1
)1
(
+
-
∑∞
=
k
k
k
Mathematica程序
小数点后10000位数字为8
2、
沃里斯(Wallis)方法
∏∞
=
?
?
?
?
?
+
?
-
?
=
?
?
?
?
?
?
??
?
?
?
?
?
??
?
?
?
?
?
?
=
1
1
2
2
1
2
2
2
7
6
5
6
5
4
3
4
3
2
1
2
2
k
k
k
k
k
Λ
π
Mathematica 程序
小数点后10000位数字为8
3、
基于arctan x 的级数
麦琴给出：239
1arctan 51arctan 44-=π；
推出π=4(239
1
arctan 51arctan 4-)
MATLAB 程序
小数点后10000位数字为8
求取的方法很多，过程类似，不再累述。
验证
·
·
·
结论：π的小数点后10000位数字为8
第五部分：
心得体会
对于π的值我们早已知晓，但是对于这个数值如何得到却不清楚。通过这次学习，我们知道原来计算π的值有如此多的方法，并且知道了级数在一些计算中的作用。通过本次实验，我们知道了运用各种方法计算π的近似值，知
道了matlab和mathematica中关于级数的运算，学习是一个枯燥的过程，有些东西是必须懂的。也就是一些基础知识吧，包括基础的操作和相关的理论知识。因此，有一本教材看一次，必定会有所收获。看了，练习，思考。MATLAB是一个实用性很强，操作相对容易，比较完善的工具软件。使用起来比较方便，通过操作可以很快看到结果，能够清晰地感觉到成功与失败，编程问题最头疼的不是编程序，而是调程序，所以在你的程序编完之后，一定要进行验证其正确性，你要尽量多的设想你的问题的复杂性，当然，要一步一步复杂，这样才能保证你的程序的适用性很强。以后虽然没有数学应用软件设计的课程了，但是对于一个学习数学的学生来说，以后用到MATLAB 和mathematica的地方还有很多，我不会因为不用学而不去学，希望能在数学建模之前，把两个软件运用灵便，为以后更强大的数学世界打下夯实的基础。
《计算方法》课内实验报告
《计算方法》实验报告 姓名： 班级： 学号： 实验日期: 2011年10月26日
一、实验题目： 数值积分 二、实验目的： 1．熟悉matlab 编写及运行数值计算程序的方法。 2．进一步理解数值积分的基础理论。 3．进一步掌握应用不同的数值积分方法求解给定的积分并给出数据结果及误差分析。 三、实验内容： 1．分别用复合梯形求积公式及复合辛普森求积公式计算积分xdx x ln 10 ? , 要求计算精度达到410-,给出计算结果并比较两种方法的计算节点数. 2．用龙贝格求积方法计算积分dx x x ?+3 021,使误差不超过510-. 3．用3=n 的高斯-勒让德公式计算积分?3 1 sin x e x ,给出计算结果. 4．用辛普森公式(取2==M N ) 计算二重积分.5 .00 5 .00 dydx e x y ? ? - 四、实验结果： 1.（1）复合梯形法: 将区间[a,b]划分为n 等份，分点n k n a b h kh a x k ,2,1,0,,=-=+=在每个区间[1,+k k x x ](k=0,1,2,···n-1)上采用梯形公式，则得 )()]()([2)()(1 11 1 f R x f x f h dx x f dx x f I n n k k k b a n k x x k k ++===∑?∑? -=+-=+ 故)]()(2)([21 1 b f x f a f h T n k k n ++=∑-=称为复合梯形公式 计算步长和划分的区间 Eps=1E-4 h1=sqrt(Eps/abs(-(1-0)/12*1/(2+1))) h1 =0.0600 N1=ceil(1/h1) N1 =17 用复合梯形需要计算17个结点。 复合梯形： function T=trap(f,a,b,n) h=(b-a)/n;
计算方法上机实验报告
. / 《计算方法》上机实验报告 班级：XXXXXX 小组成员：XXXXXXX XXXXXXX XXXXXXX XXXXXXX 任课教师：XXX 二〇一八年五月二十五日
前言 通过进行多次的上机实验，我们结合课本上的内容以及老师对我们的指导，能够较为熟练地掌握Newton 迭代法、Jacobi 迭代法、Gauss-Seidel 迭代法、Newton 插值法、Lagrange 插值法和Gauss 求积公式等六种算法的原理和使用方法，并参考课本例题进行了MATLAB 程序的编写。 以下为本次上机实验报告，按照实验内容共分为六部分。 实验一： 一、实验名称及题目： Newton 迭代法 例2.7(P38):应用Newton 迭代法求在附近的数 值解，并使其满足. 二、解题思路： 设'x 是0)(=x f 的根,选取0x 作为'x 初始近似值,过点())(,00x f x 做曲线)(x f y =的切线L ,L 的方程为))((')(000x x x f x f y -+=,求出L 与x 轴交
点的横坐标) (') (0001x f x f x x - =,称1x 为'x 的一次近似值,过点))(,(11x f x 做曲线)(x f y =的切线,求该切线与x 轴的横坐标) (') (1112x f x f x x - =称2x 为'x 的二次近似值,重复以上过程,得'x 的近似值序列{}n x ,把) (') (1n n n n x f x f x x - =+称为'x 的1+n 次近似值，这种求解方法就是牛顿迭代法。 三、Matlab 程序代码： function newton_iteration(x0,tol) syms z %定义自变量 format long %定义精度 f=z*z*z-z-1; f1=diff(f);%求导 y=subs(f,z,x0); y1=subs(f1,z,x0);%向函数中代值 x1=x0-y/y1; k=1; while abs(x1-x0)>=tol x0=x1; y=subs(f,z,x0); y1=subs(f1,z,x0); x1=x0-y/y1;k=k+1; end x=double(x1) K 四、运行结果：
太原理工大学数值计算方法实验报告
本科实验报告 课程名称：计算机数值方法 实验项目：方程求根、线性方程组的直接解 法、线性方程组的迭代解法、代数插值和最 小二乘拟合多项式 实验地点：行勉楼 专业班级： ******** 学号： ********* 学生姓名： ******** 指导教师：李誌，崔冬华 2016年 4 月 8 日
y = x*x*x + 4 * x*x - 10; return y; } float Calculate(float a,float b) { c = (a + b) / 2; n++; if (GetY(c) == 0 || ((b - a) / 2) < 0.000005) { cout << c <<"为方程的解"<< endl; return 0; } if (GetY(a)*GetY(c) < 0) { return Calculate(a,c); } if (GetY(c)*GetY(b)< 0) { return Calculate(c,b); } } }; int main() { cout << "方程组为：f(x)=x^3+4x^2-10=0" << endl; float a, b; Text text; text.Getab(); a = text.a; b = text.b; text.Calculate(a, b); return 0; } 2.割线法： // 方程求根（割线法）.cpp : 定义控制台应用程序的入口点。// #include "stdafx.h" #include"iostream"
心得体会 使用不同的方法，可以不同程度的求得方程的解，通过二分法计算的程序实现更加了解二分法的特点，二分法过程简单，程序容易实现，但该方法收敛比较慢一般用于求根的初始近似值，不同的方法速度不同。面对一个复杂的问题，要学会简化处理步骤，分步骤一点一点的循序处理，只有这样，才能高效的解决一个复杂问题。
数值计算实验报告
(此文档为word格式，下载后您可任意编辑修改！) 2012级6班###(学号)计算机数值方法 实验报告成绩册 姓名:宋元台 学号： 成绩：
数值计算方法与算法实验报告 学期： 2014 至 2015 第 1 学期 2014年 12月1日课程名称: 数值计算方法与算法专业:信息与计算科学班级 12级5班 实验编号： 1实验项目Neton插值多项式指导教师：孙峪怀 姓名：宋元台学号：实验成绩： 一、实验目的及要求 实验目的： 掌握Newton插值多项式的算法，理解Newton插值多项式构造过程中基函数的继承特点，掌握差商表的计算特点。 实验要求： 1. 给出Newton插值算法 2. 用C语言实现算法 二、实验内容 三、实验步骤(该部分不够填写.请填写附页)
1.算法分析： 下面用伪码描述Newton插值多项式的算法： Step1 输入插值节点数n，插值点序列{x(i),f(i)},i=1,2,……,n,要计算的插值点x. Step2 形成差商表 for i=0 to n for j=n to i f(j)=((f(j)-f(j-1)(x(j)-x(j-1-i)); Step3 置初始值temp=1,newton=f(0) Step4 for i=1 to n temp=(x-x(i-1))*temp*由temp(k)=(x-x(k-1))*temp(k-1)形成 (x-x(0).....(x-x(i-1)* Newton=newton+temp*f(i); Step5 输出f（x）的近似数值newton(x)=newton. 2.用C语言实现算法的程序代码 #include
c 计算器实验报告
简单计算器 姓名: 周吉祥 实验目的：模仿日常生活中所用的计算器，自行设计一个简单的计算器程序，实现简单的计算功能。 实验内容： (1)体系设计： 程序是一个简单的计算器，能正确输入数据，能实现加、减、乘、除等算术运算，运算结果能正确显示，可以清楚数据等。 (2)设计思路： 1)先在Visual C++ 6.0中建立一个MFC工程文件，名为 calculator. 2)在对话框中添加适当的编辑框、按钮、静态文件、复选框和单 选框 3)设计按钮，并修改其相应的ID与Caption. 4)选择和设置各控件的单击鼠标事件。 5)为编辑框添加double类型的关联变量m_edit1. 6)在calculatorDlg.h中添加math.h头文件，然后添加public成 员。 7)打开calculatorDlg.cpp文件，在构造函数中，进行成员初始 化和完善各控件的响应函数代码。 (3)程序清单：
●添加的public成员： double tempvalue; //存储中间变量 double result; //存储显示结果的值 int sort; //判断后面是何种运算：1.加法2.减法3. 乘法 4.除法 int append; //判断后面是否添加数字 ●成员初始化： CCalculatorDlg::CCalculatorDlg(CWnd* pParent /*=NULL*/) : CDialog(CCalculatorDlg::IDD, pParent) { //{{AFX_DATA_INIT(CCalculatorDlg) m_edit1 = 0.0; //}}AFX_DATA_INIT // Note that LoadIcon does not require a subsequent DestroyIcon in Win32 m_hIcon = AfxGetApp()->LoadIcon(IDR_MAINFRAME); tempvalue=0; result=0; sort=0; append=0; }
计算方法实验报告
计算方法实验报告（四） 方程和方程组的迭代解法 一、实验问题 利用简单迭代法，两种加速技术，牛顿法，改进牛顿法，弦割法求解习题5-1，5-2，5-3中的一题，并尽可能准确。 选取5-3：求在x=1.5附近的根。 二、问题的分析（描述算法的步骤等） （1）简单迭代法算法： 给定初始近似值,求的解。 Step 1 令i=0； Step 2 令(计算)； Step 3 如果,则迭代终止，否则重复Step 2。 （2）Aitken加速法算法 Step 1 令k=0，利用简单迭代算法得到迭代序列； Step 2 令-(计算得到一个新的序列，其中k=0,1,2…)；Step 3 如果,则迭代终止，否则重复Step 2。 (3)插值加速法算法 Step 1 令k=0，利用简单迭代算法得到迭代序列； Step 2 令+(计算得到一个新的序列，其中k=1,2,3…)； Step 3 如果,则迭代终止，否则重复Step 2。 (4)牛顿法算法
Step 1给定初始近似值； Step 2令,其中k计算得到的序列； Step 3如果,则迭代终止，否则重复Step 2。 （5）改进牛顿法的算法 Step 1给定初始近似值； Step 2令,其中k迭代计算得到的序列； Step 3如果,则迭代终止，否则重复Step 2。 （6）弦割法算法（双点弦割法） Step 1给定初始近似值，； Step 2令其中k计算得到的序列； Step 3如果,则迭代终止，否则重复Step 2。 三、程序设计 (1)简单迭代法 利用迭代公式进行迭代运算。 #include
PI值计算方法1
计算圆周率的一些公式 -|waruqi 发表于 2005-12-8 9:24:00 Machin公式 这个公式由英国天文学教授John Machin于1706年发现。他利用这个公式计算到了100位的圆周率。Machin公式每计算一项可以得到1.4位的十进制精度。因为它的计算过程中被乘数和被除数都不大于长整数，所以可以很容易地在计算机上编程实现。 还有很多类似于Machin公式的反正切公式: pi/4=arctg(1/2)+arctg(1/5)+ arctg(1/8) 1844.达塞利 = arctg(1/2)+ arctg(1/3) =2 arctg(1/3)+ arctg(1/7) =12 arctg(1/18)+8 arctg(1/57)-5 arctg(1/239) 在所有这些公式中，Machin公式似乎是最快的了。虽然如此，如果要计算更多的位数，比如几千万位，Machin公式就力不从心了。下面介绍的算法，在PC机上计算大约一天时间，就可以得到圆周率的过亿位的精度。这些算法用程序实现起来比较复杂。因为计算过程中涉及两个大数的乘除运算，要用FFT(Fast Four ier Transform)算法。FFT可以将两个大数的乘除运算时间由O(n2)缩短为O(nl og(n))。 （FFT算法不在此文讲诉）
Ramanujan公式 1914年，印度数学家Srinivasa Ramanujan在他的论文里发表了一系列共1 4条圆周率的计算公式，这是其中之一。这个公式每计算一项可以得到8位的十进制精度。1985年Gosper用这个公式计算到了圆周率的17,500,000位。 1989年，David & Gregory Chudnovsky兄弟将Ramanujan公式改良成为： 这个公式被称为Chudnovsky公式，每计算一项可以得到15位的十进制精度。1994年Chudnovsky兄弟利用这个公式计算到了4,044,000,000位。Chudnovsky 公式的另一个更方便于计算机编程的形式是： AGM(Arithmetic-Geometric Mean)算法 Gauss-Legendre公式： 初值： 重复计算：
计算方法实验报告格式
计算方法实验报告格式 小组名称: 组长姓名(班号)： 小组成员姓名(班号): 按贡献排序情况: 指导教师评语: 小组所得分数: 一个完整的实验，应包括数据准备、理论基础、实验内容及方法，最终对实验结果进行分析，以达到对理论知识的感性认识，进一步加深对相关算法的理解，数值实验以实验报告形式完成，实验报告格式如下： 一、实验名称 实验者可根据报告形式需要适当写出． 二、实验目的及要求 首先要求做实验者明确，为什么要做某个实验，实验目的是什么，做完该实验应达到什么结果，在实验过程中的注意事项，实验方法对结果的影响也可以以实验目的的形式列出． 三、算法描述（实验原理与基础理论） 数值实验本身就是为了加深对基础理论及方法的理解而设置的，所以要求将实验涉及到的理论基础，算法原理详尽列出． 四、实验内容 实验内容主要包括实验的实施方案、步骤、实验数据准备、实验的算法以及可能用到的仪器设备． 五、程序流程图 画出程序实现过程的流程图，以便更好的对程序执行的过程有清楚的认识，在程序调试过程中更容易发现问题． 六、实验结果 实验结果应包括实验的原始数据、中间结果及实验的最终结果，复杂的结果可以用表格
形式列出，较为简单的结果可以与实验结果分析合并出现． 七、实验结果分析 实验结果分析包括对对算法的理解与分析、改进与建议． 数值实验报告范例 为了更好地做好数值实验并写出规范的数值实验报告，下面给出一简单范例供读者参考． 数值实验报告 小组名称: 小组成员(班号): 按贡献排序情况: 指导教师评语: 小组所得分数: 一、实验名称 误差传播与算法稳定性． 二、实验目的 1．理解数值计算稳定性的概念． 2．了解数值计算方法的必要性． 3．体会数值计算的收敛性与收敛速度． 三、实验内容 计算dx x x I n n ? += 1 10 ，1,2,,10n = ． 四、算法描述 由 dx x x I n n ? += 1 10 ，知 dx x x I n n ?+=--101110，则
数学计算方法实验报告
数学计算方法实验报告 习题二 2.估计用二分法求方程f(x)=x3+4x2-10=0在区间[1,2]内根的近似值，为使方程不超过10时所需的二分次数。f(x k) 程序过程： function two (tolerance) a=1;b=2;counter=0; while (abs(b-a)>tolerance) c=(a+b)/2; fa=a^3+4*a^2-10;
fb=b^3+4*b^2-10; fc=c^3+4*c^2-10; if ((fa==0|fb==0)) disp(counter); elseif (fa*fc<0) b=c;counter=counter+1; elseif (fb*fc<0) a=c;counter=counter+1; elseif (fb==0) disp(counter); end end solution=(a+b)/2; disp(solution); disp(counter); 实验结果： 6.取x0=1.5,用牛顿迭代法求第三中的方程根.f(x)=x3+4x2-10=0的近似值（精确到||x k+1-x k|≦10-5,并将迭代次数与3题比较。 程序过程： function six (g) a=1.5; fa=a^3+4*a^2-10;
ga=3*a^2+8*a; b=a-fa/ga; k=1; while(abs(b-a)>g) a=b; fa=a^3+4*a^2-10; ga=3*a^2+8*a; b=a-fa/ga; k=k+1; end format long; disp(a); disp(k); 实验结果：程序结果计算结果 8.用弦割法求方程f(x)=x3-3x2-x+9=0在区间[-2,-1]内的一个实根近似值x k,|f(x k)|≦10-5. 程序过程： function eight (t) a=-2; b=-1; fa=a^3-3*a^2-a+9; fb=b^3-3*b^2-b+9; c=b-fb*(b-a)/(fb-fa); k=1; while(abs(c-b)>t) a=b; b=c; fa=a^3-3*a^2-a+9; fb=b^3-3*b^2-b+9; c=b-fb*(b-a)/(fb-fa); k=k+1; end
计算pi
一、实验目的 探索精确计算π值的方法，并且比较不同方法之间的不同之处和优缺点。掌握数值积分的辛普森公式。 二、问题描述 1. 任务1 1) 用反正切函数的幂级数展开式结合有关公式求π，若要精确到40位、50位数 字，试比较简单公式和Machin 公式所用的项数。 2) 验证公式 111=arctan arctan arctan 4 258π ++ 试试此公式右端做幂级数展开完成任务1所需要的步数。 2. 任务2 用数值积分计算π，分别用梯形法和Simpson 法精确到10位数字，用Simpson 法精确到15位数字。 3. 任务3 用Monte Carlo 法计算π，除了加大随机数，在随机数一定时可重复算若干次后求平均值，看能否求得5位精确数字？ 设计方案用计算机模拟Buffon 实验 4. 任务4 利用积分 2 0(1)!!sin !!2 n n xdx n π π-=? ，n 为奇数 推导公式 224422213352121 n n n n π=-+ ……… 用此公式计算π的近似值，效果如何？ 5. 任务5 利用学过的知识（或查阅资料），提出其他计算π的方法（先用你学过的知识证明），然后实践这种方法。 对你在实验中应用的计算π的方法进行比较分析。 6. 任务6 e 是一个重要的超越数 1e lim 1)n n n →∞=+（ 1111...2!!(1)! e e n n θ =++++++ 试用上述公式或其他方法近似计算e 。
三、问题解法 1. 任务1 1) 根据幂级数展开的相关知识，易知： 24122211(1)1n n x x x x --=-+-+-++……… 因为2 1(arctan )'1x x =+，故可以求得arctan x 的幂级数展开式为： 35 211arctan (1)3521 n n x x x x x n --=-+-+-+-……… 当x=1时， -11111--(-1)4352-1 n n π=+??++? 当叠加了十万次以后得到结果π=3.141582654…只有五位有效数字，可见其精度与效率极低。如果想要精确计算π的数值的话，非常有必要寻找改进以后的方法，这就引出了两个能够提高计算效率的公式—— 简单公式： 11=arctan arctan 4 23 π + Machin 公式： 11=4arctan arctan 45239 π- 对以上两式进行arctan 的幂级数展开可以非常快速的求得π比较精确的数值。下面比较π精确到40位和50位数字时两个公式各需要计算多少项。 用简单公式得到40位有效数字需要叠加62项： 3.141592653589793238462643383279502884197 用简单公式得到50位有效数字需要叠加79项： 3.1415926535897932384626433832795028841971693993751 用Machin 公式得到40位有效数字需要叠加27项： 3.141592653589793238462643383279502884197 用Machin 公式得到50位有效数字需要叠加35项： 3.1415926535897932384626433832795028841971693993751 从上面简单的对比可以看出Machin 公式要优于简单公式，简单公式要优于不用公式的arctan 幂级数展开。在得到相同精度的条件下，Machin 公式所需要的叠加步数要显著少于简单公式，并且在计算精度越高的情况下，优势越明显。道理很简单，因为 Machin 公式计算的收敛速度要显著快于普通公式。简单公式决定收敛速度的是12n ?? ??? ，而Machin 公式决定收敛速度的是15n ?? ???，因为15n ?? ???的收敛速度快于12n ?? ??? ，故Machin 公式计算pi 的时候收敛速度要快于普通公式。所以Machin 公式比普通公式更加精确，并且在计算高精度的时候有更大优势。 2) 根据三角函数公式有：
计算方法实验报告 拟合
南京信息工程大学实验（实习）报告 一、实验目的: 用最小二乘法将给定的十个点拟合成三次多项式。 二、实验步骤： 用matlab编制以函数为基的多项式最小二乘拟合程序，并用于对下列数据作三次多项式最小二乘拟合(取权函数wi=1) x -2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 y -2.30 -1 -0.14 -0.25 0.61 1.03 1.75 2.75 4.42 6.94 给定直线方程为：y=1/4*x3+1/2*x2+x+1 三、实验结论： 最小二乘法：通过最小化误差的平方和寻找数据的最佳函数匹配。利用最小二乘法可以简便地求得未知的数据，并使得这些求得的数据与实际数据之间误差的平方和为最小。最小二乘法还可用于曲线拟合。 一般地。当测量数据的散布图无明显的规律时，习惯上取n次代数多项式。 程序运行结果为： a = 0.9731 1.1023 0.4862 0.2238 即拟合的三次方程为：y=0.9731+1.1023x+0.4862*x2+0.2238*x3
-2.5 -2-1.5-1-0.5 00.51 1.52 2.5 -4-20246 81012 x 轴 y 轴 拟合图 离散点 y=a(1)+a(2)*x+a(3)*x.2+a(4)*x.3 结论： 一般情况下，拟合函数使得所有的残差为零是不可能的。由图形可以看出最小二乘解决了残差的正负相互抵消的问题，使得拟合函数更加密合实验数据。 优点：曲线拟合是使拟合函数和一系列的离散点与观测值的偏差平方和达到最小。 缺点：由于计算方法简单，若要保证数据的精确度，需要大量的数据代入计算。
计算方法实验报告
实验报告 一、求方程f(x)=x^3-sinx-12x+1的全部根, ε=1e -6 1、 用一般迭代法; 2、 用牛顿迭代法; 并比较两种迭代的收敛速度。 一、首先，由题可求得：12cos 3)(2 ' --=x x x f . 其次，分析得到其根所在的区间。 ① 令()0=x f ，可得到x x x sin 1123 =+-. ② 用一阶导数分析得到1123 +-x x 和x sin 两个函数的增减区间；再用二阶导数分析得到 两个函数的拐点以及凹凸区间. ③ 在直角坐标轴上描摹出01123 =+-x x 和0sin =x 的图，在图上可以看到他们的交点，然后估计交点所在的区间，即是所要求的根的区间。经过估计，得到根所在的区间为 []3,4--，[]1,0和[]4,3. 1、 一般迭代法 （1）算法步骤： 设ε为给定的允许精度，迭代法的计算步骤为： ① 选定初值0x .由()0=x f 确定函数()x g ，得等价形式()x g x =. ② 计算()0x g .由迭代公式得()01x g x =. ③ 如果ε≤-01x x ，则迭代结束，取1x 为解的近似值；否则，用1x 代替0x ，重复步骤②和步骤③. （2）程序代码： ① 在区间[]3,4--内， 代码： clc
x0=-3.5; %初值0x iter_max=100; %迭代的最大次数 ep=1e-6; %允许精度 ε k=0; while k<=iter_max %k 从0开始到iter_max 循环 x1=(sin(x0)+12*x0-1).^(1/3); %代入0x ，算出1x 的值 if abs(x1-x0) 上机7：文件操作 一、问题： 古人计算pi （3.1415926535 89793 23846……），一般是用割圆法。即用圆的内接或外切正多边形来逼近圆的周长。大约公元前1200年，中国人计算到小数点后1位；Archimedes 在公元前250年用正96边形得到pi 小数点后3位的精度；公元263年，刘徽用正3072边形计算到小数点后5位 ； 公元480年，祖冲之计算到小数点后7位；Ludolph Van Ceulen 在1609年用正262边形得到了35位 精度；1706年，John Machin 计算到小数点后100位。1874年，William Shanks 计算到小数点后707位（527位正确）。这种基于几何的算法计算量大，速度慢，吃力不讨好。 1973年Guilloud & Bouyer 用 CDC 7600计算机计算到1,001,250位；1989年Kanada & Tamurar 用 HITACHI S-820/80计算机计算到1,073,741,799位；1999年Takahashi & Kanada 用 HITACHI SR8000计算机计算到206,158,430,000位；pi 的最新计算纪录由两位日本人Daisuke Takahashi 和Yasumasa Kanada 所创造。他们在日本东京大学的IT 中心，以Gauss-Legendre 算法编写程序，利用一台每秒可执行一万亿次浮点运算的超级计算机，从日本时间1999年9月18日19:00:52起，计算了37小时21分04秒，得到了pi 的 206,158,430,208(3*236)位十进制精度，之后和他们于1999年6月27日以Borwein 四次迭代式计算 了46小时得到的结果相比，发现最后45位小数有差异，因此他们取小数点后206,158,430,000位的p 值为本次计算结果。这一结果打破了他们于1999年4月创造的68,719,470,000位的世界纪录。随着数学的发展，数学家们在进行数学研究时有意无意地发现了许多计算pi 的公式。 Machin 公式 1 15239164arg arctg tg π=- 3 5 7 211()(1)35721 n n x x x x a rctg x x n --=-+-++-- 这个公式由英国天文学教授John Machin 于1706年发现。他利用这个公式计算到了100位的pi 值。Machin 公式每计算一项可以得到1.4位的十进制精度。因为它的计算过程中被乘数和被除数都不大于长整数，所以可以很容易地在计算机上编程实现。 Bailey-Borwein-Plouffe 算法 0142 1 1 ()1681848586n n n n n n π∞ ==---++++∑ 这个公式简称BBP 公式，由David Bailey, Peter Borwein 和Simon Plouffe 于1995年共同发表。它打破了传统的p 的算法，可以计算p 的任意第n 位，而不用计算前面的n-1位。这为p 的分布式计算提供了可行性。 因此使用这个公式计算，与计算常数e （2.718282……）的公式非常类似，编程很容易实现。 二、要求： 编写程序计算常数π到小数点后任意位，并将数据写入文件中。 参考程序： #include 《数据分析方法与技术》上机实验——实验1描述性统计方法 学号： 姓名： 日期： 实验项目（一）：描述性统计方法 一、实验内容 1．实验目的 掌握常用的描述性图表展示方法的原理及操作，包括：频数分布表、分组频数表、列联表、茎叶图、箱线图、误差图、散点图等； 掌握常用的描述性统计方法的原理及操作，包括：算术平均值、中位数、众数、四分位数、极差、平均差、方差、标准差、标准分数、离散系数等。 2. 实验内容和要求 实验内容：基于标准数据集，属性描述性图表展示方法（数分布表、分组频数表、列联表、茎叶图、箱线图、误差图、散点图等），对统计指标（算术平均值、中位数、众数、极差、平均差、方差、标准差、标准分数、离散系数、偏态峰态）进行计算。 实验要求：掌握各种描述性统计指标的计算思路及其在SPSS或EXCEL环境下的操作方法，掌握输出结果的解释。 二、实验过程 1、数据集介绍 1.数据库标题:鲍鱼数据 2.该数据库共计4177行数据 3.该数据有八个属性（包含性别共有九项） 4.以下是关于属性的描述，包括属性的名称，数据类型，测量单元和一个简短的描述： Name Data TypeMeas.Description ---- --------- ----- ----------- Sex nominal M, F, and I (infant)鲍鱼宝宝 Length continuousmm Longest shell measurement最长壳 Diameter continuousmm perpendicular to length垂直长度 Height continuousmm with meat in shell有肉的壳高度 Whole weightcontinuousgramswhole abalone整个鲍鱼 Shucked weightcontinuousgramsweight of meat肉的重量 Viscera weightcontinuousgramsgut weight (after bleeding)放血后内脏重 Shell weightcontinuousgramsafter being dried弄干后重量 Rings integer +1.5 gives the age in years +1.5=年龄 5.数据的值域 怎样计算 姓名： 学号 班级：数学与应用数学4班 实验报告 实验目的：自己尝试利用Mathematica软件计算的近似值，并学会计算的近似值的方法。 实验环境：Mathematica软件 实验基本理论和方法： 方法一：数值积分法（单位圆的面积是，只要计算出单位圆的面积也就计算出了的值） 其具体内容是：以单位圆的圆心为原点建立直角坐标系，则单位圆在第一象限内的部分G是一个扇形， 由曲线（）及坐标轴围成，它的面积是，算出了S的近似值，它的4倍就是的近似值。而怎样计算扇形G的面积S的近似值呢?如图 图一 扇形G中，作平行于y轴的直线将x轴上的区间[0，1]（也就是扇形在x轴上的半径）分成n等份（n=20），相应的将扇形G分成n个同样宽度1/n的部分（）。每部分是一个曲边梯形：它的左方、右方的边界是相互平行的直线段，类似于梯形的两底；上方边界是一段曲线，因此称为曲边梯形。如果n很大，每个曲边梯形的上边界可以近似的看成直线段，从而将近似的看成一个梯形来计算它的面积；梯形的高（也就是它的宽度）h=1/n，两条底边的长分别是和，于是这个梯形面积可以作为曲边梯形面积的近似值。所有这些梯形面积的和T就可以作为扇形面积S的近似值： n越大，计算出来的梯形面积之和T就越接近扇形面积S，而4T就越接近的准确值。 方法二：泰勒级数法 其具体内容是：利用反正切函数的泰勒级数 计算。 方法三：蒙特卡罗法 其具体内容是：单位正方形的面积=1，只要能够求出扇形G 的面积S在正方形的面积中所占的比例，就能立即得到S，从而得到的值。而求扇形面积在正方形面积中所占的比例k的值，方法是在正方形中随机地投入很多点，使所投的每个点落在正方形中每一个位置的机会均等，看其中有多少个点落在扇形内。将落在扇形内的点的个数m与所投的点的总数n的比可以作为k的近似值。能够产生在区间[0，1]内均匀分布的随机数，在Mathematica中语句是 Random[ ] 产生两个这样的随机数x，y，则以（x，y）为坐标的点就是单位正方形内的一点P，它落在正方形内每一个位置的机会均等。P落在扇形内的充分必要条件是。这样利用随机数来解决数学问题的方法叫蒙特卡罗法。 实验内容、步骤及其结果分析： 问题1：在方法一中，取n=1000，通过计算图一中扇形面积计算的的近似值。 分析：图一中的扇形面积S实际上就是定积分。 与有关的定积分很多，比如的定积分 实验一——插值方法 实验学时：4 实验类型：设计 实验要求：必修 一 实验目的 通过本次上机实习，能够进一步加深对各种插值算法的理解；学会使用用三种类型的插值函数的数学模型、基本算法，结合相应软件（如VC/VB/Delphi/Matlab/JAVA/Turbo C ）编程实现数值方法的求解。并用该软件的绘图功能来显示插值函数，使其计算结果更加直观和形象化。 二 实验内容 通过程序求出插值函数的表达式是比较麻烦的，常用的方法是描出插值曲线上尽量密集的有限个采样点，并用这有限个采样点的连线，即折线，近似插值曲线。取点越密集，所得折线就越逼近理论上的插值曲线。本实验中将所取的点的横坐标存放于动态数组[]X n 中，通过插值方法计算得到的对应纵坐标存放 于动态数组[]Y n 中。 以Visual C++.Net 2005为例。 本实验将Lagrange 插值、Newton 插值和三次样条插值实现为一个C++类CInterpolation ，并在Button 单击事件中调用该类相应函数，得出插值结果并画出图像。CInterpolation 类为 class CInterpolation { public : CInterpolation();//构造函数 CInterpolation(float *x1, float *y1, int n1);//结点横坐标、纵坐标、下标上限 ~ CInterpolation();//析构函数 ………… ………… int n, N;//结点下标上限，采样点下标上限 float *x, *y, *X;//分别存放结点横坐标、结点纵坐标、采样点横坐标 float *p_H,*p_Alpha,*p_Beta,*p_a,*p_b,*p_c,*p_d,*p_m;//样条插值用到的公有指针，分别存放 i h ，i α，i β，i a ，i b ，i c ，i d 和i m }; 其中，有参数的构造函数为 CInterpolation(float *x1, float *y1, int n1) { //动态数组x1，y1中存放结点的横、纵坐标，n1是结点下标上限（即n1+1个结点） n=n1; N=x1[n]-x1[0]; X=new float [N+1]; x=new float [n+1]; y=new float [n+1]; 实验四你会用几种方法计算PI（圆周率）的值 一、问题分析 若想计算π的值，就要将跟π有关联的联系在一起，找到与π近似等价的式子，利用计算其值来得到π的值，还有对于含有π的面积、体积等关系式，可以尽量使用较规则的图形来代替进行面积、体积的求解。 二、模型建立 2.1数值积分法 找一个积分值等于π的定积分,则只要利用定积分计算出的值，就可以得到π的近似值。2.2幂级数法 利用arctanx的泰勒展开式，计算π的近似值。 当x=1时，arctan1= 2.3迭代法 1976年的迭代算法： 2.4 随机模拟法（蒙特卡罗方法） 用随机模拟求单位圆面积 向单位正方形随机投n块小石头，n很大时小石头大致均匀第分布在正方形中，如果有k块落在单位圆内，单位圆面积的近似值 三、解决问题所需的基本理论和方法 (1)对于定积分,则只要计算出的值，就可以得到π的近似值,也就是计算出与直线 y=0,x=0,x=1所围成的曲边梯形，而对于此类计算往往采用数值积分梯形公式计算。 梯形公式：将积分区间n等分将所有梯形面积加起来得到 Trapz(x)：输出数组x，输出按梯形公式x的积分(单位步长) Trapz(x,y):计算y对x的梯形积分，其中x、y定义函数关系y=f(x) (2)利用arctanx的泰勒展开式，计算π的近似值。 函数taylor用于实现Taylor级数 r=taylor(f,n,v),指定自变量v和阶数n r= taylor(f,n,v,a),指定自变量v、阶数n，计算f在a的级数 (3)级数法 由于利用arctanx的幂级数展开法的收敛较慢，可采用公式 的计算来求pi值。 (4)特殊公式（BBP） 四、设计算法、编程求解 4.1数值积分法 梯形公式Matlab代码： format long x=0:0.1:1; % x=0:0.01:1; x=0:0.02:1; x=0:0.001:1; x=0:0.0001:1; y=sqrt(1-x.^2); pi=4*trapz(x,y) MAtlab结果： 数值积分法（梯形公式）Pi值 n=10 3.104518326248318 n=100 3.140417031779045 n=500 3.141487477002141 n=1000 3.141555466911028 n=10000 3.141591477611324 4.2幂级数法Matlab代码： (1) format long syms x f=atan(x); t= taylor(f,10,x,0); % t= taylor(f,100,x,0); t= taylor(f,500,x,0); t= taylor(f,1000,x,0); t= taylor(f,10000,x,0); x=1; pi=4*eval(t) 结果： 级数法（arctanx）Pi值 n=10 3.339682539682540 n=100 3.121594652591011 n=500 3.137592669589472 n=1000 3.139592655589785 计算方法上机实验报告——拉格朗日插值问题 一、方法原理 n次拉格朗日插值多项式为：Ln(x)=y0l0(x)+y1l1(x)+y2l2(x)+…+ynln(x) n=1时，称为线性插值，L1(x)=y0(x-x1)/(x0-x1)+y1(x-x0)/(x1-x0)=y0+(y1-x0)(x-x0)/(x1-x0) n=2时，称为二次插值或抛物线插值，精度相对高些 L2(x)=y0(x-x1)(x-x2)/(x0-x1)/(x0-x2)+y1(x-x0)(x-x2)/(x1-x0)/(x1-x 2)+y2(x-x0)(x-x1)/(x2-x0)/(x2-x1) 二、主要思路 使用线性方程组求系数构造插值公式相对复杂，可改用构造方法来插值。 对节点xi(i=0,1,…,n)中任一点xk(0<=k<=n)作一n次多项式lk(xk)，使它在该点上取值为1，而在其余点xi(i=0,1,…,k-1,k+1,…,n)上为0，则插值多项式为Ln(x)=y0l0(x)+y1l1(x)+y2l2(x)+…+ynln(x) 上式表明：n个点xi(i=0,1,…,k-1,k+1,…,n)都是lk(x)的零点。可求得lk 三．计算方法及过程：1.输入节点的个数n 2.输入各个节点的横纵坐标 3.输入插值点 4.调用函数，返回z 函数语句与形参说明 程序源代码如下： 形参与函数类型 参数意义 intn 节点的个数 doublex[n]（double*x） 存放n个节点的值 doubley[n]（double*y） 存放n个节点相对应的函数值 doublep 指定插值点的值 doublefun() 函数返回一个双精度实型函数值，即插值点p处的近似函数值 #includePI的计算算法
数据分析方法与技术- 实验报告模板
数学实验：怎样计算圆周率
计算方法实验报告册
计算pi
计算方法上机实验报告——拉格朗日插值问题