Calculate the integral using Monte Carlo method

#include "stdio.h"
#include "stdlib.h"
#include "math.h"

double rnd1(double *r) //the definition body of rnd1
{ int m;
double s, u, v, p;
s=65536.0; u=2053.0; v=13849.0;
m=(int)(*r/s); *r=*r-m*s;
*r=u*(*r)+v; m=(int)(*r/s);
*r=*r-m*s; p=*r/s;
return(p);
}

double mtml(int n, double a[], double b[], double (*f)()) /*n为积分的维度,数组a保存着积分下限,
数组b保存着积分上限,f为函数指针,指向被积函数*/
{
int m, i;
double r, s, d, *x;
x=malloc(n*sizeof(double));
r=1.0; d=10000.0; s=0.0;
for (m=0; m<=9999; m++)
{
for (i=0; i<=2-1; i++)
x[i]=a[i]+(b[i]-a[i])*rnd1(&r);
s=s+(*f)(n, x)/d;
}
for (i=0; i<=n-1; i++)
s=s*(b[i]-a[i]);
free(x);
return(s);
}


main()
{
static double a[3]={1.0, 1.0, 1.0};
static double b[3]={2.0, 2.0, 2.0};
double mtmlf(int, double []);
printf("\n");
printf("s=%e\n", mtml(3, a, b, mtmlf));
printf("\n");
}

double mtmlf(int n, double x[])//mtmlf为被积函数,这里是定值1
{
int i;
double f;
f=0.0;
for(i=0; i<=n-1; i++)
f=f+x[i]*x[i];
return(f);
}
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值