修改一个蒙特卡洛计算定积分的源程序
高级算法语言期末实习作业,从若干题中抽一个题,我头晕 选了个这个。求各位大神达人们帮帮忙修改下!!!!!!!!!蒙特卡洛计算定积分的源程序!!!!!#define IM1 2147483563
#define IM2 2147483399
#define AM (1.0/IM1)
#define IMM1 (IM1-1)
#define IA1 40014
#define IA2 40692
#define IQ1 53668
#define IQ2 52774
#define IR1 12211
#define IR2 3791
#define NTAB 32
#define NDIV (1+IMM1/NTAB)
#define EPS 1.2e-7
#define RNMX (1.0-EPS) /*以上是为ran2函数做准备的*/
#include <stdlib.h>
#include<stdio.h>
#include<math.h>
#include <time.h>
double ran2(long *idum) /*产生随机数*/
{
int j;
long k;
static long idum2=123456789;
static long iy=0;
static long iv[NTAB];
double temp;
if(*idum<=0)
{
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];
}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)
return RNMX;
else
return temp;
}
double fun1(double x) /*被积函数*/
{
return(sin(x));
}
double Monte(double n,double a,double b,double (*p)(),double fMAX) /*计算蒙特卡洛因子*/
{
double x,y;
double sx=0,sy=0,vx=0,vy=0;
double s,t;
double i;
long ko1,ko2;
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++;
}
printf("\nEnter b=:");scanf("%lf",&b);
printf("\nEnter N=:");scanf("%lf",&N); /*N表示区间的个数*/
printf("\nEnter n=:");scanf("%lf",&n);
printf("\nEnter fMAX(此函数的界)=:");scanf("%lf",&fMAX);
h=(b-a)/N;
time0=clock();
else
{
if(y>=(*p)(x))
vx++;
}
}
return ((sx-vx)/n);
}
main()
{
double n,fMAX,a,b,Mont,A,B,h,i=0,N;
long time0,time1;
printf("Enter a=:");scanf("%lf",&a);
for(i=0;i<N;i++)
{
A=a+i*h;
B=A+h;
Mont+=Monte(n,A,B,fun1,fMAX)*(B-A)*2*fMAX;
}
printf("\nintegration=%lf",Mont);
time1=clock();
printf("\n耗时:%ldms",time1-time0);
getch();
}
[ 本帖最后由 jxh41321 于 2014-1-10 13:06 编辑 ]