蒙特卡罗洛法检验三重积分
这个主要是作微积分的时候突然大家发现答案不对时的一种检验方法,算得挺准,就是有点废计算机。
那么废话不多说,蒙特卡洛法其实就是一种概率方法,基本思想和三重积分的分割定义也差不多。就是切成好多小方块。
一开始先在一个大区域内撒点,也就是均匀地在一个立方体区域内取点,像异次元杀阵的网格点一样(没看过异次元杀阵的话理解成NaCl晶体模型吧○| ̄|_)
比较方便的是只要包含积分区域就好了,所以还是编起来挺方便的,题目不管规则不规则,只要找到一个超大矩形就可。
来看例题,如果有清华教科书可以看下册的161页7(5),强烈吐槽这题答案错了(其实清华教科书答案错得挺多的)。
题目描述:求∫∫∫Ω xyzdxdydz,Ω={(x,y,z) | x2+y2+z2≤4,x2+y2+(z-2)2≤4,x≥0,y≥0};
那么很显然啦,方程t=xyz,区域嘛就是两个球的交集,很容易就发现[0,2]×[0,2]×[0,2]这个矩形域包括了Ω区间。所以用蒙特卡洛,就是在区域[0,2][0,2][0,2]的区域上随机撒N个均匀分布的点,计算公式是V/N*所有在积分区域内的点的xyz乘积之和。其中V是撒点区域的体积,这里就是8(2×2×2),于是上代码:
#include <stdio.h>
long double func(long double x,long double y,long double z)
{
long double t = x*y*z; //三重积分的方程
return t;
}
int main()
{
long double x,y,z; //坐标
unsigned long long i,j,k; //分割的算子
unsigned long long N; //分割成立方体的个数N*N*N
long double sum = 0;
long double answer;
printf("请给出对坐标的切割次数,输入应为正整数,越大越准\n");
scanf("%lld",&N);
for(i = 0;i<N;i++)
{
for(j = 0;j<N;j++)
{
for(k = 0;k<N;k++)
{
x = (2.0/((double)N))*i; //把x,y,z采样点的坐标算出来
y = (2.0/((double)N))*j;
z = (2.0/((double)N))*k;
if(x*x+y*y+z*z <= 4 && x*x+y*y+(z-2)*(z-2) <= 4) //如果在区域内,就加上
{
sum += func(x,y,z);
}
}
}
}
answer = (2.0*2.0*2.0)/((double)(N*N*N))*sum; //蒙特卡洛计算
printf("%lf\n",answer);
}
运行一下看结果:实际上答案为53/60,总值之如果愿意搞成很大的数字可以算得很准,不过要等很久