通俗易懂的Monte Carlo的积分方法(三)

通俗易懂的Monte Carlo的积分方法(三)

• 考虑曾经在参加MCM时的一个多重积分的计算难题
  $$\int {\int }_D(H-z(x,y){)}^2\sqrt{1+z_x^2+z_y^2}d\sigma$$
  其中：$\frac{x^2}{a^2}+\frac{y^2}{b^2}+\frac{z^2}{c^2}=1$且$z\ge 0$。
  

• 采用Monte Carlo积分法(期望法)

  用的是面向对象的模块化编程思想

• Python代码：

  import random
  import numpy as np
  import math
  import matplotlib.pyplot as plt
  def z_value(a,b,c,x,y):
      return c*math.sqrt(1-x**2/(a**2)-y**2/(b**2))
  
  def z_dotx(a,b,c,x,y):
      return -c**2/a**2*x/z_value(a,b,c,x,y)
  
  def z_doty(a,b,c,x,y):
      return -c**2/b**2*y/z_value(a,b,c,x,y)
  
  def function_Calcul(x,y,H,a,b,c):
      if((x**2/(a**2)+y**2/(b**2))<1):
          return ((H-z_value(a,b,c,x,y))**2)*math.sqrt(1+(z_dotx(a,b,c,x,y))**2+z_doty(a,b,c,x,y)**2)
      else:
          return 0
  
  def int_Calcul(a,b,c,H,numbers):
      sum = 0
      Prop = []
      Prop_final = []
      S = 2 * a * 2 * b
      for i in range(1,numbers):
          x_value = -a+2*a*random.random()
          y_value = -b+2*b*random.random()
          sum += function_Calcul(x_value,y_value,c,a,b,c)
          Prop.append(sum/(i)*S)
      for i in range(1,numbers):
          Prop_final.append(sum/(numbers)*S)
      plt.xlabel("Numbers")
      plt.ylabel("IntergalValue")
      plt.plot(range(1,numbers),Prop,color = 'blue',label = 'Change_value' )
      plt.plot(range(1,numbers),Prop_final,color = 'red',label = 'Mean_value')
      plt.legend()
      plt.show()
      print("The final intergal value is:{:.2f}".format(Prop[numbers-2]))
  
  int_Calcul(0.5,1,2,2,10000)
  

• 计算结果：
  
  当$a=0.5,b=1,c=2,H=2$时，其最终的结果大概在$12.05$左右。