MATLAB模拟蒲丰掷针实验计算PI的值
matlab模拟蒲丰掷针实验计算pi的值
%投掷次数为n,针和直线相交的次数为m
%分析知真和直线相交的概率为p=2l/pid
%则有pi=2nl/md
clear all
clc
d=1;%设置两条平行线之间的距离,相当于在一张纸上画很多间距为d的平行线
l=0.6;%投针的长度（l n=1000000;%n为投掷次数 x=unifrnd(0,d/2,1,n);%产生n个(0,d/2)之间均匀分布的随机数，这里d/2是投针的中点到最近的平行线的距离 alpha=unifrnd(0,pi,1,n);%产生n个(0,pi)之间均匀分布的随机数，这里pi是投针到最近的平行线的角度 m=0;%记针与平行线相交的次数初始值为0 for i=1:n if x(i) m=m+1; end end p=m/n;%计算相交的频率，即相交次数比总次数 Vpi=2*l/(d*p)%从相交的频率总求的Value pi,即pi的值 %运行结果： %Vpi= % 3.140710093614099 %因每次产生的随机数都不一样，每次运算出来的pi的值都不同。n越大，越精确。 系统建模与仿真题目：Buffon实验的仿真 院系: 电子工程学院 专业：信息对抗技术 班级：021231 姓名：余颖智 学号：02123021 指导老师：刘洋 完成时间：2015年4月 西安电子科技大学 基于MATLAB的投针实验仿真 摘要 在求证圆周率的过程中经过割圆术后，出现的投针试验以求出圆周率，目前利用MATLAB数学建模的仿真实验，运用到计算机中，简化其随机实验的操作量大，运算慢等特点。不同针距相同实验量运算后得出不同的π，其针距与线间距离相等，所得值接近于π。 目录 摘要 (2) 二、实验内容 (4) 三、建模流程图 (5) 四、程序主要代码 (6) 五、运行结果 (6) 六、结论 (7) 一、实验原理 1777年法国科学家布丰提出的一种计算圆周率的方法——随机投针法，即著名的布丰投针问题。该投针实验主要有如下三个步骤：(一)取一张白纸，在上面画许多条间距为a的平行线；（二）取一根长度为l（l 三、建模流程图 四、程序主要代码 str(handles.edit1,'string'); %取得变量，定义变量，变量初始化 n = str2double(str); str = get(handles.edit2,'string'); l = str2double(str); str = get(handles.edit3,'string'); a = str2double(str); counter = 0; %变量初始化 phi = 0; frequency = 0; Pi = 0; x = unifrnd(0,a/2,1,n);%产生n个（0，a/2）之间均匀分布的随机数，这里a/2是投针的中点到最近的平行线的距离 phi = unifrnd(0,pi,1,n);% 产生n个（0，pi）之间均匀分布的随机数，这里pi是投针与最近平行线的角度 for i=1:n if x(i) 6-1 试验目的计算特征值，实现算法 试验容：随机产生一个10阶整数矩阵，各数均在-5和5之间。 （1） 用MATLAB 函数“eig ”求矩阵全部特征值。 （2） 用幂法求A 的主特征值及对应的特征向量。 （3） 用基本QR 算法求全部特征值（可用MATLAB 函数“qr ”实现矩阵的QR 分解）。 原理 幂法：设矩阵A 的特征值为12n ||>||||λλλ≥???≥并设A 有完全的特征向量系12,,,n χχχ???(它们线性无关)，则对任意一个非零向量0n V R ∈所构造的向量序列1k k V AV -=有11()lim ()k j k k j V V λ→∞ -=， 其中()k j V 表示向量的第j 个分量。 为避免逐次迭代向量k V 不为零的分量变得很大（1||1λ>时）或很小（1||1λ<时），将每一步的k V 按其模最大的元素进行归一化。具体过程如下： 选择初始向量0V ，令1max(),,,1k k k k k k k V m V U V AU k m +===≥，当k 充分大时1111,max()max() k k U V χλχ+≈ ≈。 QR 法求全部特征值： 111 11222 111 ,1,2,3,k k k k k A A Q R R Q A Q R k R Q A Q R +++==????==??=???? ??????==?? 由于此题的矩阵是10阶的，上述算法计算时间过长，考虑采用改进算法——移位加速。迭 代格式如下： 1 k k k k k k k k A q I Q R A R Q q I +-=?? =+? 计算k A 右下角的二阶矩阵() () 1,1 1,() (),1 ,k k n n n n k k n n n n a a a a ----?? ? ??? 的特征值()()1,k k n n λλ-，当()()1,k k n n λλ-为实数时，选k q 为()()1,k k n n λλ-中最接近(),k n n a 的。 程序 %求四阶线性方程组的MA TLAB程序 clear Ab=[0.001 2 1 5 1; 3 - 4 0.1 -2 2; 2 -1 2 0.01 3; 1.1 6 2.3 9 4];%增广矩阵 num=[1 2 3 4];%未知量x的对应序号 for i=1:3 A=abs(Ab(i:4,i:4));%系数矩阵取绝对值 [r,c]=find(A==max(A(:))); r=r+i-1;%最大值对应行号 c=c+i-1;%最大值对应列号 q=Ab(r,:),Ab(r,:)=Ab(i,:),Ab(i,:)=q;%行变换 w=Ab(:,c),Ab(:,c)=Ab(:,i),Ab(:,i)=w;%列变换 n=num(i),num(i)=num(c),num(c)=n;%列变换引起未知量x次序变化for j=i:3 Ab(j+1,:)=-Ab(j+1,i)*Ab(i,:)/Ab(i,i)+Ab(j+1,:);%消去过程 end end %最后得到系数矩阵为上三角矩阵 %回代算法求解上三角形方程组 x(4)=Ab(4,5)/Ab(4,4); x(3)=(Ab(3,5)-Ab(3,4)*x(4))/Ab(3,3); x(2)=(Ab(2,5)-Ab(2,3)*x(3)-Ab(2,4)*x(4))/Ab(2,2); x(1)=(Ab(1,5)-Ab(1,2)*x(2)-Ab(1,3)*x(3)-Ab(1,4)*x(4))/Ab(1,1); for s=1:4 fprintf('未知量x%g =%g\n',num(s),x(s)) end %验证如下 %A=[0.001 2 1 5 1; 3 -4 0.1 -2 2;2 -1 2 0.01 3; 1.1 6 2.3 9 4]; %b=[1 2 3 4]'; %x=A\b; %x1= 1.0308 %x2= 0.3144 %x3= 0.6267 %x4= -0.0513 几何概型的解析方法和类型例解 及蒲丰投针试验、贝特朗奇论的解析 几何 概型是高中数学新课程的新增加内容之一.部分师生在几何概型知识的理解上存在一些偏差，在有关内容的教学中有说理不清或解法错误的现象。本文依据几何概型理论知识和课例，谈谈几何概型问题的解析方法和基本事件为线段、圆、球射线等非质点问题的转化，以及蒲丰投针试验、贝特朗奇论的解析。 1 问题引入 题 1 如图所示，A 、B 两盏路灯之间长度是30米，由于光线较暗，想在其间再随意安装两盏路灯C 、D ，问A 与C ，B 与D 之间的距离都不小于10米的概率是多少？ 在第99页里和网上一些教师课件中有上述习题及如下答案。 解：：记E ：“A 与C ，B 与D 之间的距离都不小于10米”，把AB 三等分，由于中间长度为130103? =米，所以101()303 P E ==。 笔者认为：上述解答错误。错误主要原因几何概型问题的解析方法不正确。 题2 欧阳修《卖油翁》中写到：（翁）乃取一葫芦置于地，以钱覆其口，徐以杓酌油沥之，自 钱孔入，而钱不湿．可见“行行出状元”，卖油翁的技艺让人叹为观止．若铜钱是直径为4cm 的圆，中间有边长为1cm 的正方形孔，若你随机向铜钱上滴一滴油（油滴不出铜钱面），油滴的直径是0.2cm 的球，则油滴整体落入孔中的概率是 。 关于此题的教学中，部分学生和教师出现分析无方，说理不清等现象，那么几何概型应怎样分析解决呢？ 题3 苏教版必修3提到的“贝特朗奇论”：在半径为1的圆内任作一条弦，求该弦的长度长于圆的内接正三角形边长的概率。 贝特朗奇论的经典解法有如下三种。 解法一：如图1弦被其中点位置唯一确定。只有当弦的中点落在半径缩小了一半的同心圆内，其长才合乎要求。中点位置都是等可能的,则所求概率为4 1。 解法二：如图2由于对称性，可预先指定弦的方向。作垂直于此方向的直径，只有交直径于41点与43 点间的弦，其长才大于内接正三角形边长。所有交点是等可能的,则所求概率为2 1 。 图 1 系统建模与仿真 基于MATLAB的布丰实验模拟 姓名：石星宇 学号： 02123010 指导教师：刘洋 2015年4月9日 目录 基于MATLAB的布丰实验模拟 .................................................................... - 1 - 一、实验原理......................................................................................... - 1 - 二、编程模拟......................................................................................... - 1 - 1、程序流程图............................................................................... - 1 - 2、程序代码................................................................................... - 2 - 三、实验结果......................................................................................... - 2 - 基于MATLAB 的布丰实验模拟 一、实验原理 找一根铁丝弯成一个圆圈，使其直径恰恰等于平行线间的距离a 。可以想象得到，对于这样的圆圈来说，不管怎么扔下，都将和平行线有两个交点。因此，如果圆圈扔下的次数为n 次，那么相交的交点总数必为n 2。现在设想把圆圈拉直，变成一条长为a π的铁丝。显然，这样的铁丝扔下时与平行线相交的情形要比圆圈复杂些，可能有4个交点，3个交点，2个交点，1个交点，甚至于都不相交。由于圆圈和直线的长度同为a π，根据机会均等的原理（即等概率事件），当它们投掷次数较多，且相等时，两者与平行线组交点的总数期望也是一样的。这就是说，当长为a π的铁丝扔下n 次时，与平行线相交的交点总数应大致为n 2。现在转而讨论铁丝长为l 的情形。当投掷次数n 增大的时候，这种铁丝跟平行线相交的交点总数k 应当与长度l 成正比，因而有：l k λ=，式中λ是比例系数。为了求出λ来，只需注意到，对于a l π=的特殊情形，有n k 2=。于是求得a n πλ2=。代入前式就有：a m πln 2≈从而ak nl 2≈π。 二、编程模拟 1、程序流程图 参数初始化 产生位置随机数； 产生角度随机数 判断相交 1+=k k 1+=n n 是 否 判断结束 常用计算方法 1．超越方程的求解 一超越方程为 x (2ln x – 3) -100 = 0 求超越方程的解。 [算法]方法一：用迭代算法。将方程改为 01002ln()3 x x =- 其中x 0是一个初始值，由此计算终值x 。取最大误差为e = 10-4，当| x - x 0| > e 时，就用x 的值换成x 0的值，重新进行计算；否则| x - x 0| < e 为止。 [程序]P1_1abs.m 如下。 %超越方程的迭代算法 clear %清除变量 x0=30; %初始值 xx=[]; %空向量 while 1 %无限循环 x=100/(2*log(x0)-3); %迭代运算 xx=[xx,x]; %连接结果 if length(xx)>1000,break ,end %如果项数太多则退出循环（暗示发散） if abs(x0-x)<1e-4,break ,end %当精度足够高时退出循环 x0=x; %替换初值 end %结束循环 figure %创建图形窗口 plot(xx,'.-','LineWidth',2,'MarkerSize',12)%画迭代线'.-'表示每个点用.来表示，再用线连接 grid on %加网格 fs=16; %字体大小 title('超越方程的迭代折线','fontsize',fs)%标题 xlabel('\itn','fontsize',fs) %x 标签 ylabel('\itx','fontsize',fs) %y 标签 text(length(xx),xx(end),num2str(xx(end)),'fontsize',fs)%显示结果 [图示]用下标作为自变量画迭代的折线。如P0_20_1图所示，当最大误差为10-4时，需要迭代19次才能达到精度，超越方程的解为27.539。 [算法]方法二：用求零函数和求解函数。将方程改为函数 100()2ln()3f x x x =-- MATLAB 求零函数为fzero ，fzero 函数的格式之一是 x = fzero(f,x0) 其中，f 表示求解的函数文件，x0是估计值。fzero 函数的格式之二是 x = fzero(f,[x1,x2]) 1、Newdon迭代法求解非线性方程 function [x k t]=NewdonToEquation(f,df,x0,eps) %牛顿迭代法解线性方程 %[x k t]=NewdonToEquation(f,df,x0,eps) %x:近似解 %k:迭代次数 %t:运算时间 %f:原函数，定义为内联函数 ?:函数的倒数，定义为内联函数 %x0:初始值 %eps:误差限 % %应用举例： %f=inline('x^3+4*x^2-10'); ?=inline('3*x^2+8*x'); %x=NewdonToEquation(f,df,1,0.5e-6) %[x k]=NewdonToEquation(f,df,1,0.5e-6) %[x k t]=NewdonToEquation(f,df,1,0.5e-6) %函数的最后一个参数也可以不写。默认情况下，eps=0.5e-6 %[x k t]=NewdonToEquation(f,df,1) if nargin==3 eps="0".5e-6; end tic; k=0; while 1 x="x0-f"(x0)./df(x0); k="k"+1; if abs(x-x0) < eps || k >30 break; end x0=x; end t=toc; if k >= 30 disp('迭代次数太多。'); x="0"; t="0"; end 2、Newdon迭代法求解非线性方程组 function y="NewdonF"(x) %牛顿迭代法解非线性方程组的测试函数 %定义是必须定义为列向量 y(1,1)=x(1).^2-10*x(1)+x(2).^2+8; y(2,1)=x(1).*x(2).^2+x(1)-10*x(2)+8; return; function y="NewdonDF"(x) %牛顿迭代法解非线性方程组的测试函数的导数 y(1,1)=2*x(1)-10; y(1,2)=2*x(2); y(2,1)=x(2).^+1; y(2,2)=2*x(1).*x(2)-10; return; 以上两个函数仅供下面程序的测试 function [x k t]=NewdonToEquations(f,df,x0,eps) %牛顿迭代法解非线性方程组 %[x k t]=NewdonToEquations(f,df,x0,eps) %x:近似解 %k:迭代次数 %t:运算时间 %f:方程组（事先定义） ?:方程组的导数（事先定义） %x0:初始值 %eps:误差限 % %说明：由于虚参f和df的类型都是函数，使用前需要事先在当前目录下采用函数M文件定义% 另外在使用此函数求解非线性方程组时，需要在函数名前加符号“@”，如下所示 % %应用举例： %x0=[0,0];eps=0.5e-6; %x=NewdonToEquations(@NewdonF,@NewdonDF,x0,eps) %[x k]=NewdonToEquations(@NewdonF,@NewdonDF,x0,eps) %[x k t]=NewdonToEquations(@NewdonF,@NewdonDF,x0,eps) %函数的最后一个参数也可以不写。默认情况下，eps=0.5e-6 %[x k t]=NewdonToEquations(@NewdonF,@NewdonDF,x0,eps) 蒲丰投针问题 1.蒲丰简介 蒲丰有的时候翻译成布丰，是18世纪法国著名 的博物学家。他喜欢研究数学和生物学。主要的贡献 有：（1）翻译了牛顿的《流数法》，流数法按现在的 说法就叫微积分。（2）写了一本巨著，这部巨著的名 字叫《自然史》，因为他特别喜欢研究生物。这个自 然史一共有44卷，其中他生前写了36卷，后来他学 生又完成了。这本书对后来的世界有很大的影响，尤 其影响到一个人叫达尔文，所以蒲丰这个人其实是很 厉害的。 2.蒲丰投针 1777年，在蒲丰晚年的时候，他有一次举行了一 个家庭宴会。邀请了一大堆他的朋友来帮他做实验。 做什么实验呢，就“投针”。那朋友来了之后发现，就 是桌子上有很多根间距相等的平行线。然后蒲丰就说 了，给你们同样大的针，你把这些针随机扔到这个桌子上。然后宾客就随便扔吗，有可能这样，有可能 这样……，随便扔是吧，这都有可能，什么情况都 有可能。有的针就没有跟平行线相交，比如这个， 这个，这个，就没有相交，也有相交的，比如这个， 这个，这个，这是相交的，对吧，然后他就数，他 说这个针一共投了多少个呢？一共投了n =2212个。 其中与这个平行线相交的针有多少 个，数了一下有m =704个。然后他说， 我现在可以计算圆周率了，别人都不 信，他说你看我圆周率怎么算，我只 要把这两个数相除就行了。我用n 除 以m ，这个数除完了大概是3.142，这个就是圆周率了。别人说好神奇，这怎么回事儿，蒲丰说我给你解释解释这个原理是什么？其实这个原理并不复杂，我们来看一下它的原理是什么。 3. 蒲丰投针原理 (1)首先，它这个平行线是严格平行的，那平行线之间的距离是固定的，是a 。然后我随意地把一根针投上去，也许相交，也许不相交，这不一定。比如说这个针投上去了，投上去了之后，针的总长是b ，针有一个中点叫M ，对吧，这个M 到它比较近的平行线之间的距离我们设为x ，大家注意，这个是针的中点到比较近的平行线的距离是x ，所以我们应该知道x 的范围。x 的最小值就是这个终点正好落在平行线上，那最小值是0，对吧。最大值就 是针的中点正好在两条平行线中间，那最大值是a 2 ，不会再大了。因为我这个x 的定义是针的终点到比较近的平行线的距离，对吧！所以x ∈[0，a 2 ]。 (2)其次就是我想知道这个针与这个平行线的夹角是多少？令夹角为α，α的范围是什么呢，如果你完全跟这个平行线平行的话，那么这个夹角是00，对吧。如果你往上竖过来， 蒲丰氏投针问题的模拟过程，随机数发生器也是自编的，以供大家参考和提出建议。谢谢。（seed1和seed2最好选择3和5，为了使投针次数达到1000000，CVF进行如下设置Project->settings->link-> output，将stack allocations reserve:设为1000000000） program getpi implicit none real,parameter::a=5,L=4,pi=3.14159 integer::n1,i,counter=0 real,allocatable::R1(:),R2(:) real::theta,x,pi1 write(*,*) 'input the size of the array:' read(*,*) n1 allocate(R1(n1)) allocate(R2(n1)) call random(n1,R1,R2) do i=1,n1 x=a*(2*R1(i)-1) theta=pi*R2(i) if(abs(x) 概率论与数理统计实验 蒲丰投针与蒙特卡罗法 班级应数12级01班 学号2012444086 姓名张旭东 蒲丰投针与蒙特卡罗法 张旭东2012444086 （重庆科技学院数学与应用数学，重庆沙坪坝） 【摘要】通过设计一个投针实验使这个事件的概率和未知量π有关，然后通过重复实验，以频率估计概率，即可求得未知参数π的近似解。这种方法称为随机模拟法，也称为蒙特卡罗法。一般来说，实验次数越多所得的近似值就越接近真值。可以利用MATLAB来大量重复地模拟所设计的随机实验。 【关键词】随机模拟；投针实验；重复实验 1 引言 蒲丰投针问题是由法国科学家蒲丰（Buffon）在1777年提出的，它是概率中非常有代表性的问题，它是第一个用几何形式表达概率问题的例子，其结论具有很强的理论与实际意义。蒲丰针问题的解决不仅较典型的反应了集合概率的特征及处理方法，而且还可以由此领略到从“概率土壤”上开出的一朵瑰丽的鲜花——蒙特卡洛（Monte-Carlo）方法。 蒙特卡罗(Monte Carlo)方法，也称计算机模拟方法，是一种基于“随机数”的计算方法，大数定律为近年来发展迅速的随机计算机和随机模拟方法提供了理论基础。 MATLAB是一个适合多学科，具有多种工作平台的功能强大的大型软件。MATLAB已经成为线性代数、自动控制理论、数理统计、数字信号处理、时间序列分析、动态系统仿真等高级课程的进本教学工具，Matlab随机数发生器的种类丰富且用法简便。 本文介绍了利用随机模拟方法和大数定律的相关理论解决蒲丰投针问题计算π的近似值。 2 有关数学实验的有关基础 定理（贝努力大数定律） 设n μ是n 重贝努力实验中事件A 出现的次数，P 是事件A 每次实验中出现的概率，即P(A)=p,则对任意的 ε>0，有 3 实验 蒲丰投针问题 在平面上画有等距离的一些平行线，平行线间的距离为a(a>0)，向平面上随机投一长为l(l Buffon投针实验的理论证明 我们知道，当正多边形的边数无限增多时，它的极限是圆。所以“圆”这种图形可以代表弯曲得最厉害的小针。现在假定圆形小针的直径恰好与纸上两条相邻的平行线间的距离相等，那末这个圆形小针投掷下来时，不是和一条直线相交两次，就是和两条相邻的平行线相切。不管怎样，它的相交次数是2。因此，当投掷的次数为n时，碰线的次数便是2n。 现在小针的长度只有两条相邻平行线间距离的一半，所以针的长度只有上述圆形小针长度（即圆周长）的。但是可能碰线的次数是与针的长度成正比的，因此小针的可能碰线的次数k就必须满足下面的比例式： 1:(1/2π) ＝2n: k 于是就得到π＝n/k，也就是 π＝投掷总次数/碰线次数 这就是上面“投针实验”的理论根据。它又叫莆丰氏实验，在概率论中是很出名的，也可以说是近代的“统计试验法”（又叫“蒙特卡罗法”）的滥觞。 蒲丰(Buffon)投针求π 蒲丰(Buffon)投针问题：在平面上画有等距离的一些平行线，平行线间的距离为a(a>0),向平面上随机投一长为l (l 我们也可以来做这个实验，而且希望做更多次，但是投针又比较费时费力，于是，可以采用另一种设计随机实验的方法，随机模拟的办法来模拟蒲丰投针实验。从而求得π的近似值。 二、实验方法 可以采用MatLAB软件进行模拟实验，即用MATLAB编写程序来进行“蒲丰投针实验”。 1.基本原理：由于针投到纸上的时候，有各种不同的方向和位置（图a），但是，每一次投针时，其位置和方向都可以由两个量唯一确定，那就是针的中点和偏离水平的角度（图b）。 随机投针图 概率模型的随机模拟与蒲丰投针实验 第1章模拟 1.1 模拟的概念 每一个现实系统外部环境之间都存在着一定的数学的或者逻辑的关系，这些关系在系统内部的各个组成部分之间也存在。对数学、逻辑关系并不复杂的模型，人们一般都可用解析论证和数值计算求解。但是，许多现实系统的这种数学、逻辑模型十分复杂，例如大多数具有随机因素的复杂系统。这些系统中的随机性因素很多，一些因素很难甚至不可以用准确的数学公式表述，从而无法对整个系统采用数学解析法求解。这类实际问题往往可以用模拟的方法解决。 模拟主要针对随机系统进行。当然，也可以用于确定性系统。本文讨论的重点是其中的随机模拟。采用模拟技术求解随机模型，往往需要处理大批量的数据。因此，为了加速模拟过程，减少模拟误差，通常借助于计算机进行模拟，因此又称为计算机模拟。 计算机模拟就是在已经建立起的数学、逻辑模型的基础之上，通过计算机试验，对一个系统按照一定的决策原则或作业规则，由一个状态变换为另一个状态的行为进行描述和分析。 1.2 模拟的步骤 整个模拟过程可以划分为一定的阶段，分步骤进行。 (1)明确问题，建立模型。 在进行模拟之前，首先必须正确地描述待研究的问题，明确规定模拟的目的和任务。确定衡量系统性能或模拟输出结果的目标函数，然后根据系统的结构及作业规则，分析系统各状态变量之间的关系，以此为基础建立所研究的系统模型。为了能够正确反映实际问题的本质，可先以影响系统状态发生变化的主要因素建立较为简单的模型，以后再逐步补充和完善。 (2)收集和整理数据资料。 模拟技术的正确运用，往往要大量的输入数据。在随机模拟中，随机数据仅靠一些观察值是不够的。应当对具体收集到的随机性数据资料进行认真分析。确定系统中随机性因素的概率分布特性，以此为依据产生模拟过程所必需的抽样数 《计算方法》实验报告 指导教师： 学院： 班级： 团队成员： 一、题目 例2.7应用Newton 迭代法求方程210x x --=在1x =附近的数值解 k x ，并使其满足8110k k x x ---< 原理： 在方程()0f x =解的隔离区间[],a b 上选取合适的迭代初值0x ，过曲线()y f x =的点()() 00x f x ，引切线 ()()()1000:'l y f x f x x x =+- 其与x 轴相交于点：()() 0100 'f x x x f x =-，进一步,过曲线()y f x =的 点()()11x f x , 引切线 ()()()2111: 'l y f x f x x x =+- 其与x 轴相交于点：() () 1211 'f x x x f x =- 如此循环往复，可得一列逼近方程()0f x =精确解*x 的点 01k x x x ,,,,,其一般表达式为： ()() 111 'k k k k f x x x f x ---=- 该公式所表述的求解方法称为Newton 迭代法或切线法。 程序： function y=f(x)%定义原函数 y=x^3-x-1; end function y1=f1(x0)%求导函数在x0点的值 syms x; t=diff(f(x),x); y1=subs(t,x,x0); end function newton_iteration(x0,tol)%输入初始迭代点x0及精度tol x1=x0-f(x0)/f1(x0);k=1;%调用f函数和f1函数 while abs(x1-x0)>=tol x0=x1;x1=x0-f(x0)/f1(x0);k=k+1; end fprintf('满足精度要求的数值为x(%d)=%1.16g\n',k,x1); fprintf('迭代次数为k=%d\n',k); end 结果： 投针实验计算圆周率的数学分析 王向东 投针实验计算圆周率的数学证明方法，初中一般是采取假设针弯成直径等于平行线距离的方法巧妙证明。这个方法是基于不管针弯成什么形状，针上的每一个部位与平行线相交的概率相同，但这是感观上的认识，要把其中原因解释清楚不是很容易。笔者从纯数学的角度来推导这个公式。 一、投针问题的由来 1777年法国科学家布丰提出的一种计算圆周率的方法——随机投针法，即著名的蒲丰投针问题。 这一方法的步骤是： 1) 取一张白纸，在上面画上许多条间距为d 的平行线。 2) 取一根长度为()l l d <的针，随机地向画有平行直线的纸上掷n 次，观察针与直线相交的次数，记为m 3）计算针与直线相交的概率． 18世纪，法国数学家布丰和勒可莱尔提出的“投针问题”，记载于布丰1777年出版的著作中：“在平面上画有一组间距为d 的平行线，将一根长度为()l l d <的针任意掷在这个平面上，求此针与平行线中任一条相交的概率。”布丰本人证明了，这个概率是： 2l p d π=，π为圆周率。 二、投针实验的数学证明 投针这个动作是由两个事件构成的。 事件1：针投下后与平行线构成一定的夹角。 我们来分析一下针投下后与平行线之间的成某一特定夹角时的概率。 设针投下后与平行线之间的夹角为θ，则θ在0与π之间。针与平行线之间的夹角在θ到θ+θ?之间的概率为1p θ π?=，当0θ?→时，可看作针投下后与平行线之 间成某一特定夹角为θ的概率。 事件2：针投下后会在平行线垂直的方向形成一个投影，针与平行线相交等于它的垂直投影与平行线相交。这个投影的长度'l 在0到l 之间。 function [x0,k]=bisect1(fun1,a,b,ep) if nargin<4 ep=1e-5; end fa=feval(fun1,a); fb=feval(fun1,b); if fa*fb>0 x0=[fa,fb]; k=0; return; end k=1; while abs(b-a)/2>ep x=(a+b)/2; fx=feval(fun1,x); if fx*fa<0 b=x; fb=fx; else a=x; fa=fx; end end x0=(a+b)/2; >> fun1=inline('x^3-x-1'); >> [x0,k]=bisect1(fun1,1.3,1.4,1e-4) x0 = 1.3247 k = 7 >> 简单迭代法 function [x0,k]=iterate1(fun1,x0,ep,N) if nargin<4 N=500; end if nargin<3 ep=1e-5; end x=x0; x0=x+2*ep; while abs(x-x0)>ep & k 计算方法作业（作者：夏云木子） 1、help linspace type linspace 2、a1=[5 12 47;13 41 2;9 6 71];a2=[12 9;6 15;7 21];B=a1*a2, C=a1(:,1:2).*a2, D=a1.^2, E=a1(:).^2 3、a1=[5 12 47;13 41 2;9 6 71];a2=[12 9;6 15;7 21];a1(4:5,1:3)=a2.';a1([4 5],:)=a1([5 4],:);b1=a1 c1=b1(4,1),c2=b1(5,3),D=b1(3:4,:)*a2 4、a1=[5 12 47;13 41 2;9 6 71]; E=eye(3,3); S = a1 + 5*a1' - E, S1=a1^3-rot90(a1)^2+6*E 5、a1=[5 12 47;13 41 2;9 6 71];s=5;A=s-a1,B=s*a1,C=s.*a1,D=s./a1,E=a1./s 6、c=[1 2 3 4;5 6 7 8;9 10 11 12;13 14 15 16];A=c^-4,B=(c^3)^-1,C=(3*c+5*c^-1)/5 7、a=[1 i 3;9i 2-i 8;7 4 8+i];A=a.' 8、abc=[-2.57 8.87;-0.57 3.2-5.5i];m1=sign(abc),m2=round(abc),m3=floor(abc) Sign为符号函数，round表示四舍五入取整，floor表示舍去小数部分取整 9、x=[1 4 3 2 0 8 10 5]';y=[8 0 0 4 2 1 9 11]';A=dot(x,y) 10、a=[3.82 5.71 9.62];b=[7.31 6.42 2.48];A=dot(a,b),B=cross(a,b) 11、P=[5 7 8 0 1];Pf=poly(P);Px=poly2str(Pf,'x') 12、P=[3 0 9 60 0 -90];K1=polyval(P,45),K2=polyval(P,-123),K3=polyval(P,579) 13、P1=[13 55 0 -17 9];P2=[63 0 26 -85 0 105];PP=conv(P1,P2);P1P2=poly2str(PP,'x'),[Q,r]=deconv(P2,P1) 蒙特卡罗方法概述 § 8.2 引例：蒲丰投针问题 在用传统方法难以解决的问题中，有很大一部分可以用概率模型进行描述．由于这类模型含有不确定的随机因素，分析起来通常比确定性的模型困难．有的模型难以作定量分析，得不到解析的结果，或者是虽有解析结果，但计算代价太大以至不能使用．在这种情况下，可以考虑采用Monte Carlo 方法。下面通过例子简单介绍Monte Carlo 方法的基本思想． Monte Carlo 方法是计算机模拟的基础，它的名字来源于世界著名的赌城——摩纳哥的蒙特卡洛，其历史起源于1777年法国科学家蒲丰提出的一种计算圆周π的方法——随机投针法，即著名的蒲丰投针问题。这一方法的步骤是： 1) 1) 取一张白纸，在上面画上许多条间距为d 的平行线，见图8.1(1) 2) 2) 取一根长度为)(d l l <的针，随机地向画有平行直线的纸上掷n 次，观察针与直线相交的次数，记为 m 3）计算针与直线相交的概率． 由分析知针与平行线相交的充要条件是 ?sin 21≤ x 其中 π?≤≤≤≤0,2 0d x 建立直角坐标系),(x ?,上述条件在坐标系下将是曲线所围成的曲边梯形区域，见图 8．l （2）． 由几何概率知 (*)22 sin 210d l d d G g p ππ??π===?的面积的面积 4）经统计实验估计出概率,n m P ≈由(*)式即?2=?=ππd l n m Monte Carlo 方法的基本思想是首先建立一个概率模型，使所求问题的解正好是该模型的参数或其他有关的特征量．然后通过模拟一统计试验，即多次随机抽样试验（确定m 和n ），统计出某事件发生的百分比．只要试验次数很大，该百分比便近似于事件发生的概率．这实际上就是概率的统计定义．利用建立的概率模型，求出要估计的参数．蒙特卡洛方法属于试验数学的一个分支． ************************************************************************* 提示：设x 是一个随机变量，它服从区间[０，d/2]是的均匀分布，同理，?是一个随机变量，它服从区间],0[π上的均匀分布。按照某种抽样法，产生随机变量的可能取值，例如 蒲丰投针――Monte Carlo 算法 背景： 蒙特卡罗方法（Monte Carlo），也称统计模拟方法，是在二次世界大战期间随着科学技术的发展和电子计算机的发明，而被提出的一种以概率统计理论为基础的一类非常重要的数值计算方法。蒙特卡罗方法在应用物理、原子能、固体物理、化学、生态学、社会学以及经济行为等领域中得到广泛利用。 蒙特卡罗方法的名字来源于世界著名的赌城——摩纳哥的蒙特卡罗。其历史起源可追溯到1777年法国科学家蒲丰提出的一种计算圆周的方法——随机投针法，即著名的蒲丰投针问题。 问题： 设在平面上有一组平行线，间距为d，把一 根长L的针随机投上去，则这根针和平行线相交 的概率是多少？(其中L < d ) 分析：由于L < d，所以这根针至多只能与一条平行线相交。设针的中点与最近的平行线之间的距离为y，针与平行线的夹角为θ (0 ≤θ≤π)。 相交情形不相交情形 易知针与平行线相交的充要条件是： sin 2 L y xθ ≤= 由于 1 [0,],[0,] 2 y dθπ ∈∈，且它们的取值均 满足平均分布。建立直角坐标系，则针与平行线 的相交条件在坐标系下就是曲线所围成的曲边梯 形区域（见右图）。所以有几何概率可知针与平行 线相交的概率是 sin d2 2 1 2 L L p d d π θθ π π == ? Monte Carlo 方法： 随机产生满足平均分布的 y 和 θ，其中1 [0, ], [0, ]2 y d θπ∈∈，判断 y 是否在曲边梯形内。重复上述试验，并统计 y 在曲边梯形内的次数 m ，其与试验次数 n 的比值即为针与平行线相交的概率的近似值。 clear; n = 100000; L = 1; d = 2; m = 0; for k = 1 : n theta = rand(1)*pi; y = rand(1)*d/2; if y < sin(theta)*L/2 m = m + 1; end end fprintf('针与平行线相交的概率大约为 %f\n', m/n) 计算π的近似值 利用该方法可以计算 π 的近似值： sin d 22 2 2 1n L L m p d m d L d n π θθπππ?≈= =≈? 下面是一些通过蒲丰投针实验计算出来的 π 的近似值： 蒲丰投针问题的重要性并非是为了求得比其它方法更精确的π值，而是在于它是第一个用几何形式表达概率问题的例子。计算π的这一方法，不但因其新颖，奇妙而让人叫绝，而且它开创了使用随机数处理确定性数学问题的先河，是用偶然性方法去解决确定性计算的前导。 1777年法国科学家布丰提出的一种计算圆周率的方法——随机投针法，即著名的蒲丰投针问题。 投针步骤 这一方法的步骤是： 1) 取一张白纸，在上面画上许多条间距为d的平行线。 2) 取一根长度为l（l基于MATLAB的布丰投针实验仿真
《应用计算方法教程》matlab作业二
计算方法_全主元消去法_matlab程序
几何概型的解析方法和类型例解以及蒲丰投针试验、贝特朗奇论的解析
布丰投针实验模拟
(整理)matlab16常用计算方法.
matlab用于计算方法的源程序
蒲丰投针问题
蒲丰氏投针问题的模拟过程
蒲丰投针实验模拟
Buffon投针实验的理论证明
蒲丰投针及蒙特卡罗模拟电子教案模拟
计算方法上机实验报告-MATLAB
投针实验计算圆周率的数学分析
数值计算方法matlab程序
计算方法及其MATLAB实现第一章作业
蒲丰投针问题
蒲丰投针――MonteCarlo算法
布丰投针实验