蒲丰投针及蒙特卡罗模拟电子教案模拟
概率模型的随机模拟与蒲丰投针实验 第1章模拟 1.1 模拟的概念 每一个现实系统外部环境之间都存在着一定的数学的或者逻辑的关系,这些关系在系统内部的各个组成部分之间也存在。对数学、逻辑关系并不复杂的模型,人们一般都可用解析论证和数值计算求解。但是,许多现实系统的这种数学、逻辑模型十分复杂,例如大多数具有随机因素的复杂系统。这些系统中的随机性因素很多,一些因素很难甚至不可以用准确的数学公式表述,从而无法对整个系统采用数学解析法求解。这类实际问题往往可以用模拟的方法解决。 模拟主要针对随机系统进行。当然,也可以用于确定性系统。本文讨论的重点是其中的随机模拟。采用模拟技术求解随机模型,往往需要处理大批量的数据。因此,为了加速模拟过程,减少模拟误差,通常借助于计算机进行模拟,因此又称为计算机模拟。 计算机模拟就是在已经建立起的数学、逻辑模型的基础之上,通过计算机试验,对一个系统按照一定的决策原则或作业规则,由一个状态变换为另一个状态的行为进行描述和分析。 1.2 模拟的步骤 整个模拟过程可以划分为一定的阶段,分步骤进行。 (1)明确问题,建立模型。 在进行模拟之前,首先必须正确地描述待研究的问题,明确规定模拟的目的和任务。确定衡量系统性能或模拟输出结果的目标函数,然后根据系统的结构及作业规则,分析系统各状态变量之间的关系,以此为基础建立所研究的系统模型。为了能够正确反映实际问题的本质,可先以影响系统状态发生变化的主要因素建立较为简单的模型,以后再逐步补充和完善。 (2)收集和整理数据资料。 模拟技术的正确运用,往往要大量的输入数据。在随机模拟中,随机数据仅靠一些观察值是不够的。应当对具体收集到的随机性数据资料进行认真分析。确定系统中随机性因素的概率分布特性,以此为依据产生模拟过程所必需的抽样数
计算方法上机实验报告-MATLAB
《计算方法》实验报告 指导教师: 学院: 班级: 团队成员:
一、题目 例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 结果:
蒲丰投针问题
蒙特卡罗方法概述 § 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 是一个随机变量,它服从区间[0,d/2]是的均匀分布,同理,?是一个随机变量,它服从区间],0[π上的均匀分布。按照某种抽样法,产生随机变量的可能取值,例如
数值计算方法matlab程序
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> fun1=inline('(x+1)^(1/3)'); >> [x0,k]=iterate1(fun1,1.5) x0 = 1.3247 k = 7 >> fun1=inline('x^3-1'); >> [x0,k]=iterate1(fun1,1.5) x0 = Inf k =
计算方法及其MATLAB实现第一章作业
计算方法作业(作者:夏云木子) 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)
蒲丰投针――MonteCarlo算法
蒲丰投针――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 年的一天,法国科学家布丰的家里宾客满堂,原来他们是应主人的邀请前来观看一次奇特试验的。 试验开始,但见年已古稀的布丰先生兴致勃勃地拿出一张纸来,纸上预先 画好了一条条等距离的平行线。接着他又抓出一大把原先准备好的小针。然后 布丰先生宣布:“请诸位把这些小针一根一根往纸上扔吧!不过,请大家务必 把扔下的针是否与纸上的平行线相交,以及相交的次数告诉我。 客人们不知布丰先生要玩什么把戏,只好客随主意,一个个加入了试验的 行列。一把小针扔完了,把它捡起来再扔。而布丰先生本人则不停地在一旁数着、记着,如此这般地忙碌了将近一个钟头。最后,布丰先生高声宣布:“先 生们,我这里记录了诸位刚才的投针结果,共投针 2212 次,其中与平行线相交的有 704 次。总次数 2212 与相交次数 704 的比值为 3.142。”说到这里,布丰先生故意停了停,并对大家报以神秘的一笑,接着有意提高声调说:“先生们,这就是圆周率π的近似值!” 客人们一片哗然,议论纷纷,大家全都感到莫名其妙:“圆周率π?这可跟投针半点也不沾边呀!” 布丰先生似乎猜透了大家的心思,得意洋洋地解释道:“诸位,这里用的是概率的原理,如果大家有耐心的话,再增加投针的次数,还能得到π 的更精确的近似值呢。” 那么,“布丰投针实验”的依据究竟是什么呢?下面就是书中简单而巧妙 的证明。为了便于理解,我把证明过程说得稍微详细一点。 假设那组平行线的间距等于 d。如果把一个直径为 d 的铁丝圆圈,扔到平行线组上,因为它的周长等于πd,所以,不论怎样扔,每个圆圈都会与平行线有两个交点。因此,如果扔下的次数为 n,交点的总数为 m,必定有 m=2n。 还用那组平行线,不过这回把圆圈剪开拉直,变成长度为πd的直铁丝。显然,直铁丝与平行线相交的情形要比圆圈复杂,最多可能有 4 个交点,也可能有 3 个、2 个、1 个交点,也可能不相交,没有交点。不过,由于圆圈和直铁
布丰投针实验
1777年法国科学家布丰提出的一种计算圆周率的方法——随机投针法,即著名的蒲丰投针问题。 投针步骤 这一方法的步骤是: 1) 取一张白纸,在上面画上许多条间距为d的平行线。 2) 取一根长度为l(lz,x²+y²﹤z²,容易证明这两个式子即为以这3个正数为边长可以围成一个钝角三角形的充要条件,用线性规划可知满足题设的可行域为直线x+y=z与圆x²+y²=z²围成的弓形,总的可行域为一个边长为z的正方形,则可
matlab计算
基本代数部分 如何用matlab求阶乘 factorial(n)求n的阶乘 如何用matlab配方 没有发现matlab有这一命令,不过我们可以调用maple的命令,调用方法如下: 首先加载maple中的student函数库,加载方法为:maple(‘with(student)’) 然后运行maple中的配方命令,格式为: maple(’completesquare(f)’)把f配方,其中f为代数表达式或代数方程maple(’completesquare(f,x)’)把f按指定的变量x配方,其中f同上maple(’completesquare(f,{x,y,...})’)把f按指定的变量x,y,...配方 maple(’completesquare(f,[x,y,...])’)把f按指定的变量x,y,...配方, 如何用matlab进行多项式运算 (1)合并同类项 syms 表达式中包含的变量 collect(表达式,指定的变量) (2)因式分解 syms 表达式中包含的变量factor(表达式) (3)展开
syms 表达式中包含的变量 expand(表达式) (4)化简 syms 表达式中包含的变量simplify(表达式) (5)变量替换 syms 表达式和代换式中包含的所有变量subs(表达式,要替换的变量或式子,代换式)我们也可在matlab中调用maple的命令进行多项式的运算,如下: 如何在matlab中调用maple 方法1: maple(’maplestatement’) 其中maplestatement 是完整的maple语句,由一条或几条命令组成,必须符合maple 的语法 方法2: maple(‘function’,arg1, arg2,…) 其中function为maple中的函数名称,arg1, arg2,…是函数function所用的参数。 如何用matlab进行分式运算 发现matlab只有一条处理分式问题的命令,其使用格式如下: [n,d]=numden(f)把符号表达式f化简为有理形式,其中分子和分母的系数为整数且分子分母不含公约项,返回结果n为分子,d为分母。注意:f必须为符号表达式
用递推公式计算定积分matlab版
用递推公式计算定积分 实验目的: 1.充分理解不稳定的计算方法会造成误差的积累,在计算过程中会导致误差的迅速增加,从而使结果产生较大的误差。 2.在选择数值计算公式来进行近似计算时,应学会选用那些在计算过程中不会导致误差迅速增长的计算公式。 3.理解不稳定的计算公式造成误差积累的来源及具体过程; 4.掌握简单的matlab语言进行数值计算的方法。 实验题目: 对n=0,1,2,…,20,计算定积分: 实验原理: 由于y(n)== – 在计算时有两种迭代方法,如下: 方法一: y(n)=– 5*y(n-1),n=1,2,3, (20) 取y(0)== ln6-ln5 ≈0.182322 方法二: 利用递推公式:y(n-1)=-*y(n),n=20,19, (1) 而且,由= * ≤≤*=
可取:y(20)≈*()≈0.008730. 实验内容: 对算法一,程序代码如下: function [y,n]=funa() syms k n t; t=0.182322; n=0; y=zeros(1,20); y(1)=t; for k=2:20 y(k)=1/k-5*y(k-1); n=n+1; end y(1:6) y(7:11) 对算法二,程序代码如下: %计算定积分; %n--表示迭代次数; %y用来存储结果; function [y,n]=f(); syms k y_20;
y=zeros(21,1); n=1; y_20=(1/105+1/126)/2; y(21)=y_20; for k=21:-1:2 y(k-1)=1/(5*(k-1))-y(k)/5; n=n+1; end 实验结果: 由于计算过程中,前11个数字太小,后9个数字比较大,造成前面几个数字只显示0.0000的现象,所以先输出前6个,再输出7—11个,这样就能全部显示出来了。 算法一结果: [y,n]=funa %先显示一y(1)—y(6) ans = 0.1823 -0.4116 2.3914 -11.7069 58.7346