在Python中数值求解二重积分

Question

我正在尝试使用矩形面积来数值求解二重积分，其中积分看起来只是两个变量（theta 和 phi）的函数。

这是我迄今为止的尝试：

from math import*

import numpy as np

N = 100000

T = 100000

j = 1

m = 1

lphi = 0.0

uphi = pi

ltheta = pi/2.0

utheta = 0.0

sum2 = 0.0

def fa(theta,phi):
  fac = sqrt(1.0 - ((sin(theta))**2)*((cos(phi))**2))
  top = -np.sign(cos(phi))*fabs(sin(i)*cos(phi))
  fangular = (1.0/fac)*((pi/2.0) - atan(top/fac)) - 1.0
  return fangular


while m < T:
  thetaa = ltheta + m
  dphi = fabs((uphi - lphi)/T)
  dtheta = fabs((utheta - ltheta)/N)
  sum1 = 0.0
  while j < N:
      phia = lphi + j*dphi
      fun = fa(thetaa,phia)
      Area1 = fun*dphi
      sum1 += Area1
      j +=1
  Area2 = sum1*di
  sum2 += Area2
m+=1
print sum2

我只是有点担心第二个积分（涉及 theta）没有更新。

Answer 1

在 phi 循环开始之前，内积分 (sum1) 应重置为 0。目前，它将在 theta 迭代中不断累积值。

m = 1
lphi = 0.0
uphi = pi
ltheta = pi/2.0
utheta = 0.0 
sum2 = 0.0

def fa(theta,phi):
   fac = sqrt(1.0 - ((sin(theta))**2)*((cos(phi))**2))
   top = -np.sign(cos(phi))*fabs(sin(i)*cos(phi))
   fangular = (1.0/fac)*((pi/2.0) - atan(top/fac)) - 1.0
   return fangular

while m < T:
   thetaa = ltheta + m
   dphi = fabs((uphi - lphi)/T)
   dtheta = fabs((utheta - ltheta)/N)
   sum1 = 0.0
   while j < N:
      phia = lphi + j*dphi
      fun = fa(thetaa,phia)
      Area1 = fun*dphi
      sum1 += Area1
      j +=1

 Area2 = sum1*di
 sum2 += Area2
 m+=1
 print(sum2)