蒙特卡罗方法在积分计算中的应用

例：计算定积分 S ( x sin x )dx
2 2.5
8.4
主程序： EXTERNAL F DOUBLE PRECISION F，S，A，B DATA A，B/2.5D0，8.4D0/ CALL FMTCL（A，B，F，S） WRITE（*，*） WRITE （*，10）S 10 FORMAT（1X，‘S=’，D13.6） WRITE（*，*） END
子程序FMTCL.FOR
SUBROUTINE FMTCL(A,B,N,F,X,S) DOUBLE PRECISION A(N),B(N),F,S,R,X(N),K REAL M,NRND1 F=1.0D0 M=10000.0 K=10000.0D0 S=0.0D0 IF (M+1.0.NE.1.0)THEN M=M-1.0 DO 20 I=1,N X=A+(B-A)*NRND1(R) CONTINUE S=S+F(N,X)/K GOTO 10 ENDIF DO 30 I=1,N S=S*(B(I)-A(I)) END
形参说明
N：整型变量，输入参数，积分的重数。 A,B： 均为双精度实型一维数组，长度为N，输入参数， 积分的下限值和上限值。 F：双精度实型函数子程序名,输入参数。用于计算被积函 数值f(x1,x2,…,xn)。在主程序中必须用外部语句及类型 说明语句对相应遥实参进行说明。 DOUBLE PRECISION FUNCTION F（N,X） 其中：X为双精度实型一维数组，长度为N ；用于存放N 个自变量值。 S：双精度实型变量，输出参数。返回积分值。
例：计算定积分
S
2 2 2
1 1 1

2 2 ( x12 x2 x3 )dx1dx2dx3
主程序： EXTERNAL F DIMENSION X(3)，A(3)，B(3) DOUBLE PRECISION X，F ， S，A，B DATA A，B/3*1.0D0，3*2.0D0/ N=3 CALL FMTCL（A,B,N,F,X,S） WRITE（*，*） WRITE （*，10）S 10 FORMAT（1X，‘S=’，D13.6） WRITE（*，*） END
蒙特卡罗方法在积分计算中的应用
计算多重积分是蒙特卡罗方法的重要应用 领域之一。本章着重介绍计算定积分的蒙特卡 罗方法的各种基本技巧，而这些技巧在粒子输 运问题中也是适用的。
1. 蒙特卡罗方法求积分
蒙特卡罗方法求积分的一般规则如下：任何一个 积分，都可看作某个随机变量的期望值，因此，可以 用这个随机变量的平均值来近似它。 蒙特卡洛方法的基本思想： 当问题可以抽象为某个确定的数学问题时，应当首先 建立一个恰当的概率模型，即确定某个随机事件A或随机 变量X，使得待求的解等于随机事件出现的概率或随机变 量的数学期望值。然后进行模拟实验，即重复多次地模拟 随机事件A或随机变量X。最后对随机实验结果进行统计平 均，求出A出现的频数或X的平均值作为问题的近似解。这 种方法也叫做间接蒙特卡洛模拟。
计算一维积分的蒙特卡洛方法
说明:设定积分为
S f ( x)dx
a
b
取0到1之间均匀分布的随机数序列 ξ i(i=1,2,…,m),并令
xi a (b a)i
只要m足够大,则有
ba m S f ( xi ) m i 1
形参说明
A,B： 均为双精度实型变量，输入参数，积分的下限和 上限。 F,双精度实型函数程序名,输入参数。用于计算被积函数 值f(x)。在主程序中必须用外部语句及类型说明语句对 相应遥实参进行说明。 DOUBLE PRECISION FUNCTION F（X） 其中：X为双精度实型变量，自变量值；函数名F返回双 精度实型被积函数值。 S：双精度实型变量，输出参数。返回积分值。
10
20
30
被积函数子程序
DOUBLE PRECISION FUNCTION F(N,X) DOUBLE PRECISION X(N),F F=0.0D0 DO 10 I=1,N F=F+X(I)*X(I) END
10
计算结果:S=.697043D+01
多重积分的计算, Monte Carlo 方法:
由上述第二种表述, 物理量的平均值的计算归结为一个多重积分的计算, 设系统 由100个粒子构成, 每个粒子有6个自由度, 所以需要计算600重积分. 现在考虑在 每一维取10个点, 总共有10600个点. 假设计算机每秒可以计算1012个点, 计算 这个积分需要10588秒!!!!! 问题比这个更严重! 如果如上述方式取点, 则积分区域的内点数为 8600, 在总的点 中所占比例为 (8/10)600=0.714£ 10-58, 也就是说, 取的点基本上都在表面上!!!
子程序FMTCL.FOR
SUBROUTINE FMTCL(A,B,F,S) DOUBLE PRECISION A,B,F,S,R,X,K REAL M,NRND1 F=1.0 M=10000.0 K=10000.0D0 S=0.0D0 IF (M+1.0.NE.1.0)THEN M=M-1.0 X=A+(B-A)*NRND1(R) S=S+F(X)/K GOTO 10 ENDIF S=S*(B-A) END
取0到1之间均匀分布的随机数点列 i i i (1 ,2 ,n ) (i=1,2,…,m),并令
x a j (bj a j )
(i ) j
(i ) j
(j=1,2,…,n)
只要m足够大,则有
n 1m (i ) (i ) S f ( x1( i ) , x2 xn ) * (b j a j ) m i 1 j 1
10
被积函数子程序
DOUBLE PRECISION FUNCTION F(X) DOUBLE PRECISIONX F=X*X+SIN(X) END 计算结果:S=.191553D+03
计算多维积分的蒙特卡洛方法
说明:设定积分为
S
b1 b2
a1 a2

bn
anwenku.baidu.com
f ( x1, x2 xn )dxdxdx