vegas

说明：这里只是对子程序的算法和如何分层做详细介绍，对打印输出和其他一些东西都做略过。还有这个子程序中涉及的变量很多，而且变量的意义也变化，所以下面涉及到的变量我都会在画流程图前加以适当的说明。

1．子程序所涉及的变量及其的意义：

有用输入参数：

FXN：被积函数。 ACC：所要求的精度。

NDIM：被积函数的变量的个数。 NCALL：投点的总数。 ITMX：要求的叠代的次数。

NPRN：打印的要求（具体要求略过）。 子程序内部的参数及其意义： TI（AVGI）：数值积分的结果。 IT：叠代的次数。

X（I，J）：各随机变量的分布，J表示各个变量。 ERR（SD）：数值积分的误差。 ALPH0：变化区间的参数。 RAND（20）：[0，1]的随机数列。

X（J）：用于积分的坐标等。

2．主要的程序模块和各模块的功能

总的流程图：

1．初始化设定各个参数 2．初始化内部的分布 3．初始化叠代值 4．计算各个随机点的位置X（J），计算各个随机点的函数值。 5．计算积分值和积分的误差 是 6．误差和叠代次数是否符合要求 否8．打印输出结果 7．按误差改变各区间的大小，重新获得一个变量的分布 结束

3．各模块的流程图

1）初始化设定各个参数

该部分的主要工作就是给各参数输入初始值，下面简述各个初始化的参数以及各参数的意义或在以下何模块里要用到。 ALPH=1.5; 这是在7中用到，作为变换分布的参数。 CALLS=NCALL； 投点的总数。 XND=ND=50； 区间分为50份。 NDM=ND-1=49； 为以下处理方便。 IT=0； 叠代初始为零。 SI=0； SI2=0； SWGT=0； SCHI=0；

SCALLS=0； 以上五个变量均是叠代参数，在5中用。 2）初始化内部分布：

变量的意义：

RC：分层的间隔。 NDIM：变量的个数。

XI（I，J）：变量的区间的分布。

XND： 区间分层数。 XNDM=XND-1。 DR： 求和号。

初始化的分布是一个均匀分布，实现的流程图如下：

RC=1/XND J=1 N J<=NDIM Y XI(ND,J)=1 DR=0 J++ I<=NDM Y N DR=DR+RC XI(I,J)=DR I++ CONTINUE

3) 初始叠代参数：

这不仅仅是初始化的过程了，而是大的循环的一部分，每改变一次变量的分布，就需要初始化叠代参数一次。下面简述各个初始化的参数以及各参数的意义或在以下何模块里要用到。

IT=IT+1； 叠代的次数。 TI=0； 积分值清零。

SFUN2=0； 在4，5中求误差时要用。

D（I，J）=0； 对D矩阵清零，在4，5，7中要用。

4）计算各个随机点的位置X（J），计算各个随机点的函数值。

这里主要是算各个随机点的位置X（J）和各个随机点的函数值，并为下面的求误差作好准备。总的流程图如下：

N JJ<=NCALL Y 计算各个随机点的位置X（J）CONTINUE 计算各个随机点的函数值。JJ++

下面分别看计算各个随机点的位置X（J）和计算各个随机点的函数值的流程图：

A）各个随机点的位置X（J） 变量定义：

NDIM： 变量的个数。 XI（I，J）： 变量的分布。 WGT： 权重。 流程图：

WGT=1 产生随机数列RAND（N） J<=NDIM CONTINUE XN=RAND（J）*XND+1 取XN的整数部分到数列IA（J）让XIM1=0 IA（J）>1 XIM1=XI（IA（J）-1，J）表示区间的距离X0=XI（IA（J），J）-XIM1。 X（J）=XIM1+（XN-IA（J））*X0。 权重WGT=WGT*X0*XND。 J++

这里是用随机数经过RAND（J）*XND+1的变换值的整数部分IA（J）选择了分层区间，用其小数部分来选择了坐标点在该分层区间的具体位置。该部分得到X（J），和该点的权重。 B）计算各个随机点的函数值

计算带权重的函数值： FUN=FXN（X）*WGT/CALLS 对下面求误差要用的值： FUN2=FUN*FUN SFUN2=SFUN2+FUN2 积分值： TI=TI+FUN J<=NDIM Y N 把D矩阵用FUN2填充： J++，IAJ=IA（J） D（IAJ，J）= D（IAJ，J）+FUN2CONTINUE

这里求出了带权重的函数值FUN，以及下面积分和求误差要用的几个量：SFUN2，FUN2，D阵，TI。其中D阵对第7个模块也很重要。 D矩阵的矩阵元是对第J个变量的第I个区间内所有投点的FUN2的求和。

5）计算积分值和积分的误差

该部分是计算了在以上的投点和分布下的积分值和积分的误差值。注意这个积分值和积分误差是所有叠代的平均值。 具体的流程图如下：

在4中求得的本次分布下的积分值的平方 TI2=TI*TI； 本次积分值TI的误差： TSI=(SFUN2*CALLS−TI2)/(CALLS−1) 上式相当于： CALLSTSI=(∑i=1FUN−2CALLS∑i=1FUN2)/[(CALLS−1)CALLS]权重： WGT=TI2/TSI2以前所有的叠代值之和 SI=SI+TI*WGT SI2=SI2+TI2 SWGT=SWGT+WGT SCHI=SCHI+TI2*WGT 所有叠代平均积分值：AVGI=SI/SWGT；所有叠代误差的平均： SD=SWGT*IT/SI2=∑1/TSI2 ERR=SD*100/AVGI

该部分计算出了所有叠代的平均积分值AVGI，所有叠代的相对误差平均ERR，在本次分布下的积分TI，本次积分的误差TSI 。 6）判断误差是否符合要求

该部分只是一个判断语句：

SD

(误差)>ACC 且ITAVGI

[其中第一个是表示精度是否符合要求，第二个是表示叠代次数是否满足要求]

如果成立就进入7，否则就可以按照要求打印输出了，即进入8。 7）按误差改变各区间的大小，重新获得一个变量的分布

这部分是VEGAS的精华所在，其中的有些算式其意义也不是很明确，只能作为定义给出。该部分的主要工作就是改变X（I，J）的分布，对每个变量的分层区间作了一个调整，其调整的具体依据是前面4中得到的D（IAJ，J）矩阵。具体操作先用下面的主流程图表示：

对D矩阵进行平滑操作，得到新的D矩阵和新D矩阵各列的和X（J）N J<=NDIM Y 根据前面的D矩阵和X（J）来构造R（I）和R（I）的和J变量的RC。回到模块3根据RC和R（I）来重新给出J变量的分层区间X（I，J）。J++

下面来仔细看以上三部分的具体流程图：

A） 对D矩阵进行平滑操作，得到新的D矩阵和新D矩阵各列的

和X（J） 流程图：

N J<=NDIM Y给变量初值： X0=D（1，J）；XN=D（2，J）；CONTINUE 给D矩阵和X（J）初值：D（1，J）=（X0+XN）/2 X（J）=D（1，J） N I<=NDIM Y对D矩阵调整：（该循环中I初值为2） D（I，J）=X0+XN，X0=XN，XN=D（I+1，J）D（I，J）=（D（I，J）+XN）/3 效果就是新D阵中D（I，J）是原来D阵中 D（I-1，J），D（I，J），D（I+1，J）和的平均。求X（J）：X（J）=X（J）+D（I，J） I++求最后一项的D（ND，J）和X（J）终值 D（ND，J）=（XN+X0）/2 X（J）=X（J）+D（ND，J） J++

B）根据前面的D矩阵和X（J）来构造R（I）和R（I）的和J变量的RC。 流程图如下：

初值设定：RC=0 I<=ND YN R(I)清零R（I）=0CONTINUE Y D（I，J）<=0 N 用下面的方法求R（I）： X0=X（J）/D（I，J） R(I)=(X0−1ALPH)X0*lnX0 求RC，RC是对R（I）的求和：RC=RC+R（I） I++

这里ALPH是定义的一个参数，上述的公式不清楚如何得来。 C）根据RC和R（I）来重新给出J变量的积分区间X（I，J）。 这部分重新定义了每个积分区间的长度，从而给出了新的X（I，J），

要用到前面的值RC，R（I），XND（是总的分的积分区间的个数）。流程图如下：

让RC均分：RC=RC/XND给定初值： K=0，XN=0，DR=0，I=0 开始寻找区间： K=K+1；DR=DR+R（K）X0=XN，XN=XI（K，J）Y RC>DR N 区间的标号：I=I+1 ，DR=DR-RC， 区间的断点： XIN（I）=XN-（XN-X0）*DR/R（K）NI<=NDM Y记入最终的分布X（I，J）：XI（I，J）=XIN（I） I++ X（I，J）最后一行的值（类比于模块2的X（I，J）的最后一行值）： X（ND，J）=1 CONTINU

控制新的分布的量是RC和R（I），刚开始时由B）可知RC是所有R（I）的和，然后把RC取为其原来的XND分之一，再由R（I）的累加得到DR，当DR超过RC时，取R（I）的标记I值取原来分布的X（I，J）值来得到XN值，而X（I-1，J）作为X0值，最后得到新。这就是变换得到新分布的XI（I，J）==XN-（XN-X0）*DR/R（K）的过程。 8）打印输出结果

这一部分就是按要求即输入参数NPRN[打印的要求（具体要求略过）]来打印输出，这里不作详述了。

总上前面的内容，已经基本把VEGAS程序的思想说清楚了。

吴佳俊 PB03203076