| 网站首页 | 业界新闻 | 小组 | 威客 | 人才 | 下载频道 | 博客 | 代码贴 | 在线编程 | 编程论坛
欢迎加入我们,一同切磋技术
用户名:   
 
密 码:  
共有 533 人关注过本帖
标题:修改一个蒙特卡洛计算定积分的源程序
只看楼主 加入收藏
jxh41321
Rank: 1
等 级:新手上路
帖 子:2
专家分:0
注 册:2014-1-10
收藏
 问题点数:0 回复次数:0 
修改一个蒙特卡洛计算定积分的源程序
高级算法语言期末实习作业,从若干题中抽一个题,我头晕 选了个这个。求各位大神达人们帮帮忙修改下!!!!!!!!!蒙特卡洛计算定积分的源程序!!!!!
#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 编辑 ]
搜索更多相关主题的帖子: 蒙特卡洛 include 源程序 神达 
2014-01-10 13:03
快速回复:修改一个蒙特卡洛计算定积分的源程序
数据加载中...
 
   



关于我们 | 广告合作 | 编程中国 | 清除Cookies | TOP | 手机版

编程中国 版权所有,并保留所有权利。
Powered by Discuz, Processed in 0.066777 second(s), 8 queries.
Copyright©2004-2024, BCCN.NET, All Rights Reserved