(08063063)(林豪)(Monte-Carlo方法在定积分中的应用)

数学与统计学院课程设计报告课程：数值分析题目：Monte-Carlo方法在数值积分中的应用年级：2008专业：数学与应用数学学号：08063063姓名：林豪指导教师：宁娣2010年12月12日数学与统计学院本科课程设计Monte-Carlo方法在数值积分中的应用【摘要】本文分析了Monte-Carlo方法在数值积分计算中的优缺点，发现该方法的误差阶仅为1On，为了改善误差精度，提出了改进的Monte-Carlo方法——“类矩形”法和“类梯形”法，其误差阶分别为1On和21On，并给出了两种方法的一般计算步骤，最后通过实例验证了其误差精度的确如上所述。【关键词】数值积分；Monte-Carlo；类矩形法；类梯形法；误差阶一、Monte-Carlo方法简介Monte-Carlo方法又称随机模拟方法或统计试验方法，它是一类通过随机变量的统计试验，随机模拟求解数学物理，工程技术问题近似解的数值方法。随着科学及计算机技术的发展，该方法已在求解数值积分，积分方程，微分方程，非线性方程组问题，计算物理，大型系统可靠性分析等方面得到了广泛的应用。最常见的计算定积分的Monte-Carlo方法是平均值法，该方法计算简单，适用范围广，但在某些要求精度非常高的场合不能取得理想的效果。为了提高精度，本文结合矩形求积公式和梯形求积公式，给出了Monte-Carlo计算积分的一类改进算法。二、Monte-Carlo的原始平均值法1、Monte-Carlo法的计算步骤已知()fx在[,]ab上可积，试用平均值法计算定积分()baIfxdx。基于上述Monte-Carlo的思想步骤如下：1）随机产生n个服从[,]ab上均匀分布的独立随机变量(1,2,,)iin；2）计算1()niibaIfn，并将其作为I的近似值，即II。2、误差分析由中心极限定理知，若(1,2,,)iin相互独立、同分布且数学期望及标准差0存在，则当n充分大时，随机变量IIYn服从正太分布(0,1)N，即对任意的0t，220212xttPYtPIIedxn这表明，用平均值法计算定积分的收敛速度较慢，在概率意义下的误差仅为1On，为了提高精度，下面对平均值法做如下两种改进。三、平均值法的改进1、“类矩形”的Monte-Carlo方法简介先将区间[,]abn等分，再在每个区间上各产生一个随机点，然后由着n个随机点类似于矩形公式构造计算公式，次称为“类矩形”的Monte-Carlo方法，如此处理计算量几乎没有变化，但误差阶却提高到1On。2、“类矩形”的计算步骤为方便见，仍以计算()baIfxdx为例，其计算步骤如下：1）随机产生n个相互独立且服从[0,1]上均匀分布的随机变量序列(1,2,,)iin；2）作变换(1),1,2,iibaaiinn将i映射到子区间1[(),()][,],1,2,,iiabaabaabinnn3）计算1()niibaIfn并将其作为I的近似值。四、平均值法的进一步改进1、“类矩形”的改进——“类梯形”法如上，通过将区间等分，在每个子区间上任取一个随机点，类似于矩形公式计算，得到了改进的“类矩形’法。进一步地，若先将[,]abn等分，再在每个子区间上各产生两个随机点，然后类似于梯形公式构造计算公式，便可得到改进的“类梯形”法。其计算精度可进一步提高，误差阶可提高到21On2、计算步骤1）随机产生2n个相互独立且服从[0,1]上均匀分布的随机变量序列，并两两分组得：212,,(1,2,,)iiin；2）做变换2121(22)2iibaain22(21)2iibaain将212,ii分布映射到子区间121[(),()]2iiabaabann211[(),()]2iiabaabann其中1,2,,in。3）在每个等分子区间[(1),]babaaiainn上利用212,ii两点类似于梯形公式构造“类梯形”公式。212()()()2iiffSban来近似(1)()aihaihfxdx4）计算2121()()2niiffbaIn将其作为I的近似值。五、实例分析以20xIedx为例，分别用三种方法计算该积分。得出结果用下表列出：表1：三种方法结果对比nMonte-Carlo“类矩形”“类梯形”准确值100.779189830.893413130.880806190.886219051000.899911660.885627270.886223140.8862190510000.877536690.886246410.886219160.88621905100000.885652680.886254840.886219050.88621905六、结束语上述数值试验表明，直接用原始的平均值法计算定积分，410个节点的计算量已经很可观了，但计算结果只有两位有效数字，而选取同样的节点数，计算了几乎不变，“类矩形”法就达到了4位有效数字，而“类梯形”法则达到了8位有效数字，恰好与上述论述中的误差阶的估计是一致的，从而验证了该方法的高效性。参考文献[1]李庆杨，王能超，易大义．数值分析[M]．武汉：华中科技大学，2010.[2]张平文，李铁军.数值分析[M].北京：北京大学出版社，2006.[3]王明兰，叶恒青.计算近似积分的样本均值Monte-Carlo方法[J].华南师范大学学报（自然科学版），2001,1：50-52.附录:一、用Monte-Carlo法计算函数定积分的MATLAB源程序：functionMC_method(a,b,N)%a=input('请输入区间左断点（a）=');%b=input('请输入区间右断点（b）=');%N=input('请输入节点数（N）=');x=a+(b-a)*rand(1,N);S=subs('exp(-x^2)',x);INT=sum(S)*(b-a)/Nvpa(int('exp(-x^2)','x',a,b))二、用“类矩形”法计算函数定积分的MATLAB源程序：functionrectangle_method(a,b,N)%a=input('请输入区间左断点（a）=');%b=input('请输入区间右断点（b）=');%N=input('请输入节点数（N）=');s=rand(1,N);deta=0:N-1;x=(s+deta)*(b-a)/N+a;S=subs('exp(-x^2)',x);INT=sum(S)*(b-a)/Nvpa(int('exp(-x^2)','x',a,b))三、用“类梯形”法计算函数定积分的MATLAB源程序：functiontrapezia_method(a,b,N)%a=input('请输入区间左断点（a）=');%b=input('请输入区间右断点（b）=');%N=input('请输入节点数（N）=');s1=rand(1,N);s2=rand(1,N);deta1=0:N-1;deta2=1:N;x1=(s1+2*deta1)*(b-a)/2/N+a;x2=(s2+deta1+deta2)*(b-a)/2/N+a;S1=subs('exp(-x^2)',x1);S2=subs('exp(-x^2)',x2);INT1=sum(S1)*(b-a)/2/N;INT2=sum(S2)*(b-a)/2/N;INT=INT1+INT2vpa(int('exp(-x^2)','x',a,b))