Matlab计算圆周率

格式：ppt
大小：718.50 KB
文档页数：30

下载文档原格式

拉格朗日插值多项式积分求圆周率近似Matlab实现

Lagrange 插值多项式积分求圆周率近似摘要：公式1：y1=4/(1+x^2) 公式2：y2=4*sqrt(1-x^2) 分别对公式1、公式2求其拉格朗日插值多项式，再对其求0-1上的定积分来求圆周率π的近似值，并在Matlab 中通过画图来比较两个所求得的值与真实值π的偏差。

Lagrange 插值多项式：)()(l )(L 0i n i ni x f x x ∑==其中 )())(())(()())(())(()(l 11101110i n i i i i i i i n i i x x x x x x x x x x x x x x x x x x x x x -⋯⋯--⋯⋯---⋯⋯--⋯⋯--=+-+-)(i x f 为函数在i x 处的函数值，)(x L n 为Lagrange 插值多项式。

Matlab 实现：clc;clear;a=0;b=1;n=input('Enter a number n:'); %将0-1分割成n 节点，即n-1段 X=zeros(1,n); %用来放置节点x 的值P=zeros(1,n); %用来放置节点x 对应的函数值y1 Q=zeros(1,n); %用来放置节点x 对应的函数值y2 x=0;h=(b-a)/(n-1); %h 为步长for i=1:ny1=4/(1+x^2);y2=4*sqrt(1-x^2);X(i)=x;P(i)=y1;Q(i)=y2;x=x+h;endX;P;Q; %通过循环对X、P、Q进行赋值syms s;l=1;z1=0;z2=0;for j=1:1for k=2:nl=l*(s-X(k))/(X(j)-X(k));endz1=z1+l*P(j);z2=z2+l*Q(j);endfor j=2:nl=1;for k=1:j-1l=l*(s-X(k))/(X(j)-X(k));endfor k=j+1:nl=l*(s-X(k))/(X(j)-X(k));endz1=z1+l*P(j); %通过循环求的函数y1的Lagrange插值多项式z1 z2=z2+l*Q(j); %通过循环求的函数y2的Lagrange插值多项式z2 endI1=int(z1,s,0,1); % z1对s在0-1上求定积分I1=eval(I1) %用小数形式表示I1I2=int(z2,s,0,1); % z2对s在0-1上求定积分I2=eval(I2) %用小数形式表示I2x=3.10:0.0001:3.20;y0=pi;y1=I1;y2=I2;plot(x,y0,'r') %红线为圆周率π的真实值hold onplot(x,y1,'g') %绿线为公式1所求值hold onplot(x,y2,'b') %蓝线为公式2所求值运行结果：从图中可以看出，当n=6时绿线很接近红线即圆周率π的真实值，而蓝线则偏离较远，当n=11时，绿线基本与红线重叠，而蓝线相对之前来说也减小偏差。

自己动手计算圆周率

自己动手计算圆周率圆周率的计算历程圆周率是一个极其著名的数。

从有文字记载的历史开始，这个数就引进了外行人和学者们的兴趣。

作为一个非常重要的常数，圆周率最早是出于解决有关圆的计算问题。

仅凭这一点，求出它的尽量准确的近似值，就是一个极其迫切的问题了。

事实也是如此，几千年来作为数学家们的奋斗目标，古今中外一代又一代的数学家为此献出了自己的智慧和劳动。

回忆历史，人类对π 的认识过程，反映了数学和计算技术开展情形的一个侧面。

π 的研究，在一定程度上反映这个地区或时代的数学水平。

德国数学史家康托说：“历史上一个国家所算得的圆周率的准确程度，可以作为衡量这个国家当时数学开展水平的指标。

〞直到19世纪初，求圆周率的值应该说是数学中的头号难题。

为求得圆周率的值，人类走过了漫长而曲折的道路，它的历史是饶有趣味的。

我们可以将这一计算历程分为几个阶段。

〔1〕实验时期通过实验对π 值进行估算，这是计算π 的的第一阶段。

这种对π 值的估算根本上都是以观察或实验为根据，是基于对一个圆的周长和直径的实际测量而得出的。

在古代世界，实际上长期使用π ＝3这个数值。

在我国刘徽之前“圆径一而周三〞曾广泛流传。

我国第一部《周髀算经》中，就记载有圆“周三径一〞这一结论。

在我国，木工师傅有两句从古流传下来的口诀：叫做：“周三径一，方五斜七〞，意思是说，直径为1的圆，周长大约是3，边长为5的正方形，对角线之长约为7。

这正反映了早期人们对圆周率π 和√2 这两个无理数的粗略估计。

东汉时期官方还明文规定圆周率取3为计算面积的标准。

后人称之为“古率〞。

在我国东、西汉之交，新朝王莽令刘歆制造量的容器――律嘉量斛。

刘歆在制造标准容器的过程中就需要用到圆周率的值。

为此，他大约也是通过做实验，得到一些关于圆周率的并不划一的近似值。

现在根据铭文推算，其计算值分别取为3.1547，3.1992，3.1498，3.2031比径一周三的古率已有所进步。

人类的这种探索的结果，当主要估计园田面积时，对生产没有太大影响，但以此来制造器皿或其它计算就不适宜了。

(完整版)MATLAB常用函数总结,推荐文档

MATLAB 常用函数总结Matlab 的内部常数pi 圆周率exp(1)自然对数的底数ei 或j虚数单位Inf 或 inf无穷大Matlab 的常用内部数学函数指数函数exp(x)以e 为底数log(x)自然对数，即以e 为底数的对数log10(x)常用对数，即以10为底数的对数对数函数log2(x)以2为底数的x 的对数开方函数sqrt(x)表示x 的算术平方根绝对值函数abs(x)表示实数的绝对值以及复数的模sin(x)正弦函数cos(x)余弦函数tan(x)正切函数cot(x)余切函数sec(x)正割函数三角函数（自变量的单位为弧度）csc(x)余割函数反三角函数asin(x)反正弦函数acos(x)反余弦函数atan(x)反正切函数acot(x)反余切函数asec(x)反正割函数acsc(x)反余割函数sinh(x)双曲正弦函数cosh(x)双曲余弦函数tanh(x)双曲正切函数coth(x)双曲余切函数sech(x)双曲正割函数双曲函数csch(x)双曲余割函数asinh(x)反双曲正弦函数acosh(x)反双曲余弦函数atanh(x)反双曲正切函数acoth(x)反双曲余切函数asech(x)反双曲正割函数反双曲函数acsch(x)反双曲余割函数求角度函数atan2(y,x)以坐标原点为顶点，x轴正半轴为始边，从原点到点（x，y）的射线为终边的角，其单位为弧度，范围为（，]gcd(a,b)两个整数的最大公约数数论函数lcm(a,b)两个整数的最小公倍数排列组合函数factorial(n)阶乘函数，表示n的阶乘real(z)实部函数imag(z)虚部函数复数函数abs(z)求复数z的模angle(z)求复数z 的辐角，其范围是（，]conj(z)求复数z 的共轭复数ceil(x)表示大于或等于实数x 的最小整数floor(x)表示小于或等于实数x 的最大整数求整函数与截尾函数round(x)最接近x 的整数max([a ，b ，c ，．．．])求最大数最大、最小函数min([a ，b ，c ，．．])求最小数符号函数sign(x)Matlab 中的数学运算符a+b 加法 a./b 数组右除a-b 减法 a.\b 数组左除a*b 矩阵乘法a^b 矩阵乘方a.*b 数组乘法 a.^b 数组乘方a/b 矩阵右除-a负号a\b矩阵左除' 共轭转置.'一般转置Matlab 的关系运算符 ==等于<小于>大于<=小于或等于>=大于或等于~=不等于如何用matlab求阶乘factorial（n）　求n的阶乘如何用matlab进行多项式运算（1）合并同类项 syms 表达式中包含的变量 collect(表达式，指定的变量)（2）因式分解　syms 表达式中包含的变量 factor(表达式)（3）展开syms 表达式中包含的变量 expand(表达式)（4）化简syms 表达式中包含的变量simplify(表达式)　　如何用matlab进行复数运算　a+b*i 或 a +b*j表示复数a+bi 或 a+bjreal（z）求复数z的实部imag（z）求复数z的虚部abs（z）求复数z的模angle（z）求复数z的辐角，conj（z）求复数z的共轭复数exp（z）复数的指数函数，表示e^z如何用Matlab求集合的交集、并集、差集和补集　union(A,B)求集合A和B的并集intersect(A,B)求集合A和B的交集setdiff(A,B)求集合A和B的差集A-Bsetdiff(U,A)求集合A关于全集U的补集如何用matlab排序sort（v）　将向量v的元素从小到大排列（升序排列）sort(v,dim,’descend or ascend’)当dim=1时矩阵按列排序，descend or ascend用来控制升序还是降序当dim=2时矩阵按行排序，descend or ascend用来控制升序还是降序如何用Matlab求极限（1）极限：syms xlimit(f(x), x, a)求f(x)关于x趋于a时的极限（2）单侧极限：左极限：syms x limit(f(x), x, a,’left’)求f(x)关于x趋于a时的左极限右极限：syms x limit(f(x), x, a,’right’)求f(x)关于x趋于a时的右极限如何用Matlab求导数diff('f(x)')　diff('f(x)','x') 求f(x)关于x的导数或者：syms x diff(f(x))syms x diff(f(x), x)如何用Matlab求高阶导数如何用Matlab求高阶导数diff('f(x)'，n) diff('f(x)','x'，n)求f(x)关于x的n阶导数syms x diff(f(x)，n)syms x diff(f(x), x，n)如何用Matlab求不定积分int('f(x)') int ('f(x)','x')求f(x)关于x的积分syms x int(f(x))syms x int(f(x), x)如何用Matlab求定积分、广义积分int('f(x)',a,b) int ('f(x)','x',a,b)求f(x)关于x的积分，区间为a到b syms x int(f(x),a,b)syms x int(f(x), x,a,b)如何用Matlab展开级数syms x taylor(f(x), x, n,)a如何在Matlab中进行积分变换syms s tlaplace( f(t), t, s ) 拉普拉斯变换ilaplace( F(s), s, t ) 拉普拉斯变换的逆变换　syms t ωfourier( f(t), t, ω)傅立叶变换ifourier( F(ω), ω, t ) 傅立叶变换的逆变换　syms n zztrans( f(n), n, z) Z变换iztrans( F(z), z, n ) Z变换的逆变换　如何用Matlab解微分方程dsolve('微分方程'，'自变量')dsolve('微分方程'，'初始条件或边界条件'，'自变量') dsolve('D2x+2*x+x=sin(t)','x(0)=1','Dx(0)=1','t')如何用matlab求多变量函数的极限　以两个变量为例说明，多于两个变量的函数极限可以依次类推。

MATLAB数学实验

实验三圆周率的计算学号：姓名：XX一、实验目的1. 本实验涉及概率论、定积分、三角函数等有关知识，要求掌握计算π的三种方法及其原理。

2. 学习和掌握数学软件MATLAB 的使用方法。

二、实验内容圆周率是一个极其驰名的数。

从有文字记载的历史开始，这个数就引起了外行人和学者们的兴趣。

作为一个非常重要的常数，圆周率最早是出于解决有关圆的计算问题。

仅凭这一点，求出它的尽量准确的近似值，就是一个极其迫切的问题了。

事实也是如此，几千年来作为数学家们的奋斗目标，古今中外一代又一代数学家为此献出了自己的智慧和劳动。

回顾历史，人们对π的认识过程，反映了数学和计算技术发展情形的一个侧面。

π的研究，在一定程度上反映这个地区或时代的数学水平。

德国数学家康托说：“历史上一个国家所算的圆周率的准确程度，可以作为衡量这个国家当时数学发展水平的指标。

”直到19世纪初，求圆周率的值还是数学中的头号难题。

1. 圆周率的计算方法古人计算圆周率，一般是用割圆法。

即用圆的内接或外切多边形来逼近圆的周长。

Archomedes 用正96边形得到35位精度；刘徽用正3072边形得到5位精度；Ludolph V an Ceulen 用正2^62边形得到了35位精度。

这种基于几何的算法计算量大，速度慢，吃力不讨好。

随着数学的发展，数学家们在进行数学研究时有意无意得发现了许多计算圆周率的公式。

下面挑选一些经典的常用公式加以介绍。

除了这些经典公式外，还有很多其他公式和由这些经典公式衍生出来的公式，就不一一列举了。

1) Machin 公式2391a r c t a n451a r c t a n 16-=π ()121...753arctan 121753--++-+-=--n x x x x x x n n 这个公式由英国天文学教授John Machin 于1706年发现。

他利用这个公式计算到100位的圆周率。

Machin 公式每计算一项可以得到1.4位的十进制精度。

matlab实验九 π的近似计算

实验九π的近似计算【实验目的】1.了解圆周率π的计算历程。

2.了解计算π的割圆术、韦达公式、级数法、拉马努金公式、迭代法。

3.学习、掌握MATLAB软件有关的命令。

【实验内容】利用韦达（VieTa）公式2=π计算π的近似值。

【实验准备】1.割圆术2.韦达（VieTa）公式3.利用级数计算π4.拉马努金（Ranmaunujan）公式5.迭代方法6.π的两百位近似值【实验重点】1. 圆周率的计算历程【实验难点】1. 圆周率的各种计算公式【实验方法与步骤】练习1 用刘徽的迭代公式1162062626224,32,1n n n n x x s x x ++=--== 计算π的近似值。

相应的MATLAB 代码为>>clear;>>x=1;>>for i=1:30>>x=vpa (sqrt(2-sqrt(4-x^2)),15)%计算精度为15位有效数字 >>S=vpa(3*2^i*x,10)>>end计算可得x = .517638********* S = 3.105828541x = .261052384440103 S = 3.132628613x = .130806258460286 S = 3.139350203x = .654381656435522e-1 S = 3.141031951x = .327234632529735e-1 S = 3.141452472x = .163622792078742e-1 S = 3.141557608x = .818120805246955e-2 S = 3.141583892x = .409061258232818e-2 S = 3.141590463x = .204530736067660e-2 S = 3.141592106x = .102265381402739e-2 S = 3.141592517x = .511326923724832e-3 S = 3.141592619x = .255663463951308e-3 S = 3.141592645x = .127831732236766e-3 S = 3.141592651 x = .639158661510219e-4 S = 3.141592653 x = .319579330795908e-4 S = 3.141592653 x = .159789665403054e-4 S = 3.141592654 x = .798948327021645e-5 S = 3.141592654 x = .399474163511619e-5 S = 3.141592654 x = .199737081755909e-5 S = 3.141592654 x = .998685408779670e-6 S = 3.141592654 x = .499342704389851e-6 S = 3.141592654 x = .249671352194927e-6 S = 3.141592654 x = .124835676097464e-6 S = 3.141592654 x = .624178380487320e-7 S = 3.141592654 x = .312089190243660e-7 S = 3.141592654 x = .156044595121830e-7 S = 3.141592654 x = .780222975609150e-8 S = 3.141592654 x = .390111487804575e-8 S = 3.141592654 x = .195055743902288e-8 S = 3.141592654 x = .975278719511453e-9 S = 3.141592654练习2 用韦达公式2=π计算π的近似值。

matlab 圆周率

nl n 2 , 特别取针的长度 l d / 2 时， md m
本实验涉及的方法
3、蒙特卡罗法——蒲丰(Buffon)掷针法
x表示针的中点与最近的一条平行线间的距离, 为针与平行线的夹角.则0 x d 2 ,0 , 确定x 平面内一个矩形 .如图：
平面内一个矩形确定的夹角为针与平行线一条平行线间的距离表示针的中点与最近的针与平行线相交本实验涉及的方法3蒙特卡罗法蒲丰buffon掷针法针与平行线相交时针在平行线矩形区域内所组成的区域面积与矩形区域面积比即为针与直线相交的概率
实验三
圆周率的求法
问
题
• 设半径为1的圆的周长与圆的直径比为A; • 设半径为2的圆的周长与圆的直径比为B. • 请问A与B哪个大?
计算历程的几个阶段
3.分析方法时期 17世纪出现了数学分析，这锐利的工具使得许多初等数学束手无策的问题迎刃而解。 π 的计算历史也随之进入了一个新的阶段。
计算历程的几个阶段
1593年，韦达给出
1671年，苏格兰数学家詹姆斯.格雷戈里公开了他发现的公式： 3 5 7
arctan x x x3 x5 x7 (1 x 1)

故
2 nl md
本实验涉及的方法
3、蒙特卡罗法——蒲丰(Buffon)掷针法
• 随机生成一个针的中点与最近一条平行直线间的距离x,满足0≤x≤d/2; • 随机生成针与平行线的夹角a,满足0≤a≤ ； • 判断x≤L/2*sin[a],若成立，则针与直线相交，否则不相交。
本实验涉及的方法
计算历程的几个阶段
圆周长大于内接正四边形而小于外切正四边形
计算历程的几个阶段
在我国，首先是由数学家刘徽得出较精确的圆周率。他采用割圆术。

matlab指令大全

1、运算符：＋：加，－：减， *：乘， /：除， \：左除 ^：幂，‘：复数的共轭转置，（）：制定运算顺序。

2、常用函数表：sin( ) 正弦（变量为弧度）Cot( ) 余切（变量为弧度）sind( ) 正弦（变量为度数）Cotd( ) 余切（变量为度数）asin( ) 反正弦（返回弧度）acot( ) 反余切（返回弧度）Asind( ) 反正弦（返回度数）acotd( ) 反余切（返回度数）cos( ) 余弦（变量为弧度）exp( ) 指数cosd( ) 余弦（变量为度数）log( ) 对数acos( ) 余正弦（返回弧度）log10( ) 以10为底对数acosd( ) 余正弦（返回度数）sqrt( ) 开方tan( ) 正切（变量为弧度）realsqrt( ) 返回非负根tand( ) 正切（变量为度数）abs( ) 取绝对值atan( ) 反正切（返回弧度）angle( ) 返回复数的相位角atand( ) 反正切（返回度数）mod(x,y) 返回x/y的余数sum( ) 向量元素求和3、其余函数可以用help elfun和help specfun命令获得。

4、常用常数的值：pi 3.1415926…….realmin 最小浮点数，2^-1022i 虚数单位realmax 最大浮点数，（2－eps）2^1022j 虚数单位Inf 无限值eps 浮点相对经度＝2^-52NaN 空值1、!dir 可以查看当前工作目录的文件。

!dir& 可以在dos状态下查看。

2、who 可以查看当前工作空间变量名， whos 可以查看变量名细节。

3、功能键：功能键快捷键说明方向上键 Ctrl+P 返回前一行输入方向下键 Ctrl+N 返回下一行输入方向左键 Ctrl+B 光标向后移一个字符方向右键 Ctrl+F 光标向前移一个字符Ctrl+方向右键 Ctrl+R 光标向右移一个字符Ctrl+方向左键 Ctrl+L 光标向左移一个字符home Ctrl+A 光标移到行首End Ctrl+E 光标移到行尾Esc Ctrl+U 清除一行Del Ctrl+D 清除光标所在的字符Backspace Ctrl+H 删除光标前一个字符 Ctrl+K 删除到行尾Ctrl+C 中断正在执行的命令4、clc可以命令窗口显示的内容，但并不清除工作空间。

(完整)MatLab常用函数大全,推荐文档

求向量x的最小值函数是min(x)，用法与max(x)完全相同。
13、求矩阵的最大值和最小值
求矩阵A的最大值的函数有三种调用格式，分别是：
（1）max(A)：返回一个行向量，向量的i个元素是矩阵A的第i列的最大值。
（2）[y,u]=max(A)：返回行向量y和u，y纪录A的每列的最大值，u纪录每列最大值的行号。
factor(s)：对符号表达式s分解因式。
expand(s)：对符号表达式s进行展开。
例如：
syms x y;
s1=x^3-6*x^2+11*x-6
s1 =
x^3-6*x^2+11*x-6
factor(s1)
ans =
(x-1)*(x-2)*(x-3)
s2=(x-y)*(x+y)
s2 =
(x-y)*(x+y)
findsym(s)
ans =
x, y
findsym(5*x+2)
ans =
x
findsym(a*x+b*y+c)%符号变量c不会出现在结果中
ans =
a, b, x, y
29、符号表达式四则运算
符号表达式的加、减、乘、除和幂运算可分别由函数symadd、symsub、symmul、symdiv和sympow来实现。例如
对多项式求导数的函数是：
p=polyder(p1)：求多项式p1的导函数。
p=polyder(p1,p2)：求多项式p1和p2乘积的导函数。
[p,q]=polyder(p1,p2)：求多项式p1和p2之商的导函数，p、q是导函数的分子、分母。
例：求有理分式的导函数。
命令如下：
p1=[1,-1];

x=e^x用简单迭代法matlab

x=e^x用简单迭代法matlab篇一：在MATLAB中,我们可以使用简单迭代法来求解函数`f(x) = e^x`的导数`f"(x)`,即`f"(x) = e^x - 1`.下面是一个简单的迭代公式:```x0 = 0;for i = 1:1000x = x0 + (x - x0) / (i - 1);end```这个迭代公式使用了一个递归算法,每次将当前`x`值减去上一个`x`值,并除以迭代次数`i - 1`,直到迭代次数达到1000为止。

在每次迭代中,`x`值都会增加`e^x`的值,因此`f"(x)`的值也会增加`e^x - 1`的值。

我们可以通过计算`f"(x)`的值来得到`e^x`在`x`处的导数。

下面是MATLAB代码来求解`f"(x)`:```x = 0:0.1:10;f = e^x;f" = f";```在这个例子中,我们使用一个长度为10的`x`列表来求解`f"(x)`,并将结果保存到变量`f"`中。

使用这个迭代公式和`f"`的值,我们可以计算`e^x`在`x`处的导数,即`f(x)`的斜率。

这个斜率可以用来进行加速收敛的迭代计算,例如在优化问题中。

下面是MATLAB代码来计算`e^x`在`x`处的导数:```x = 0:0.1:10;f = e^x;f" = @(x) f"(x);```这个代码定义了一个函数`f"`,它使用迭代公式和`f`的值来计算`f"(x)`的值。

这个函数会使用一个循环来迭代`x`的值,并计算`f"(x)`的值。

这样,我们就可以使用这个函数来加速收敛,并得到更快的迭代结果。

总的来说,简单迭代法是一种常用的方法来求解函数的导数,特别是在数值计算中。

它简单易行,但收敛速度可能较慢,需要根据具体情况进行调整。

篇二：在本文中,我们将介绍如何使用简单迭代法在 MATLAB 中计算圆周率的近似值。

圆周率的计算方法及其MATLAB编程实现

辨
１Ｉ－．尊执术
应用研究
Hale Waihona Puke 圆周率的计算方法及其罗洪斌钱伟马梅玲
编程实现
（北方民族大学数学与信息科学学院宁夏银川７５００２１）
摘要：介绍圆周率的起源以及在各个时期的主要的近似计算方法的数学原理及代表人物，然后详细地给出古典方法、莱布尼兹（Ｌｅｉｂｎｉｚ）＆式、数值积分法，三种近似计算方法的数学原理及其推导过程，并结合现代数学计算软件ＭＡｎＡＢ编程实现算法过程，给出其运算过程及结果。并对三种计算方法进行具体分析，得到其相关结论。关键词：圆周率近似计算数学原理ＭＡＴＬＡＢ程序实现
ｓｙｍｓｋ
可推导出：９６ｓｉｎ０＜万＜９６ｔａｎ０ＭＡＴＬＡＢ程序实现：
后得出：
７ｒ
三４一ｔａｎ薹（、）２ｋ１
化简后可得
丌：４（一）１
（２
．
３）
ｓｉｎＯ＜＜ｔａｎＯ，０
逐渐增９，ｋ的值，圆周率的精确度会越来越高。ＭＡＴＬＡＢ程序实现：
中图分类号：ＴＰ３１２文献标识码：Ａ文章编号：１００７．９４１６（２０１６）０８ — ００５７．０２
Ｉｏｒｍａｔｌｏｎｇ第十六位希腊字母，欧拉于１８世纪提倡用丌代表圆周率。圆周率最早被发现于一块古巴比伦的石板上，此时的丌是一个常数。ｘ＝ｔａｎ（ｐｉ／９６）ｙ＝９６＊ｘ在《几何原本中７ｒ以一个小数被记载，这是圆周率最早的书面记载。此时圆周率具体的定义被给出，即圆周长与其直径之比，我国的得：９６ｔａｎ０＝３．１４２７１４５９９６４５３６８古算经周髀算经也有记载。为了得到圆周率更加精确的近似值，得：３．１４１０３１９５０８９０５０９＜ｒ￣＜３．１４２７１４５９９６４５３６８１．２莱布ｆ￣＇（Ｌｅｉｂｎｉｚ）＆式从古至今，一位又一位数学家奉献出了自己毕生的精力与心血。由于数学计算的发展，圆周率的近似计算主要经历实验时期、几何法用级数法求圆周率在分析法时期十分流行，这个时期产生了莱时期、分析法时期、三大时期。下面具体介绍三个时期圆周率计算方布尼兹公式、马琴公式、傅里叶级数等计算方法，但主要介绍莱布尼法。兹公式。莱布尼兹公式主要用幂级数的推到而得到的，其原理如下：实验时期：首先，函数的麦克劳林（ｃｏｌｉｎＭａｄａｕｒｉｎ）展开式为由于文化发展滞后，科技有局限性。因此此时的圆周率都是通过测量得出，常用：３来代替圆周率，最早见于欧几里得的（（几何原１一＋一… ＋ｆ一１１＋．．．＋１本，印度的百道梵书》也有记载。几何法时期： ∑（一１），一ｌ＜ｘ＜１这个时期数学家们用严谨的科学方法来计算圆周率，在古埃及将上式中的用替换得草书（约公元前１７００年）中提到石＝（４／３）３．１６０４。先驱者阿基米德通过计算数学理论算的圆周率。他从单位圆开始，通过迭代算法和ｌ－Ｘ２＋Ｘ４－－＋ｃ－１＋．－．两侧数值逼近算得（３＋（１０／７１））＜丌＜（３＋（１／７）），这种方法被命名为古典方法。我国的刘徽利用“ 割圆术 ” 求出了更加精确的圆周率，（－１／，－ｌ＜ｘ＜ｌ（２．・１ ” 算出万＝３．１４１６。继刘徽之后，祖冲之３．１４１５９２６＜丁ｒ＜３．１４１５９２７，将的精确度提高到８位，保持世界记录九百多年［１】。将（２．１）式两边同时积分可得：分析法时期：这随着各种方法理论增多，人们开始使用无穷级数或无穷连乘积来计算７ｒ。自１５９３年韦达公式面世之后，沃利斯、ＪｏｈｎＭａｃｈｉｎ、ＳｒｉｎｉｖａｓａＲａｍａｎｕｊａｎ、ＦａｂｒｉｃｅＢｅ￣ａｒｂ等先后给出了更好的计算ｆ（）＿１），－ｌ＜ｘ＜ｌ（２．２）方法。他们将圆周率的精确度一次又一次的提升到新的高度。

1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。

接着有多种表达式出现。如沃利斯1650 年给出：

22446688 2 133 4 55 7 7
1706年，英国天文学教授John Machin 利用
x x x n1 x arctan x x 1 3 5 7 2n 1
发现了下面的公式
1 1 16 arctan 4 arctan , 5 239

1
0
e
x2
dx
抛物线方法 y=inline('exp(-x.^2)'); quad(y,0,1) ans = 0.746826
梯形方法
条件： f ( x)在[a, b]上连续或有有限个第一类间断点 . 公式推导： b f ( x )dx s1 s2 sn
aHale Waihona Puke yax1s2
当区间划分为n(n>1)等分时
o
a
x1
x2 x3
x4
x5
b x

b a
n 1 n x k 1 x k h f ( x )dx S n ( f ( x0 ) 2 f ( xk ) f ( xn ) 4 f ( )) , 6 2 k 1 k 1 ba h , xk a kh k 0,1,2, , n n
分析法时期
• 这一时期人们开始摆脱求多边形周长的繁难计算，利用无穷级数或无穷连乘积来算 π 。 • 1593年，韦达给出
2 2 2 2 2 2 2 2 2 2
这一不寻常的公式是 π 的最早分析表达式。甚至在今天，这个公式的优美也会令我们赞叹不已。它表明仅仅借助数字2，通过一系列的加、乘、除和开平方就可算出 π 值。
将积分区间 [0， 1]分成n等份，在每一个小区间上，选取中点为 i , 编写下面的程序计算的近似值.
n=50; %定义等分积分区间数，可以更改 i=0:1/n:1; s=0; for k=1:length(i)-1 s=s+(1/(1+((i(k)+i(k+1))/2)^2))*1/n; end 4*s
1 1 32 1 256 64 n 0 1024 4n 1 4n 3 10n 1
n
64 4 4 1 . 10n 3 10n 5 10n 7 10n 9
从而，大大降低了圆周率近似值的计算量.
No Yes
clear;a=0;b=1; f=inline('4/(1+x*x)'); t1=(b-a)/2*(f(a)+f(b)); er=1;n=1; while er>1.0e-6 h=(b-a)/n; s=0; for i=1:n s=s+f(a+i*h-h/2); end t2=(t1+h*s)/2; er=abs(t2-t1);
2.圆周率的数值积分计算方法 1 1 1 1 0 1 x 2 dx 4 40 1 x 2 dx
表1给出的是不同等分区间数计算出的圆周率的近似值.
等分区间数n 10 20 50 100 500 1000 5000
圆周率的近似值
3.14242598500110 3.14180098689309 3.14162598692300 3.14160098692312 3.14159298692312 3.14159273692313 3.14159265692313
1 4 2 1 1 n . 8n 4 8n 5 8n 6 n 0 16 8n 1

该公式的最大优点在于：经后来人将该公式变形后打破了传统的计算方法，可以直接计算圆周率的任意第n位数，而不是先计算前面的n-1位数．
1997 年， Fabrice Bellard 发表了一个比 BBP 算法更快的公式
s3
s4
f ( x i 1 ) f ( x i ) ( ( xi xi 1 )) o 2 i 1
n
s1
x2 x3
b x
当区间划分为n等分时
b
n 1 h f ( xi ) f (b)) , a f ( x )dx Tn 2 ( f (a ) 2 i 1 ba 其中 h , xk a kh k 1,2,, n 1 n
fprintf('t=%.6f,r=%.6f\n',t2,er);
输出结果：STOP
n=2*n; t1=t2;
end
条件：分析：
抛物线方法
y
s4
f ( x )在[a, b]上连续或有有限个第一类间断点 .
梯形法是对每个子区间用梯形面积近似曲线下面积累加而成；而抛物线法是对每个子区间考虑其中点，用三点决定的抛物线下面积来近似.
举例利用复化梯形算法求Pi的近似值.
y
取
4
运用复化梯形算法
ba T1 ( f (a ) f (b)) 2
1
0
1
x
n 4 1 h f (a ih )) 0 1 x 2 dx T2n 2 (Tn h 2 i 1
输入初值：
n 1 h T2 n (Tn h f (a ih )) 2 2 i 1
1989年，David 和 Gregory Chudnovsky 发表了下面的公式
1 ( 1)n (6n)! 13591409 545140134n 12 , 3 3 3 n n 0 ( 3n)!( n! ) 640320 2
并在1994年计算到了4044000000位．它的另一种形式是
——trapz(x,y)

b a
ba f ( x )dx T1 ( f (a ) f (b)) 2
复化梯形方法
y
a
x1
当区间等距划分为n个子区间时

b a
n 1 h f ( x )dx Tn ( f (a ) 2 f ( xi ) f (bo )) , 2 i 1 ba h , xk a kh k 0,1,2, , n n
.
1 1 ( 1) 4 arctan 1 4(1 ). 3 5 2n 1
n 1
1 1 ( 1) 4 arctan 1 4(1 ). 3 5 2n 1
编写下面的程序: n=10; %选择展开式的次数 s=0; digits(22); %定义计算过程中的精度 for k=1:n s=s+4*(-1)^(k+1)/(2*k-1); end vpa(s,20) %定义显示精度为20位
426880 10005 . (6n)!(545140134 n 13591409) 3 3n ( n ! ) ( 3 n )! ( 640320 ) n 0
1995 年，由 David Bailey,Peter Borwein 和 Simon Plouffe 共同发表了下面的圆周率计算公式（简称BBP公式）
1.圆周率的幂级数计算方法
例1
利用arctan x的Maclaurin展开式，计算的近似值.
1 2 4 1 x x 2 1 x
3 5
解
(1) n 1 x 2 n
n 1
.
x x ( 1) x arctan x x 3 5 2n 1
2 n1
数值积分简介
在高等数学中有一类积不出的积分，如

1
0
e
x2
dx

2
1 1 x
4
0
dx

（概率积分）
（椭圆积分）
抛物线法利用数值积分法梯形法
MATLAB命令 •梯形方法 ——trapz(x,y) •抛物线方法——quad(f,a,b) 如求积分的近似值梯形方法输入： x=0:0.1:1; y= exp(-x.^2); trapz(x,y) 输出： ans = 0.746211
在中国
• 祖冲之: 在刘徽研究的基础上，进一步地发展，经过既漫长又烦琐的计算，一直算到圆内接正 24576边形，而得到一个结论： • 3.1415926 ＜ π ＜ 3.1415927 同时得到π 的两个近似分数：约率为22／7；密率为355／113。
• 他算出的 π 的8位可靠数字，不但在当时是最精密的圆周率，而且保持世界记录九百多年。以致于有数学史家提议将这一结果命名为“祖率”。
实验时期
• 基于对一个圆的周长和直径的实际测量而得出的。 • 在古代世界，实际上长期使用 π ＝3这个数值。 • 最早见于文字记载的有基督教《圣经》中的章节，其上取圆周率为3。这一段描述的事大约发生在公元前950年前后。
几何法时期

真正使圆周率计算建立在科学的基础上，首先应归功于阿基米德。他是科学地研究这一常数的第一个人，是他首先提出了一种能够借助数学过程而不是通过测量的、能够把 π 的圆周长大于内接正多边值精确到任意精度的方法。形周长而小于外切正多边由此，开创了圆周率计算形周长．的第二阶段。据说阿基米德用到了正 96边形才算出他的值域。
o
a
x1s2
s3
s1
x2 x3
b x
y0
y1
y2
S 3
x3 x2 x x3 ( f ( x2 ) 4 f ( 2 ) f ( x3 )) 6 2
h/ 2
h/ 2
h S ( y0 4 y1 y2 ) 6
y
S 3 x x3 h ( f ( x2 ) 4 f ( 2 ) f ( x3 )) 6 2
x2 x3
x4