定积分的MonteCarlo计算方法的实现

《计算机高级语言》

认知实习报告

课题名称：定积分的MonteCarlo计算方法的实现

指导老师：王玉兰小组成员：200907020302曾颖超200907020301李海全200807020111向慧

实习时间：2010.12.27——2010.01.08摘要：本次实习设计并实现定积分的MonteCarlo计算方法，采用VC++6.0开发完成。定积分的MonteCarlo计算方法的实现共26页，您现在浏览的是第1页!

一、设计任务与要求1、设计内容2、设计要求二、MonteCarlo算法引导1、MonteCarlo算法的概念2、MonteCarlo算法的例子三、MonteCarlo思想计算定积分的实现1、定积分的定义2、MonteCarlo算法的积分原理四、详细设计及实现五、总结六、参考资料附录：源程序定积分的MonteCarlo计算方法的实现共26页，您现在浏览的是第2页!一、设计任务与要求1、设计内容针对某具体函数f(x)形式，利用MonteCarlo算法实现积分。定积分的MonteCarlo计算方法的实现共26页，您现在浏览的是第3页!二、MonteCarlo算法引导

1、MonteCarlo算法的概念

非形式化地说，MonteCarlo算法泛指一类算法。在这些算法中，要求解的问题是某随机事件的概率或某随机变量的期望。这时，通过“实验”方法，用频率代替概率或得到随机变量的某些数字特征，以此作为问题的解。定积分的MonteCarlo计算方法的实现共26页，您现在浏览的是第4页!三、Monte-Carlo思想计算定积分的实现1、定积分的定义2、MonteCarlo积分法的原理

定积分的MonteCarlo计算方法的实现共26页，您现在浏览的是第5页!如果无限接近于（亦即）时，上述和式无限趋近于常数，那么称该常数为函数在区间上的定积分。记为：其中成为被积函数，叫做积分变量，为积分区间，积分上限，积分下限。定积分的MonteCarlo计算方法的实现共26页，您现在浏览的是第6页!四、详细设计及实现根据蒙特卡洛算法的概念，将落在数学函数f(x)里面的随机点收集起来并计算出数量，这样可以计算出蒙特卡洛因子。然后就可计算出函数f(x)的积分。我们把积分区间较长的函数进行分段求解，然后再求和。希望以此来提高精确性和缩短时间。定积分的MonteCarlo计算方法的实现共26页，您现在浏览的是第7页!nN积分误差（取正）时间/毫秒10610.566960.02495723810620.58220.00971746510640.592270.00035393110680.5909450.0009721,813106160.595310.0033933,566106320.5951810.0032647,179定积分的MonteCarlo计算方法的实现共26页，您现在浏览的是第8页!六总结用蒙特卡洛方法计算定积分具有普遍意义。根据以上试验我们了解到，对于求解一些非常规的问题我们有了新思路，即一个不具随机性的事件可以通过一定的方法用随机事件来模拟或逼近。定积分的MonteCarlo计算方法的实现共26页，您现在浏览的是第9页!附录：源程序#defineIM12147483563#defineIM22147483399#defineAM(1.0/IM1)#defineIMM1(IM1-1)#defineIA140014#defineIA240692#defineIQ153668#defineIQ252774#defineIR112211#defineIR23791#defineNTAB32#defineNDIV(1+IMM1/NTAB)定积分的MonteCarlo计算方法的实现共26页，您现在浏览的是第10页!doubleran2(long*idum)/*产生随机数*/{intj;longk;staticlongidum2=123456789;staticlongiy=0;staticlongiv[NTAB];doubletemp;

if(*idum<=0)定积分的MonteCarlo计算方法的实现共26页，您现在浏览的是第11页!k=(*idum)/IQ1;*idum=IA1*(*idum-k*IQ1)-k*IR1;if(*idum<0)*idum+=IM1;k=idum2/IQ2;idum2=IA2*(idum2-k*IQ2)-k*IR2;if(idum2<0)idum2+=IM2;j=iy/NDIV;iy=iv[j]-idum2;iv[j]=*idum;if(iy<1)iy+=IMM1;if((temp=AM*iy)>RNMX)returnRNMX;elsereturntemp;}定积分的MonteCarlo计算方法的实现共26页，您现在浏览的是第12页!ko1=ko2=1;for(i=0;i<n;i++){t=ran2(&ko1);s=ran2(&ko2);x=a+(b-a)*t;y=(-1)*fMAX+2*fMAX*s;if(y>=0){if(y==0)sy++;if(y<=(*p)(x))sx++;}定积分的MonteCarlo计算方法的实现共26页，您现在浏览的是第13页!printf("\nEnterb=:");scanf("%lf",&b); printf("\nEnterN=:");scanf("%lf",&N);/*N表示区间的个数*/printf("\nEntern=:");scanf("%lf",&n);printf("\nEnterfMAX(此函数的界)=:");scanf("%lf",&fMAX);

h=(b-a)/N;time0=clock();定积分的MonteCarlo计算方法的实现共26页，您现在浏览的是第14页!2、设计要求1)学习MonteCarlo计算数值的积分方法；2)就某个具体的被积函数f(x)的定积分，设计算法步骤;3)编写程序实现算法；4）对程序进行调试和测试，使用有解析解的定积分以便计算对比；定积分的MonteCarlo计算方法的实现共26页，您现在浏览的是第15页!2、MonteCarlo算法的例子如图(1)所示，在一个1平方米的正方形木板上，随意画一个圈，求返个圈的面积。图(1)假设我手里有一支飞镖，我将飞镖掷向木板。木板，我们假定每一次都能掷在木板上，不会偏出木板，但每一次掷在木板的什么地方，是完全随机的。即，每一次掷飞镖，飞镖扎进木板的任何一点的概率的相等的。返样，我们投掷多次，例如100次，然后我们统计返100次中，扎入不规则图形内部的次数，假设为k，那么，我们就可以用k/100*1近似估计不规则图形的面积，例如100次有32次掷入图形内，我们就可以估计图形的面积为0.32平方米从上述可以看出，Monte-Carlo算法区别不确定性算法，它的解不一定是准确或正确的，其准确或正确性依赖于概率和统计，但在某些问题上，当重复实验次数越够大时，可从很大概率上（返个概率是可以在数学上证明的，但依赖于具体问题）确保解的准确或正确性，所以，我们可以根据具体的概率分析，设定实验的次数，从而将误差或错误率降到一个可容忍的程度。定积分的MonteCarlo计算方法的实现共26页，您现在浏览的是第16页!1.定积分的定义如图2，定积分就是求函数在区间[a,b]中图线下包围的面积，即y=0,x=a,x=b,y=f(x)所包围的面积。一般地，设函数在区间上连续，用分点将区间等分成个小区间，每个小区间长度为（），在每个小区间上取一点，作和式：定积分的MonteCarlo计算方法的实现共26页，您现在浏览的是第17页!2、MonteCarlo积分法的原理

求定积分的直观意义就是求面积，所以，用MonteCarlo求积分的原理就是通过模拟统计方法求解面积。即通过向特定区域随机产生大量点，然后统计点落在函数区域内的频率，以此频率估计面积，从而得到积分值。定积分的MonteCarlo计算方法的实现共26页，您现在浏览的是第18页!五测试被积函数：f(x)=sin(x),x∈[0,20]积分精确解为：0.591917把区间分为N段，随机点个数为n个以下是测试数据定积分的MonteCarlo计算方法的实现共26页，您现在浏览的是第19页!10710.5986440.0067272,30810720.5960520.0041354,60710740.5977180.0058019,01110780.59310.00118317,847107160.5903590.00155835,352定积分的MonteCarlo计算方法的实现共26页，您现在浏览的是第20页!参考文献WilliamH.Press等著，傅祖芸等译，C数值算法（第二版）电子工业出版社，2004.01同济大学数学系编，高等数学（第六版）高等教育出版社，2007.06孙淑霞李思明刘焕君编著，C/C++程序设计（第三版）电子工业出版社，2009.09定积分的MonteCarlo计算方法的实现共26页，您现在浏览的是第21页!#defineEPS1.2e-7#defineRNMX(1.0-EPS)/*以上是为ran2函数做准备的*/#include<stdlib.h>#include<stdio.h>#include<math.h>#include<time.h>定积分的MonteCarlo计算方法的实现共26页，您现在浏览的是第22页!

{if(-(*idum)<1)*idum=1;else*idum=-(*idum);for(j=NTAB+7;j>=0;j--){k=(*idum)/IQ1;*idum=IA1*(*idum-k*IQ1)-k*IR1;if(*idum<0)*idum+=IM1;if(j<NTAB)iv[j]=*idum;}iy=iv[0];}定积分的MonteCarlo计算方法的实现共26页，您现在浏览的是第23页!doublefun1(doublex)/*被积函数*/{return(sin(x));}doubleMonte(doublen,doublea,doubleb,double(*p)(),doublefMAX)/*计算蒙特卡洛因子*/{doublex,y;doublesx=0,sy=0,vx=0,vy=0;doubles,t;doublei;longko1,ko2;定积分的MonteCarlo计算方法的实现共26页，您现在浏览的是第24页!else{if(y>=(*p)(x))vx++;}}return((sx-vx)/n);}main(){doublen,fMAX,a,b,Mont,A,B,h,i=0,N;longtime0,time1;printf("Entera=:");scanf("%lf",&a);定积分的MonteCarlo计算方法的实现共26页，您现在浏览的是第2