Matlab计算圆周率
拉格朗日插值多项式积分求圆周率近似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 圆周率
本实验涉及的方法
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常用函数大全,推荐文档
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];
- 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