西安交大数学实验五matlab计算圆周率pi-精选文档
2021年数学实验报告圆周率的计算mathematica
数学试验汇报试验序号: 2 日期: 月日试验结果汇报及试验总结:一、数值积分法计算π因为单位圆半径为1, 它面积等于π, 所以只要计算出单位圆面积, 就算出了π。
在坐标轴上画出以圆点为圆心, 以1为半径单位圆, 则这个单位圆在第一象限部分是一个扇形, 而且面积是单位圆1/4, 于是, 我们只要算出此扇形面积, 便能够计算出π。
而且单位精度可能会影响计算结果, 下面将给出不一样n计算所得结果并讨论差异。
1.当n=1000时命令:n=1000;y[x_]:=4/(1+x*x);s1=(Sum[y[k/n],{k,1,n-1}]+(y[0]+y[1])/2)/n;s2=(y[0]+y[1]+2*Sum[y[k/n],{k,1,n-1}]+4*Sum[y[(k-1/2)/n],{k,1,n}])/( 6*n);Print[{N[s1,20],N[s2,30],N[Pi,30]}];结果以下:2.当n=5000时命令:n=5000;y[x_]:=4/(1+x*x);s1=(Sum[y[k/n],{k,1,n-1}]+(y[0]+y[1])/2)/n;s2=(y[0]+y[1]+2*Sum[y[k/n],{k,1,n-1}]+4*Sum[y[(k-1/2)/n],{k,1,n}]) /(6*n);Print[{N[s1,20],N[s2,30],N[Pi,30]}];运行结果:3.当n=10000时命令:n=10000;y[x_]:=4/(1+x*x);s1=(Sum[y[k/n],{k,1,n-1}]+(y[0]+y[1])/2)/n;s2=(y[0]+y[1]+2*Sum[y[k/n],{k,1,n-1}]+4*Sum[y[(k-1/2)/n],{k,1,n}])/( 6*n);Print[{N[s1,20],N[s2,30],N[Pi,30]}];Plot[{4(1-x*x)},{x,0,1}]运行结果:4. 结果分析: 当数值积分法得到 近似值为3.338328,能够看出, 用这种方法计算所得到 值是相当正确, n 越大, 计算出来扇形面积近似值就越靠近 正确值。
圆周率的近似计算
西安交通大学数学实验报告实验题目：圆周率的近似计算填写报告日期：2013年6月3日一、实验目的1.学会用MATLAB 软件对圆周率进行数值积分计算，求出近似值；2.学会用蒙特卡罗法求解几何体体积；二、实验任务1.对P89页示例4，利用矩形、和抛物线形方法求近似值，并与书上的方法进行比较，看哪种方法收敛速度快？2.完成P92页练习5第2题.三、问题分析与求解分析：（1）梯形公式：依次取每个极小区间中点为i ξ，利用积分式∑⎰=→=+ni i f dx x xi 10max102)(11lim ξ i x ∆计算圆周率π的近似值 （2）矩形公式：依次取每个区间某一端点i x ，利用计算圆周率π的近似值。
(3)抛物线公式：利用计算圆周率π的近似值。
1()d ()Δ,nbi i ai f x x f x x =≈∑⎰2221()d ()d i i nbx ax i f x x f x x-==∑⎰⎰2. 蒙特卡罗法求解冰淇淋体积：在区域{(x,y,z)｜x 2 + y 2 ≤ 1,0 ≤x,0≤y, 0≤ z ≤2}内随机取点，对落在冰淇淋体第一卦限内的点进行计数，则落在冰淇淋内的点数m 与落在长方体区域内的点数n 的比值等于体积的比值.故V （冰淇淋）= nm8。
（1）圆周率的数值积分计算矩形方法程序： n=5000 i=0:1/n:1; s=0;for k=1:length(i)s=s+(1/(1+i(k)^2))*1/n; end 4*s 运行结果：抛物线方法程序：n=5000;i=0:1/n:1;s=0;for k=2:length(i)/2s=s+((1/(1+i(2*k-2)^2))+(4/(1+i(2*k-1)^2))+(1/(1+i(2*k)^2))); end8/(6*n)*s运行结果：与书中方法对比，易得抛物线形公式收敛速度最快。
（2）蒙特卡罗法求冰淇淋体积程序：cs=0n=500000for i=1:na=rand(1,3);if a(1)^2+a(2)^2<=4*a(3)-(2*a(3))^2&a(1)^2+a(2)^2<=(2*a(3))^2;cs=cs+1;endendV=8*cs/n运行结果：评论：使用蒙特卡罗法计算冰淇淋的体积，存在较大的随机性，实验结果与理论值相差较大。
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 龙贝格求圆周率
matlab 龙贝格求圆周率龙贝格求圆周率是一种计算圆周率的方法，它是一种数值积分方法，采用复合梯形法和复合辛普森法相结合的方法求圆周率。
该方法通过将积分区间分段，采用递归算法求解，从而获得更高的精确度。
该方法常常用于数值计算和科学计算中，能够快速准确地得到圆周率的近似值。
在Matlab中，可以通过编写相应的程序实现龙贝格求圆周率的计算。
Matlab中提供了许多数值计算函数，可以方便地实现该方法。
以下是一个计算圆周率的示例程序：% 龙贝格求圆周率clear all;clc;format long;f = @(x) 4 ./ (1 + x.^2); % 定义被积函数a = 0;b = 1; % 积分区间tol = 1e-10; % 精度要求H = b - a;% 初始化T(1,1) = H / 2 * (f(a) + f(b));n = 1;% 迭代计算while 1n = n + 1;H = H / 2;x = a + H : 2*H : b;T(n,1) = 1/2 * T(n-1,1) + H * sum(f(x)); % 复合梯形公式 for m = 2 : nk = 2^(2*(m-1)); % 计算k值T(n,m) = (k*T(n,m-1) - T(n-1,m-1)) / (k-1); % 复合辛普森公式endif abs(T(n,n-1) - T(n,n)) < tol % 迭代停止条件break;endend% 输出结果pi_approx = T(n,n);fprintf('圆周率的近似值为：%f\n', pi_approx);请注意，本文中不含任何网址、超链接和电话信息，仅介绍数学计算方法和Matlab编程技巧。
如何计算π的值(MATLAB)
如何计算π的值1、蒙特卡罗(Monte Carlo)法思想:取一正方形A,以A 的一个顶点为圆心,A 的边长为半径画圆,取四分之一圆(正方形内的四分之一圆)为扇形B 。
已知A 的面积,只要求出B 的面积与A 的面积之比B AS k S =,就能得出B S ,再由B 的面积为圆面积的四分之一,利用公式2=S R π圆即可求出π的值。
因此,我们的目的就就是要找出k 的值。
可以把A 与B 瞧成就是由无限多个点组成,而B 内的所有点都在A 内。
随机产生n 个点,若落在B 内的有m 个点(假定A 的边长为1,以扇形圆心为坐标系原点。
则只要使随机产生横纵坐标x 、y 满足221x y +≤的点,就就是落在B 内的点),则可近似得出k 的值,即m k n=,由此就可以求出π的值。
程序(1):i=1;m=0;n=1000;for i=1:na=rand(1,2);if a(1)^2+a(2)^2<=1m=m+1;endendp=vpa(4*m/n,30)程序运行结果:p =3、14000002、泰勒级数法思想:反正切函数arctan x 的泰勒级数展开式为:35211arctan (1)3521k k x x x x x k --=-+-⋅⋅⋅+-+⋅⋅⋅- 将1x =代入上式有1111arctan11(1)43521n n π-==-+-⋅⋅⋅+--、 利用这个式子就可以求出π的值了。
程序(2):i=1;n=1000;s=0;for i=1:ns=s+(-1)^(i-1)/(2*i-1);endp=vpa(4*s,30)程序运行结果:p =3、1494126当取n 的值为10000时,就会花费很长时间,而且精度也不就是很高。
原因就是1x =时,arctan1的展开式收敛太慢。
因此就需要找出一个x 使得arctan x 收敛更快。
若取12x =,则我们只有找出α与4π的关系,才能求出π的值。
令1arctan 2α=,4πβα=-, 根据公式tan tan tan()1tan tan αβαβαβ++=-有1tan 3β=,则有11arctan arctan 423π=+。
21——matlab——用蒙特卡罗方法计算pi的值
用蒙特卡罗( Monte Carlo ) 投点法计算 的值程序：n=10000;m=0; a=2;for i=1:nx=rand(1)*a/2;y=rand(1)*a/2;if x^2+y^2<=(a/2)^2m=m+1;endendfprintf('蒙特卡罗投点法计算的pi为：%f\n',4*m/n)结果：蒙特卡罗投点法计算的pi为：3.117200蒙特卡罗投点法计算的pi为：3.132400蒙特卡罗投点法计算的pi为：3.165600蒙特卡罗投点法计算的pi为：3.135200蒙特卡罗投点法计算的pi为：3.146800蒙特卡罗投点法计算的pi为：3.140400蒙特卡罗投点法计算的pi为：3.114800蒙特卡罗投点法计算的pi为：3.117200 蒙特卡罗投点法计算的pi为：3.154800 蒙特卡罗投点法计算的pi为：3.140400 蒙特卡罗投点法计算的pi为：3.121200 蒙特卡罗投点法计算的pi为：3.134400 蒙特卡罗投点法计算的pi为：3.122800 蒙特卡罗投点法计算的pi为：3.127600 蒙特卡罗投点法计算的pi为：3.147200 蒙特卡罗投点法计算的pi为：3.145200 蒙特卡罗投点法计算的pi为：3.158400 蒙特卡罗投点法计算的pi为：3.147600 蒙特卡罗投点法计算的pi为：3.147600 蒙特卡罗投点法计算的pi为：3.146400 蒙特卡罗投点法计算的pi为：3.112000 蒙特卡罗投点法计算的pi为：3.180000 蒙特卡罗投点法计算的pi为：3.120000 蒙特卡罗投点法计算的pi为：3.153600 蒙特卡罗投点法计算的pi为：3.144000 蒙特卡罗投点法计算的pi为：3.148000 >>。
matlab使用级数求pi的值的程序
[matlab使用级数求pi的值的程序]在数学和计算机科学领域，级数是一种非常重要的概念，它常常被用来进行数值计算和数学建模。
其中，π（pi）的计算就是一个著名的例子。
π是一个无理数，其精确值无法用有限的小数、分数或代数式来表示。
人们常常使用级数来逼近π的值。
在这篇文章中，我将和大家共享如何利用matlab来使用级数求π的值的程序，深入讨论该程序的原理和实现方式，并回顾整个计算过程，以便读者更深入地理解这一数学计算的背后原理。
1. 原理和思路在计算π的值时，可以利用数学中的级数公式来逐步逼近π的值。
其中，有一种著名的级数公式就是莱布尼茨级数公式，即：π/4 = 1 - 1/3 + 1/5 - 1/7 + 1/9 - ...根据该级数公式，我们可以不断地增加级数的项数，来逼近π/4的值，最终得出π的近似值。
2. 算法和程序接下来，我们将利用matlab来实现这个级数求π的值的程序。
我们可以定义一个变量sum来表示级数的和，然后利用一个循环来不断更新sum的值。
在每一轮循环中，我们可以根据级数的奇偶性来确定每一项的正负号，并将其加到sum中。
具体的matlab代码如下：```matlabsum = 0;n = 1;precision = 1e-6;term = 1;while abs(term) > precisionterm = (-1)^(n-1) * 1/(2*n-1);sum = sum + term;n = n + 1;endpi_value = 4 * sum;disp(['计算得到的π的近似值为：', num2str(pi_value)]);```在上述程序中，我们设定了一个精度precision，当每一项的绝对值小于precision时，我们认为级数已经收敛，可以结束循环。
我们将得到的sum乘以4，得到π的近似值。
3. 实际运行和结果当我们在matlab中运行上述程序时，会得到一个近似值为3.141591的π。
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数学实验报告
数学实验报告院系：西安交通大学软件学院软件工程系；班级：软件11；项目：MATLAB软件与基础数学实验；指导教师：张芳；日期：2012年6月11日星期一；学生姓名：贺翔；学号：2111601006；题目【一】在同一坐标系下画出函数y=sin x, y=cos x, y=0.2e0.1x sin (0.5x)和y=0.2e0.1x cos(0.5x)在区间[0，2pi]的曲线图，并对该图进行修饰。
(1)解题思路：首先按步长赋值法生成x向量，则生成相应函数值向量；然后运用plot命令，再添加网格或者其他修饰等。
(2)算法设计：x=0:0.07*pi:2*pi;y1=sin(x);y2=cos(x);y3=0.2.*exp(0.1.*x).*sin(0.5.*x);y4=0.2.*exp(0.1.*x).*cos(0.5.*x);plot(x,y1,'r--',x,y2,'k:',x,y3,'g.',x,y4,'b+','linewidth',3,'markersize',5); grid;xlabel('variable\it{x}')ylabel('variable\it{y}')title('four cruves')text(2.6,0.7,'sin(x)')text(3.5,0.3,'0.2.*exp(0.1.*x).*sin(0.5.*x)')text(5.8,0.8,'cos(x)')text(4.1,-0.4,'0.2.*exp(0.1.*x).*cos(0.5.*x)')(3)结果截图：题目【二】某农夫有一个半径10m的圆形牛栏，长满了草。
他要将一头牛拴在牛栏边界的栏桩上，但只让牛吃到一半草，问栓牛鼻的绳子应为多长？(1)解题思路：设R 为牛栏的半径，而栓牛绳长为r; 则根据数学公式：S=12R 2·4arcsin(r 2R )+ 12r 2·2arccos(r 2R )-2×12r √R 2−r 24；以及令S=12πR 2，即可解出方程的解。
测量圆周率
if C(i)<1,cnt=cnt+1; end
end
利用面积比 =P= ,
可以知道pi值为
mid_pi(j)=cnt/numel(A)*4;
%随机数每次测量10次求平均:结束%
res_pi=vpa(sum(sum(mid_pi))/10,7)
1,随机事件发生的概率以计算圆周率pi
.m文件:cau_pi1(num)
定义随机矩阵行数变量num,
%随机数每次测量10次求平均:开始%
for j=1:10;
产;2*rand(num,400);B=-1+2*rand(num,400);
然后再以A,B中矩阵的对应位置为x,y坐标,依次判断产生的随机数是否在意(0,0)为圆心,1为半径的圆内,如果再圆内,则计数相加
由于num参数不确定,然后统计次数少了不具代表性,可以让num递增来观察效果,写成函数:
cau_pi2()
开头:for num=250:250:20000;
结尾:t=num/250;%num变量数
per_pi(t)=vpa((res_pi(t)-3.14159)/3.141592,6);%num对误差影响
西安交通大学数学实验报告(用MATLAB绘制二维、三维图形)
实验报告（二）完成人：L.W.Yohann注：本次实验主要学习了用MATLAB绘制二维、三维图形的基本命令、图形的标识与修饰以及用符号函数绘图，在学习完成后小组对52页的上机练习题进行了程序编辑和运行。
1.绘制数列变化趋势图.解：在编辑窗口输入：n=1:100;an=(1+1./n).^n;plot(n,an,'r*')grid并保存，命名为lab1；在命令窗口中输入lab1,得：2.绘制数列变化趋势图.解：在编辑窗口输入：n=1:0.1:50;an=n.^(1./n);plot(n,an,'r*')grid并保存，命名为lab2；在命令窗口中输入lab2,得：3.绘制函数在无定义点处的变化趋势.解：在编辑窗口输入：x=-10:0.05:10;y=sin(x)./x;plot(x,y,'r*')grid并保存，命名为lab3；在命令窗口中输入lab3,得：4.在同一坐标系中画出函数及其Taylor多项式的图像解：y=sinx在编辑窗口输入：syms xf=sin(x);T6=taylor(f,x);T8=taylor(f,x,'Order',8);T10=taylor(f,x,'Order',10);T12=taylor(f,x,'Order',12);fplot([T6 T8 T10 T12 f])xlim([-8 8])grid onlegend('approximation of sin(x) up to O(x^6)',...'approximation of sin(x) up to O(x^8)',...'approximation of sin(x) up to O(x^{10})',...'approximation of sin(x) up to O(x^{12})',...'sin(x)','Location','Best')title('Taylor Series Expansion')并保存，命名为lab4sin；在命令窗口中输入lab4sin,得：y=exp(x)在编辑窗口输入：syms xf=exp(x);T6=taylor(f,x);T8=taylor(f,x,'Order',8);T10=taylor(f,x,'Order',10);T12=taylor(f,x,'Order',12);fplot([T6 T8 T10 T12 f])xlim([-8 8])grid onlegend('approximation of exp(x) up to o(x^6)',...'approximation of exp(x) up to o(x^8)',...'approximation of exp(x) up to o(x^{10})',...'approximation of exp(x) up to o(x^{12})',...'exp(x)','Location','Best')title('Taylor Series Expansion')并保存，命名为lab4exp；在命令窗口中输入lab4exp,得：5.符号函数绘图.注：在matlab r2010b 和matlab r2019b中对绘制函数图像的输入方法有不同的要求，故此类题分两个版本来求解。
matlab 圆周率的近似计算 实验报告
开放性数学实验报告（2016 / 2017学年第 2学期）题目：基于MATLAB的圆周率近似计算专业通信工程学生姓名杨坤冯著豪周李鑫班级学号 B******** B******** B********指导教师赵礼峰指导单位南京邮电大学理学院日期 2017/5/20MATLAB圆周率的近似计算B16011115 杨坤 B16011110 冯著豪 B16011124周李鑫摘要：圆周率（Pi）是圆的周长与直径的比值，一般用希腊字母π表示，是一个在数学及物理学中普遍存在的数学常数。
π也等于圆形之面积与半径平方之比。
是精确计算圆周长、圆面积、球体积等几何形状的关键值。
在分析学里，π可以严格地定义为满足sin x = 0的最小正实数x。
计算圆周率一直是很多人的追求。
在电子计算机还没有发明的时候就有很多先贤用各种方法计算了圆周率的近似值最著名的应该是祖冲之，他计算出了圆周率的位数达到了小数点后七位。
该记录在世界范围内保持了八百年。
之后圆周率的计算进入了分析法时期，这一时期人们开始利用无穷级数或无穷连乘积求π，摆脱可割圆术的繁复计算。
无穷乘积式、无穷连分数、无穷级数等各种π值表达式纷纷出现，使得π值计算精度迅速增加。
在分析法的基础上，电子计算机的出现使得圆周率的计算精度大幅提高。
计算圆周率已经成为评判超级计算机的性能指标的项目之一。
如今个人计算机的性能也达到了一个极高的程度。
学习使用计算机计算圆周率可以帮助我们更好地学习matlab同时对数学也会有更深的理解。
关键词：圆周率计算；投点法；定积分计分法；幂级数；韦达公式一、问题分析计算圆周率有很多方法，不同方法之间自然也有好坏之分。
在强大的计算机性能的支持下，我们能使用不同的方法计算圆周率并且感受不同方法孰优孰劣。
首先我们需要了解不同的计算方法是怎么计算圆周率的，然后使用matlab编写代码帮助我们实现算法，计算出圆周率。
二、实验方法1.投点法：投点法，顾名思义就是通过投点计算圆周率。
西安交通大学数学实验报告(MATLAB求解开普勒方程和方程求根)
实验报告（五）完成人：L.W.Yohann注：本次实验主要学习了用MATLAB求解开普勒方程和方程求根的问题，了解学习了用fzero命令、二分法、Newton迭代法、一般迭代法求解方程，以及学习了非线性方程组的求解问题，完成后，小组对第90页的上级练习题进行了程序编辑和运行。
1.绘图并观察函数零点的分布.解：在编辑窗口输入：f=inline('x-0.5*sin(x)-1');fplot(f,[0,1])grid存盘后运行得2. 利用fzero 命令求解方程.解：在编辑窗口输入：f=inline('x-0.5*sin(x)-1');c=fzero(f,[1,2])存盘后运行得c =1.49873. 用二分法求解方程.求解（1）方程x^2-2=0在（0，2）内的近似根；（2）圆x^2+y^2=2与曲线y=e^-x 的两个交点；（3）方程∫t 21+t 2x 0dt =12的近似根. （1）解：在编辑窗口输入：00.10.20.30.40.50.60.70.80.91-1-0.9-0.8-0.7-0.6-0.5-0.4-0.3f=inline('x^2-2');x1=0;x2=2;while abs(x1-x2)>10^(-5)x3=(x1+x2)/2;if f(x3)==0break;elseif f(x1)*f(x3)>0x1=x3;else f(x2)*f(x3)>0;x2=x3;endendx0=x3存盘后运行得x0 =1.4142（2）解：在编辑窗口输入：f=inline('(x^2)*exp(2*x)+1-2*exp(2*x)');x1=0;x2=2;x5=-2;x6=0;while abs(x1-x2)>10^(-5)x3=(x1+x2)/2;if f(x3)==0break;elseif f(x1)*f(x3)>0x1=x3;else f(x2)*f(x3)>0;x2=x3;endendwhile abs(x5-x6)>10^(-5)x7=(x5+x6)/2;if f(x7)==0break;elseif f(x5)*f(x7)>0x5=x7;else f(x6)*f(x7)>0;x6=x7;endendx0=x3x4=x7存盘后运行得x0 =1.3922x4 =-0.3203（3）解：在编辑窗口输入：clear;clc;syms t xf1=(t^2)/(1+t^2);f2=int(f1,t,0,x);%¼ÆËã²»¶¨»ý·Öf=inline('x - atan(x)-0.5');x1=-5;x2=5;while abs(x1-x2)>10^(-5)x3=(x1+x2)/2;if f(x3)==0break;elseif f(x1)*f(x3)>0x1=x3;else f(x2)*f(x3)>0;x2=x3;endendx0=x3存盘后运行得x0 =1.47504.用Newton迭代法求解方程求解：x=0.5sinx+1的近似根；解：在编辑窗口输入：f=inline('x-0.5*sin(x)-1');df=inline('1-0.5*cos(x)');d2f=inline('0.5*sin(x)');a=1;b=2;dlt=1.0e-5;if f(a)*d2f(a)>0x0=a;elsex0=b;endm=min(abs(df(a)),abs(df(b)));k=0;while abs(f(x0))>m*dltk=k+1;x1=x0-f(x0)/df(x0);x0=x1;fprintf('k=%d x=%.5f\n',k,x0); end存盘后运行得k=1 x=1.54858k=2 x=1.49933k=3 x=1.498705.求解非线性方程组.试求非线性方程组{2x12−x1x2−5x1+1=0x1+3lgx1−x22=0的解，初值如下：（1）x0=[1.4,−1.5]（2）x0=[3.7,2.7]解：在编辑窗口输入：function f=group5(x)f=[2*x(1)^2-x(1)*x(2)-5*x(1)+1;x(1)+3*log10(x(1))-x(2)^2];(1)：输入：[f,fval]=fsolve('group2',[1.4,-1.5]) 运行得f =1.4589 -1.3968fval =1.0e-011 *0.0759-0.6178(2)：输入：[f,fval]=fsolve('group2',[3.7,2.7])运行得f =3.4874 2.2616fval =1.0e-006 *0.0059-0.20126.解决实际问题.为了在海岛I与某城市C之间铺设一条地下光缆，每千米光缆铺设成本在水下部分使C1万元，在地下部分使C2万元，为使得该光缆的总成本最低，光缆的转折点P（海岸线上）应该取在何处？如果实际测得海岛I与城市C之间的水平距离l=30km,海岛距海岸线垂直距离h1=15km，城市距海岸线垂直距离h=10km,C1=3000万元/km,C2=1500万元/km,求P点的坐标（误差<10−3km）.解：在编辑窗口输入：f=inline('(3000*x)/(x^2 + 225)^(1/2) + (750*(2*x - 60))/((x - 30)^2 + 100)^(1/2)'); x1=5;x2=10;while abs(x1-x2)>10^(-3)x3=(x1+x2)/2;if f(x3)==0break;elseif f(x1)*f(x3)>0x1=x3;else f(x2)*f(x3)>0;x2=x3;endendfprintf('x=%.5f',x3) 存盘后运行得x=7.69104>>。
圆周率的计算-MATLAB
n0
1 16n
4 8n
1
2 8n
4
1 8n
5
1 8n
6
.
该公式的最大优点在于：经后来人将该公式变 形后打破了传统的计算方法，可以直接计算圆周率 的任意第n位数，而不是先计算前面的n-1位数．
1997年，Fabrice Bellard发表了一个比BBP算 法更快的公式
1
1 n
6边形
12边形
24边形
圆... 3.142857..., (n 96)
71
7
刘徽魏晋数学家
所谓“割圆术”，是用圆内接正多边形的面积去无限 逼近圆面积以求取圆周率的方法。在刘徽看来， “割 之弥细，所失弥少，割之又割，以至于不可割，则与 圆台体而无所失矣。”
公式推导：
y
b f ( x)dx
a
s1 s2
sn
n
(
i 1
f ( xi1 ) 2
f ( xi )( xi
xi1 )) o
s4 s3 a x1s2 s1 x2 x3
bx
当区间划分为n等分时
——trapz(x,y)
b a
f (x)dx [
f
(a) 2
f (b)
n i 1
ba f (xi )] n , xi
慢
！ 能否算得更快一点，更精确一点
arctan1 arctan1
4
2
3
4arctan1 arctan 1
4
5
239
(Machin公式)
三、数值积分法
基本思想
1 0
1
1 x2
dx
4
4
11 0 1 x2
MATLAB源码求圆周率pi阿基米德法
MATLAB 源码求圆周率pi 阿基米德法%% The Computation of Pi by Archimedes% Copyright 2010, Bill McKeeman, Dartmouth College%% Abstract% It is famously known that Archimedes approximated $\pi$ by computing the % perimeters of many-edged regular polygons, one polygon inside the circle % and one outside. This presentation recapitulates and explains % Archimedes' computation. The surprise to me was how many "tweaks" % Archimedes applied at various stages of an otherwise systematic % approach.%%format longhexpi = imread('diagrams/HexPi.jpg');h=image(hexpi); axis off;%% The General Plan% The circumference of a circle of radius $R$ is $2 \pi R$ which defines % $\pi$. The perimeter of the inner hexagon is $6R$ which implies $6R < 2 % \pi R$, thus providing a lower bound $3 < \pi$. Computing the perimeter % of the outer hexagon supplies an upper bound. Increasing the number of % edges tightens the bounds. Thus thought Archimedes.%% The Context% The year was earlier than 200 BCE. Archimedes, who lived in Sicily, knew % of the work in Alexandria, where Euclid had published his "Elements" some % years before.%%% Greek mathematics used a kind of decimal integer. For example % ${_\prime}\alpha\nu\kappa\delta{^\prime}$ means 1424. Mathematicians % could express rational numbers, but did not have the concept of thedigit % 0, the decimal point or the means to represent irrational numbers such as % $\pi$ and square roots. Irrational numbers were approximated by rational % upper and lower bounds.%%% Finally, algebra was unavailable, so words were used to describe % computations and algorithms. The flavor of ancient presentation is % provided by Archimedes' statement "The surface of any sphere is four % times its greatest circle." We would write $A = 4 \pi R^2$. %% The Computation of $\pi$% The result presented by Archimedes translates to%%% $3 \frac{1}{7} > \pi > 3 \frac{10}{71}$%%% which appears, using the numbers of Archimedes' time, as%%% $\gamma^{\prime}\;\zeta^{\prime\prime} > \pi >% \gamma^{\prime}\;\iota^{\prime}\;o\alpha^{\prime\prime}$%% Advice to the Reader% What follows is a detailed rendition of the steps taken by Archimedes. % While following the logic, it may help the reader to look ahead to the % Summary to get the big picture of all the computations on one page. %% The Outer n-gonouter = imread('diagrams/archimedesOuter.jpg');h=image(outer); axis off;%%% Given the edge length $e_n$ of any regular n-gon enclosing a circle, one % can compute the edge length $e_{2n}$ of the enclosing regular 2n-gon % based on the equations%%% $\begin{array}{rcl}% x_n/R &=& (e_n - e_{2n})/e_{2n} \\% (x_n+R)/R &=& e_n/e_{2n} \\% (x_n+R)/e_n &=& R/e_{2n} \\% x_n^2 &=& R^2+e_n^2% \end{array}$%%% The first equation is a consequence of similar triangles. The second and % third equations are algebraic manipulations of the first. The fourth % equation is an application of the Pythagorean Theorem.%%% The last two equations are combined into an iteration.%%% $\begin{array}{rcl}% R/e_{2n} &=& (x_n+R)/e_n \\% &=& (R + \sqrt{R^2+e_n^2})/e_n \\% &=& R/e_n + \sqrt{(R/e_n)^2+1}% \end{array}$%%% Let $a_n \stackrel{\rm def}{=} R/e_n$.% The result is a recursive formula for $a_{2n}$.%%% $a_{2n} = a_n + \sqrt{a_n^2 + 1}$%%% The value of $a_n$ more than doubles at each step, implying that each % edge $e_n$ is less than one half the length of its predecessor at each % step, which is consistent with one's intuitive understanding of the % process.%%% Archimedes carried out the computation up to the 96-gon, starting with % the 30,60,90 degree triangle bounded by sides%%% $R, x_6, e_6$%%% which has known dimensions. In MATLAB, for the unit circle, R = 1;alpha_6 = 2*pi/12; % 30 degrees, half angle of the hexagon x_6 =R/cos(alpha_6); % R*sqrt(4/3), hypotenusee_6 = x_6*sin(alpha_6); % R*sqrt(1/3), half edge of the hexagon fprintf('e_6 = %.6g*R = R*sqrt(1/3)\n', e_6);%%% Since we know the half-edge $e_6$ for $R=1$, trying the formula is % straightforward.a_6 = R/e_6;a_12 = a_6 + sqrt(a_6^2 + 1);a_24 = a_12 + sqrt(a_12^2 + 1);a_48 = a_24 + sqrt(a_24^2 + 1);a_96 = a_48 + sqrt(a_48^2 + 1);e_96 = R/a_96;%%% The values of $a_n$ are more than doubling; the values $n/a_n$ are slowly % decreasing and converging to $\pi$ as expectedfprintf('n ');fprintf('%7d ', 6, 12, 24, 48, 96);fprintf('\na_n ');fprintf('%.6g ', a_6, a_12, a_24, a_48, a_96);fprintf('\nn/a_n ');fprintf('%.6g ', 6/a_6, 12/a_12, 24/a_24, 48/a_48, 96/a_96);fprintf('... %.6g\n', pi);%%% Unfortunately for Archimedes, he did not have double precision floating % point so he, and therefore we, still have some work to do. %% Hand Calculation% It is simpler for hand calculation if the intermediate results are % improper fractions with reasonably large integer parts. This is achieved % by multiplying the iteration formula by an arbitrary integer constant % $c$, the value of which can be chosen later.%%% $\begin{array}{rcl}% c R/e_{2n} &=& c (x_n+R)/e_n \\ % &=& (c R + c\sqrt{R^2+e_n^2})/e_n \\ % &=& c R/e_n + \sqrt{(c R/e_n)^2+c^2} %\end{array}$%%% Let $a_n \stackrel{\rm def}{=} c R/e_n$ and% $b_n \stackrel{\rm def}{=} \sqrt{(c R/e_n)^2+c^2}$.% The result is new formulas for $b_n$ and $a_{2n}$.%%% $\begin{array}{rcl}% b_n &=& \sqrt{a_n^2 + c^2} \\% a_{2n} &=& a_n + b_n% \end{array}$%%% Two mutually recursive formulas are provided because (1) $b_n$ is % irrational and must be replaced by a rational bound at each step and(2) % the edge of the inner n-gon is expressed in terms of $b_n$, not $a_n$. %% MATLAB fractions% This presentation now switches to "Greek mode," that is, variable % precision integers (vpi) and fractions (fr), provided by John D'Errico % and Ben Petschel on the MathWorks File Exchange. Except for type % conversions putting values into the vpi and fr domains, their packages % are non-intrusive in the code that follows. Thank you to MathWorks, John % and Ben.%% Bounds for Irrational Numbers% Archimedes states without explanation%%% $\frac{265}{153} < \sqrt 3 < \frac{1351}{780}$%%% This statement and other similar ones are easily checked without taking % square roots in MATLAB (and by Archimedes' readers).assert(fr(265/153)^2 < 3 && 3 < fr(1351/780)^2);%% Square Roots of Rationals% Heath, in his book "The Works of Archimedes," devotes 19 pages to % speculations by famous mathematicians on how Archimedes took square % roots. At all times Archimedes had to make a tradeoff between the better % accuracy of larger denominators against the computational labor involved % dealing with larger integers. The speculation I prefer isthat % Archimedes used the known iteration for square root with a good starting % value to get a very accurate fraction, then picked among thecontinued % fraction convergents to get suitable values of lesser accuracy. %%% The following iteration for square root of $a$,%%% $x \leftarrow \frac{1}{2}(x + a/x)$%%% can be restated for rational radicand $a/b$ as%%% $\begin{array}{rcl}% x/y &\leftarrow& \frac{1}{2}(x/y+a y/b x) \\% & \leftarrow & (b x^2 + a y^2)/2 b x y% \end{array}$%%% It is not hard to pick a good-enough starting point for an improper % fraction. If the integer has a odd number of digits and leading digit % $D$, pick $S = 1,2\; \mbox{or}\; 3$ so that $S^2 <= D$. If the integer % part has an even number of digits and leading digits $DE$, pick $S = % 1,\ldots 9$ so that $S^2 <= DE$. Then form the initial value $x$ by % appending to digit $S$ a $0$ for every pair of unexamined digits. That % is, for $349450$, $DE=34, S=5, x=500, y=1$.%%% Carrying out four steps of the computation for input $3$ gives rational % approximations for the $\sqrt 3$.n = 4; % the number of steps to takea = vpi(3);b = vpi(1); % the radicand a/bx = 1; y = 1; % starting values for the rational resultfor i=1:n % all values are now vpit = b*x^2 + a*y^2;y = 2*b*x*y;g = gcd(t,y); % known as Euclid's algorithmx = t/g; % exact divisiony = y/g;fprintf('%.16f %6d/%d\n', double([fr(x)/fr(y), x, y]));endfprintf('...\n%.16f sqrt(%d/%d)\n', ...sqrt(double(a/b)), double(a), double(b));%% Getting the Continued Fraction% Taking the most accurate of the four approximations above, num = vpi(18817); den = 10864;cf = vpi([]); % cf coefficients are always integerswhile den > 0 % compute continued fractionif den == 1 % choose long form of cfif num == 1r = 1; % finish with digit 1den = 0;else % force one more digitr = num - 1;num = 1;endelse % grab a digitr = num/den; % integer divisiont = num - r*den; % new num/dennum = den;den = t;endcf(end+1) = r;end%%% gives the continued fraction for $\frac{18817}{10864}$.fprintf('%d ', double(cf));%%% Evaluating the continued fraction gives a series of convergents, % alternating greater or less than the input fraction. As can be seen in % the printout, they are good approximations to $\sqrt 3$. %% fprintf('x/y (x/y)^2+-err\n'); % table titlen = numel(cf);for i = 1:n % evaluate continued fractionif i == 1 % startup caseh(i) = cf(i);k(i) = 1;elseif i == 2 % startup caseh(i) = cf(i)*h(i-1) + 1;k(i) = cf(i)*k(i-1);elseif i > 2 % rest of cfh(i) = cf(i)*h(i-1) + h(i-2);k(i) = cf(i)*k(i-1) + k(i-2);end% express result as L+-r/K^2 (need double for fprintf)H=double(h(i)); K=double(k(i)); L=round(H^2/K^2); r=H^2-L*K^2;conv = [sprintf('%d/%d', H, K), ' '];if r == 0fprintf('%s %d\n', conv(1:16), L);elseif r < 0fprintf('%s %d-%d/%d^2\n', conv(1:16), L, abs(r), K);elsefprintf('%s %d+%d/%d^2\n', conv(1:16), L, r, K);endend%%% Perhaps it was from this list that Archimedes chose %%% $\frac{1351}{780} > \sqrt 3 > \frac{265}{153}$%%% Since there is not a unique choice, we assume that Archimedes picked % values that would lead to the elegant result for $\pi$. Mathematicians % contemporary to Archimedes mention that later Archimedes did get a more % accurate result.%%% In any case the objective here is $\pi$, not square roots, so from this % point on I shall just report the choices Archimedes made for square % roots.%% The Upper Bound%%% For the upper bound computation, Archimedes chose the initial values c = 153;a_6 = c*R/e_6; % 265.0038b_6 = c*2; % which is precisea_6 = floor(a_6); % choose rational lower boundfprintf('a_6 = %d, b_6 = %d, c = %d\n', a_6, b_6, c);%%% There are two values for $a_6$. The first result comes from the geometry % of the triangle; the second is a chosen rational approximation. In this % case the chosen $a_6$ is just a little low. Further choices for $a_n$ % will be consistently low. The consequence is that rational estimate for % $a_{96}$ is less than the fully accurate irrational value of $a_{96}$. % The value of $e_{96}$ is therefore greater than the fully accurate % half-edge length of the enclosing 96-gon. Finally $96\times e_{96}$ is % greater than the half-perimeter of the enclosing 96-gon which is itself % greater than $\pi$.%% Iterating% Following the reasoning above, each time $b_n$ is computed, it must be % replaced by a rational bound that is smaller than the computed irrational % value. The amount of the difference will be small if the values $b_n$ are % well chosen. Considering Archimedes chose them, we should expect the % best. Nevertheless, we will check. There are three values given by % Archimedes for $b_{12} \ldots b_{48}$.assert(b_6^2 <= (c*R/e_6)^2 + c^2);a_12 = a_6 + b_6; % 571%%% Archimedes chooses $\;b_{12} = 591\frac{1}{8}% < \sqrt{571^2+153^2}% = \sqrt{349450}% = 591.143\ldots$b_12 = fr(591+1/8); assert(b_12^2 < a_12^2+c^2);a_24 = a_12 + b_12; % 1162 1/8%%% Archimedes chooses $b_{24}=1172\frac{1}{8}% < \sqrt{1162\frac{1}{8}^2+153^2} % = \sqrt{1373943\frac{33}{64}} % = 1172.153\ldots$b_24 = fr(1172+1/8); assert(b_24^2 < a_24^2+c^2);a_48 = a_24 + b_24; % 2334 1/4%%% Archimedes chooses $b_{48}=2339\frac{1}{4}% < \sqrt{2334\frac{1}{4}^2+153^2} % = \sqrt{5472132\frac{1}{16}} %= 2339.259\ldots$b_48 = fr(2339+1/4); assert(b_48^2 < a_48^2+c^2);a_96 = a_48 + b_48; % 4673 1/2%%% Archimedes makes one more approximation, replacing the final (ugly) % computed upper bound by one that is less tight but more elegant. pi_hi = fr(3+1/7) % Chosen by Archimedes for publication pi_est1 = 96*c/a_96 % Computed by Archimedespi % MATLAB's best%%% Using MATLAB, we can check that the values are ordered as expected. assert(pi_hi > pi_est1 && pi_est1 > pi);%% The Lower Bound% There is another diagram (in Heath), and another set of geometric% equations for the lower bounds which leads to the same recursive % formulas. One should not be surprised. In both cases the n-gons are the % same; it is the scale that differs.%%% This time the chosen rational values of $a_n$ must be greater than the % true values, causing the computed half-perimeter $96\times e_{96}$ to be % less than the actual half-perimeter of the inner 96-gon. %% % The initial values area_6 = 1351; b_6 = 1560; c = 780; assert(a_6^2+c^2 > b_6^2); %% Iterating% The values chosen for $b_n$ therefore must also be greater thanthe % irrational true values. Repeating the iteration, and in this case also % removing some common factors:a_12 = a_6 + b_6; % 2911%%% Archimedes chooses $b_{12}=3013\frac{3}{4}% > \sqrt{2911^2+780^2}% = \sqrt{9082321}% = 3013.689\ldots$b_12 = fr(3013+3/4); assert(b_12^2 > a_12^2+c^2);a_24 = a_12 + b_12; % 5924 3/4%%% Archimedes removes the common factor $\frac{13}{4}$ from $a_{24}$ and % $c$.a_24 = a_24*(4/13); % 1823c = c*(4/13); % 240%%% Archimedes chooses $b_{24}=1838\frac{9}{11} % > \sqrt{1823^2+240^2} % = \sqrt{3380929} % = 1838.730\ldots$ b_24 = fr(1838+9/11);assert(b_24^2 > a_24^2+c^2); a_48 = a_24 + b_24; % 3661 9/11 %%% Archimedes removes the common factor $\frac{40}{11}$ a_48 =a_48*(11/40); % 1007c = c*(11/40); % 66%%% Archimedes chooses $b_{48}=1009\frac{1}{6} % > \sqrt{1007^2+66^2} % = \sqrt{1018405} % = 1009.161\ldots$ b_48 = fr(1009+1/6);assert(b_48^2 > a_48^2+c^2); a_96 = a_48 + b_48; % 2016 1/6 %%% Archimedes chooses $b_{96}=2017\frac{1}{4} % >\sqrt{2016\frac{1}{6}^2+66^2}% = \sqrt{4069284\frac{1}{36}} % = 2017.247\ldots$ b_96 =fr(2017+1/4); assert(b_96^2 > a_96^2+c^2);pi_lo = fr(3+10/71); % Chosen by Archimedes for publication pi_est2 = 96*c/b_96; % Computed by Archimedes %%% The five values aredisp(pi_hi);disp(pi_est1);disp(pi);disp(pi_est2);disp(pi_lo);%%% The difference between the published bounds is disp(pi_hi - pi_lo) %%% Check that the five values are in decreasing order assert(pi_lo < pi_est2 && pi_est2 < pi && pi < pi_est1 && pi_est1 < pi_hi);%% QED% Therefore, $3\frac{1}{7} > \pi > 3\frac{10}{71}$.%%% Archimedes might have written: $OE\Delta$%% Summary% Here are all of the intermediate results in a table. Most of the action % is in the choosing of $b_n$. Recall that the chosen values of $b_n$ for % the outer case must be less than the computed values and conversely % greater for the inner case.%%% $\begin{array}{lllcllrl}% & \rm outer &&&& \rm inner && \rm remark \\ % a & b & c &\mathbf{n} & a & b & c \\ % 265 & 306 & 153 & \mathbf{6} & 1351 & 1560 & 780 & \rm start \\ % 571 & \sqrt{571^2+153^2} & 153 & \mathbf{12} % & 2911 & \sqrt{2911^2+780^2}\makebox[2pt]{} & 780 \\ % &591.1430\ldots &&&& 3013.6889\ldots && \rm compute \\ % & 591\frac{1}{8} &&&& 3013\frac{3}{4} && \rm choose \\ % &&& \mathbf{24}% & 5924\frac{3}{4}\makebox[5pt]{} && 780% & \mbox{factor out}\;\frac{13}{4} \\ % 1162\frac{1}{8} &\sqrt{1162\frac{1}{8}^2+153^2} & 153 & % & 1823 & \sqrt{1823^2+240^2} & 240 \\ % & 1172.1534\ldots &&&& 1838.7303\ldots && \rm compute \\ % & 1172\frac{1}{8} &&&& 1838\frac{9}{11} && \rm choose % \end{array}$%% $\begin{array}{lllcllrl}% && & \mathbf{48}% & 3661\frac{9}{11} && 240 & \mbox{factor out}\;\frac{40}{11} \\ % 2334\frac{1}{4} & \sqrt{2334\frac{1}{4}^2+153^2} & 153 % && 1007 &\sqrt{1007^2+66^2} & 66 \\ % & 2339.2589\ldots &&&& 1009.1605\ldots && \rm compute \\ % & 2339\frac{1}{4} &&&& 1009\frac{1}{6} && \rm choose \\ % 4673\frac{1}{2} & (b\;\mbox{not needed}) & 153 & \mathbf{96} % &2016\frac{1}{6} & \sqrt{2016\frac{1}{6}^2+66^2} & 66 \\ % &&&&&2017.2467\ldots && \rm compute \\ % &&&&& 2017\frac{1}{4} && \rm choose % \end{array}$%% $\begin{array}{cc}% \rm outer\;result & \rm inner\;result \\ % 96\times 153/4673\frac{1}{2}=3\frac{1335}{9347}=3.1428265982\ldots & % 96\times 66 / 2017\frac{1}{4}=3 \frac{1137}{8069}=3.1409096542\ldots% \end{array}$%% References% My principal references are "The Works of Archimedes" (1897), "The Method % of Archimedes Recently Discovered by Heiberg" (1912) and "A History of % Greek Mathematics" (1921) all by Sir Thomas Heath. The 1921 book is % still in print; Dartmouth Library dug wonderful dusty old tomes from % storage for the rest. I also watched an enjoyable lectureon Archimedes % in "Great Thinkers, Great Theorems" by Professor William Dunham (Great % Courses series) and read some entries in Wikipedia.。
(2021年整理)基于MATLAB的π
基于MATLAB的π编辑整理：尊敬的读者朋友们：这里是精品文档编辑中心，本文档内容是由我和我的同事精心编辑整理后发布的，发布之前我们对文中内容进行仔细校对，但是难免会有疏漏的地方，但是任然希望（基于MATLAB的π）的内容能够给您的工作和学习带来便利。
同时也真诚的希望收到您的建议和反馈，这将是我们进步的源泉，前进的动力。
本文可编辑可修改，如果觉得对您有帮助请收藏以便随时查阅，最后祝您生活愉快业绩进步，以下为基于MATLAB的π的全部内容。
基于MATLAB 的4πQPSK 的仿真实现（刘海亮 通信1201 1230440103）摘 要在以前的数字蜂窝系统中，往往采用FSK 、ASK 、PSK 等调制方式。
随着数字蜂窝系统的发展，对调制和数字蜂窝系统的技术要求越来越高， 许多优秀的调制技术应运而生,其中π/4QPSK 技术是无线通信中比较突出的一种二进制调制方法本文概述了π/4QPSK 的调制解调原理及其所实现的功能。
并通过MATLAB 编程对系统在相同条件下，对比了加噪声和没有噪声的情况进行了比较，并画出了它的眼图，及已调制信号的时域波形和频谱。
并通过眼图分析其性能。
关键字：π/4QPSK;调制解调原理；MATLAB 编程；眼图1引言无线通信在现代社会中起着举足轻重的作用。
从日常生活到航空航天,从工商业运作到军事领域，无线通信得到了越来越广泛的应用。
现代数字调制技术的发展，使得传输速率和频谱的利用率进一步得到提高，功率更加节省。
在相同的码元速率下，多进制系统的信息传输速率显然比二进制系统高，但信息速率的提高是以牺牲功率为代价的。
显然增大码元宽度，就会增加码元的能量，同时也减少了由于信道特性引起的码间串扰等。
恒包络调制适用于限带非线性信道中，能有效地防止非线性引起的幅频效应，节省功率,提高频谱的利用率.多进制调制和恒包络调制这两种技术结合在一起能取得更好的调制效果。
为了使基带信号更好的利用信道进行传输，必须使代表的原始信号经过调制，而调制技术的好坏影响频谱资源的利用和通信性能的好坏。
数学建模——π的计算
数学建模——π的计算实验报告实验目的1. 掌握用MATLAB软件求π值的方法，并对结果作初步分析；2.掌握用MATLAB软件实现蒙特卡洛法的过程。
实验内容运用数值积分法、泰勒级数法、模拟蒲丰实验和随机整数互素法计算π值，取n =（k=3, 4, …,7），记录所得结果和程序运行时间，并作简要分析。
实验过程一、数值积分法单位圆的面积等于以单位圆的圆心为原点建立直角坐标系，单位圆在第一象限内的部分G 是一个扇形，其面积S = π/4只需计算出S 的近似值，它的4倍就是π的近似值1、梯形公式法将扇形G 分为n 个宽度相等的部分Gk（1 ≤k ≤n）Gk：曲边梯形左右边界平行，上方边界为曲线n→∞Gk→梯形Gk 梯形面积：等价于π取 n =10k （k=3, 4, …,7），按上述方法，通过计算扇形面积计算π的近似值。
程序Clearticn=1e 3(3, 4, …,7); x=linspace(0,1,n+1); y=(1-x.^2).^0.5;pai=(0.5*y(1)+0.5*y(n+1)+sum(y(2:n)))*4/n tocG G2、辛普森（Simpson）公式法仍用分点xi = a + i(b-a)/n (1≤ i ≤ n-1) 将区间[a, b]分成n 等份，直线x = xi (1≤ i ≤ n-1) 将曲边梯形分成n 个小曲边梯形，再作每个小区间[xi-1, xi]的中点将第i个小曲边梯形的上边界y=f(x) (xi-1≤x≤xi) 近似地看作经过这三点的抛物线段，则可求得：其中于是得到S;即程序：clearticn=1e3(3,4,5,6,7);x=linspace(0,1,n+1);x2=x(1:n)+1/(2*n);y=4./(1+x.^2);y2=4./(1+x2.^2);pai=((y(1)+y(n+1))+2*sum(y(2:n))+4*sum(y2))/(6*n) toc二、泰勒级数法考虑反正切函数的泰勒级数取x =1，n =1e3(3,4,5,6,7)程序：cleartick=1e3(3,4,5,6,7);x=1;n=1:k;pai=sum((-1).^(n-1).*x.^(2.*n-1)./(2.*n-1))toc结果显示：花费时间长！准确度差！！原因：x =1时得到的arctan1的展开式收敛太慢！Maqin公式法提高收敛速度：x << 1随着指数的增加，x的幂快速接近于0，泰勒级数就会快速收敛取2k-1 = 63，其误差已经非常小，表明收敛速度非常快有再考虑收敛更快的：/4=arctan 1/2+arctan 1/3程序：cleartick=1e3(3,4,5,6,7);n=1:k;x=1/2;at_1=sum((-1).^(n-1).*x.^(2.*n-1)./(2.*n-1));x=1/3;at_2=sum((-1).^(n-1).*x.^(2.*n-1)./(2.*n-1));pai_1=4*at_1+4*at_2toc/4=4arctan 1/5-arctan 1/239程序：cleartick=1e3(3,4,5,6,7);n=1:k;x=1/5;at_5=sum((-1).^(n-1).*x.^(2.*n-1)./(2.*n-1));x=1/239;at_239=sum((-1).^(n-1).*x.^(2.*n-1)./(2.*n-1));pai=16*at_5-4*at_239toc三、蒙特卡罗（Monte Carlo）法在数值积分中，我们利用求1/4单位圆的面积来得到π/4，进而得到π 1/4单位圆是一个扇形G，它是边长为1的单位正方形G1的一部分。
蒙特卡洛方法求圆周率的matlab程序
蒙特卡洛方法是一种利用随机抽样来估计数学问题的数值解的方法。
在数值积分和求解难以解析的概率统计问题时，蒙特卡洛方法经常能够取得比较好的结果。
在本文中，我将详细介绍如何使用蒙特卡洛方法来求解圆周率，并给出相应的MATLAB程序。
1. 蒙特卡洛方法求解圆周率的原理蒙特卡洛方法求解圆周率的原理是基于统计学中的随机抽样原理。
我们知道，圆的面积公式为S=πr^2，而圆的半径r=1。
通过在一个边长为2的正方形区域内随机散布大量的点，我们可以通过统计正方形内部与圆内部的点的比例来估计圆的面积，从而得到圆周率的近似值。
2. MATLAB程序编写步骤我们需要生成大量的随机点，这些点需要均匀分布在正方形区域内。
我们统计这些点中有多少落在了圆的内部。
通过统计得到的比例，我们可以计算出圆的面积，从而得到圆周率的估计值。
下面给出蒙特卡洛方法求解圆周率的MATLAB程序：``` MATLABfunction pi_estimate = monte_carlo_pi(n)% n为随机点的数量count_inside_circle = 0;for i=1:nx = 2*rand()-1; % 生成-1到1之间的随机数x坐标y = 2*rand()-1; % 生成-1到1之间的随机数y坐标if x^2 + y^2 <= 1 % 判断点是否落在圆的内部count_inside_circle = count_inside_circle + 1;endendpi_estimate = 4 * count_inside_circle / n; % 计算圆周率的估计值end```3. 程序使用说明通过调用上述的MATLAB函数monte_carlo_pi，传入随机点的数量n，即可得到圆周率的估计值。
n越大，估计值越接近真实值。
一般来说，n的取值在几万到几百万之间时，可以得到比较准确的结果。
下面给出一个调用例子：``` MATLABn = 1000000; % 随机点数量为100万pi_estimate = monte_carlo_pi(n); % 调用函数求解圆周率disp(['使用', num2str(n), '个随机点，得到的圆周率的估计值为：', num2str(pi_estimate)]);```4. 结论蒙特卡洛方法是一种有效的数值计算方法，在求解圆周率等复杂数学问题时具有一定的优势。
matlab中表示圆周率的固定变量
matlab中表示圆周率的固定变量Matlab是一种高级数学软件，被广泛应用于科学、工程和金融等领域。
在Matlab中，圆周率是一个重要的常数，它被用于许多数学和工程计算中。
本文将介绍Matlab中表示圆周率的固定变量。
一、圆周率的定义圆周率是一个数学常数，通常用希腊字母π表示。
它的定义是：圆的周长与直径的比值。
即π=周长/直径。
因为圆的周长和直径都是一定的，所以π也是一个固定的数值。
在数学中，π是一个无理数，其小数点后的数字是无限的，且没有规律可循。
二、Matlab中表示圆周率的方法在Matlab中，表示圆周率的常量是pi。
这个常量是预定义的，因此不需要定义或初始化。
只需要在Matlab中输入pi，即可获得π的值。
例如：pians =3.1416Matlab中的pi是一个浮点数，它的值同样是无理数。
Matlab中的pi可以用于各种数学和工程计算中，例如圆的面积和周长的计算。
三、Matlab中圆面积和周长的计算在Matlab中，可以使用pi常量来计算圆的面积和周长。
以下是计算圆面积和周长的代码：% 计算圆的半径r = 5;% 计算圆的面积area = pi * r^2;% 计算圆的周长circumference = 2 * pi * r;% 输出结果disp(['圆的半径为：', num2str(r)]);disp(['圆的面积为：', num2str(area)]);disp(['圆的周长为：', num2str(circumference)]);运行以上代码，输出结果如下：圆的半径为：5圆的面积为：78.5398圆的周长为：31.4159以上代码中，首先定义了圆的半径r为5，然后使用pi常量计算了圆的面积和周长。
最后使用disp函数输出结果。
四、Matlab中的圆的绘制在Matlab中，可以使用plot函数绘制圆。
以下是绘制圆的代码： % 计算圆的半径r = 5;% 计算圆的坐标theta = linspace(0, 2*pi, 100);x = r * cos(theta);y = r * sin(theta);% 绘制圆plot(x, y);axis equal;运行以上代码，将绘制一个半径为5的圆。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
3
5
7
2 n 1
并利用这个公式计算到了圆周率的100位.
1914年，印度数学家Srinivasa Ramanujan 表了下面的公式：
9801 ( 4 n )! ( 1103 26390 n ) 22 4 n 4 4 n 4 ( n ! ) 99 n 0
发
在1985年，Gosper用这个公式计算到了圆周率 的17500000位．
在中国
• 祖冲之: 在刘徽研究的基础上，进一步地发展， 经过既漫长又烦琐的计算，一直算到圆内接正 24576边形，而得到一个结论： • 3.1415926 ＜ π ＜ 3.1415927 同时得到π 的两个近似分数：约率为22／7； 密率为355／113。
• 他算出的 π 的8位可靠数字，不但在当时是最精 密的圆周率，而且保持世界记录九百多年。以致 于有数学史家提议将这一结果命名为“祖率”。
1 .圆周率 的幂级数计算方法
例1
利用 arctan x 的 Maclaurin 展开式 ， 计算 的近似值 .
解
1 2 4 n 12 n 1 x x (1 ) x . 2 1 x
xx ( 1 )x arctan x x . 35 2 n 1
1989年，David 和 Gregory Chudnovsky 发表 了下面的公式
n 1 ( 1 ) ( 6 n )! 13591409 5451401 n 12 , 3 3 3 n n ( 3 n )! ( n ! ) 2 0 640320
并在1994年计算到了4044000000位．它的另一 种形式是
在中国
• 刘徽：公元263年前后，刘徽提出著名的 “割圆 术”求出了比较精确的圆周率。他发现：当圆内 接正多边形的边数不断增加后，多边形的周长会 越来越逼近圆周长，而多边形的面积也会越来越 逼近圆面积。于是，刘徽利用正多边形面积和圆 面积之间的关系，从正六边形开始，逐步把边数 加倍：正十二边形、正二十四边形，正四十八边 形……，一直到正三○七二边形，算出圆周率等 于三点一四一六，将圆周率的精度提高到小数点 后第四位。
426880 10005 . ( 6 n )! ( 545140134 n 13591409 ) 3 3 n ( n ! ) ( 3 n )! ( 640320 ) n 0
2019 年 ， 由 David Bailey,Peter Borwein 和 Simon Plouffe 共同发表了下面的圆周率计算公式 （简称BBP公式）
数学实验五
圆周率 π 的近似 计 算
主 讲:魏 平
１.圆周率π的计算历程
• 所谓“圆周率”是指一个圆的周长与其直径的比值。 古今中外，许多人致力于圆周率的研究与计算。为 了计算出圆周率的越来越好的近似值，一代代的数 学家为这个神秘的数贡献了无数的时间与心血。。 • 回顾历史，人类对 π 的认识过程，反映了数学和计 算技术发展情形的一个侧面。 π 的研究，在一定程 度上反映这个地区或时代的数学水平。德国数学家 康托说：“历史上一个国家所算得的圆周率的准确 程度，可以作为衡量这个国家当时数学发展水平的 指标。” • 直到19世纪初，求圆周率的值应该说是数学中的头 号难题。为求得圆周率的值，人类走过了漫长而曲 折的道路。
1 1 32 1 256 64 1024 n1 4 n3 10 n1 4 n 0
n
64 4 4 1 . 10 n 3 10 n 5 10 n 7 10 n 9
从而，大大降低了圆周率近似值的计算量.
1 2 1 1 4 . n 8 n 1 8 n 4 8 n 5 8 n 6 16 n 0
该公式的最大优点在于：经后来人将该公式变 形后打破了传统的计算方法，可以直接计算圆周率 的任意第n位数，而不是先计算前面的n-1位数．
2019 年， Fabrice Bellard 发表了一个比 BBP 算 法更快的公式
Hale Waihona Puke 分析法时期• 这一时期人们开始摆脱求多边形周长的繁难 计算，利用无穷级数或无穷连乘积来算 π 。 • 1593年，韦达给出
2 2 2 2 2 2 2 2 2 2
这一不寻常的公式是 π 的最早分析表达式。甚至 在今天，这个公式的优美也会令我们赞叹不已。它 表明仅仅借助数字2，通过一系列的加、乘、除和 开平方就可算出 π 值。
接着有多种表达式出现。如沃利斯1650 年给出：
22446688
2 1 3 345 577
1706年，英国天文学教授John Machin 利用
xxx x n 1 arctan x x 1 357 2 n 1
发现了下面的公式
1 1 16 arctan 4 arctan , 5 239
实验时期
• 基于对一个圆的周长和直径的实际测量而 得出的。 • 在古代世界，实际上长期使用 π ＝3这个数 值。 • 最早见于文字记载的有基督教《圣经》中 的章节，其上取圆周率为3。这一段描述的 事大约发生在公元前950年前后。
几何法时期
真正使圆周率计算建立在 科学的基础上，首先应归 功于阿基米德。他是科学 地研究这一常数的第一个 人，是他首先提出了一种 能够借助数学过程而不是 通过测量的、能够把 π 的 圆周长大于内接正多边 值精确到任意精度的方法。 形周长而小于外切正多边 由此，开创了圆周率计算 形周长． 的第二阶段。 据说阿基米德用到了正 96边形才算出他的值域。
3 5
n 1 2 n 1
1 1 ( 1 ) 4 arctan 1 4 ( 1 ). 3 5 2 n 1
n 1
1 1 ( 1 ) 4 arctan 1 4 ( 1 ). 3 5 2 n 1