边界点a和b内的蒙特卡洛积分的Fortran代码

Question

我知道蒙特卡罗模拟是通过绘制随机点并计算曲线外部和曲线内部之间的比率来估算面积。

我已经很好地计算了pi的值，假设曲线半径为1。

这是代码

program pi
implicit none

integer :: count, n, i
real :: r, x, y
count = 0
n=500
CALL RANDOM_SEED
DO i = 1, n
 CALL RANDOM_NUMBER(x)
 CALL RANDOM_NUMBER(y)
 IF (x*x + y*Y <1.0) count = count + 1
END DO
r = 4 * REAL(count)/n
print *, r
end program pi

但为了找到整合，教科书说应用相同的想法。但是如果我想找到它的集成，我就迷失了如何编写代码

f(x)=sqrt(1+x**2) over a = 1 and b = 5

在半径为1之前，我确实假设点在x * 2 + y ** 2的条件下落入，但我怎样才能解决上面的问题？

任何帮助都非常有帮助

Answer 1

我将首先编写代码，然后解释：

Program integral
implicit none
real f
integer, parameter:: a=1, b=5, Nmc=10000000   !a the lower bound, b the upper bound, Nmc the size of the sampling (the higher, the more accurate the result)
real:: x, SUM=0

do i=1,Nmc                  !Starting MC sampling
  call RANDOM_NUMBER(x)     !generating random number x in range [0,1]
  x=a+x*(b-a)               !converting x to be in range [a,b]
  SUM=SUM+f(x)              !summing all values of f(x). EDIT: SUM is also an instrinsic function in Fortran so don't call your variable this, I named it so, to illustrate its purpose
enddo

print*, (b-a)*(SUM/Nmc)     !final result of your integral
end program integral

function f(x)           !defining your function
  implicit none
  real, intent(in):: x
  real:: f

  f=sqrt(1+x**2)
end function f

那么发生了什么：

积分可以写成。哪里：

（这个g（x）是[a，b]中变量x的均匀概率分布。我们可以把积分写成：

哪里。

所以，最后，我们得到的积分应该是：

因此，您所要做的就是生成[a，b]范围内的随机数，然后计算此x的函数值。然后这么做很多次（Nmc次），并计算总和。然后用Nmc除以找到平均值，然后乘以（b-a）。这就是代码的作用。

互联网上有很多东西。 here's一个可视化它非常好的例子

编辑：第二种方式，这与Pi方法相同：

Nin=0                    !Number of points inside the function (under the curve)
do i=1,Nmc
  call random_number(x)
  call random_number(y)
  x=a+x*(b-a)
  y=f_min+y(f_max-f_min)
  if (f(x)<y) Nin=Nin+1
enddo
print*, (f_max-f_min)*(b-a)*(real(Nin)/Nmc)

所有这些，你可以将它包含在一个外部do循环中，将（f_max-f_min）（b-a）（实数（Nin）/ Nmc）加起来，最后打印出它的平均值。对于这个例子，你所做的实际上是创建一个从a到b（x维度）和从f_min到f_max（y维度）的封闭框，然后对该区域内的点进行采样并计算函数中的点（ Nin）。显然，您必须知道[a，b]范围内函数的最小值（f_min）和最大值（f_max）。或者你可以为你的f_min f_max使用任意低/高值，但是你会浪费很多积分，你的错误会更大。