问题描述
我正在尝试评估这个 5 倍积分
其中 a
、w
、t
、x1
、x2
是常量,y1(x)
、y2(x)
是非依赖于 y
的 x
的恒定积分边界。我可以评估它们,然后使用 nquad
的 scipy.integrate
函数计算积分,如下所示:
from scipy.integrate import nquad
from math import sqrt
# UNIT
a = 2.0 # unitless
w = 100e-9 # meter
t = 20e-9 # meter
delta = 20e-9 # meter (used to calculate x1,x2,y1,y2)
def integrand(y_prime,z_prime,y,x,z):
num = (-1/2)*((x-a*w/2) + sqrt(3)*(y-y_prime))
den = ((x-a*w/2)**2 + (y-y_prime)**2 + (z-z_prime)**2)**(3/2)
return num/den
# integration bounds of x
x1 = (1/2)*(a*w + 3*delta)
x2 = (1/2)*(w*(a + sqrt(3)) + 3*delta)
# integration bounds of y
def y1(x):
m1 = -1/sqrt(3)
n1 = (1/2)*(w*(a/sqrt(3) + 1) + delta*(1 + sqrt(3)))
return m1*x + n1
def y2(x):
m2 = sqrt(3)
n2 = (1/2)*(w*(1 - a*sqrt(3)) + delta*(1 - 3*sqrt(3)))
return m2*x + n2
integral = nquad(integrand,[[-w/2,w/2],# y' bounds
[0,t],# z' bounds
lambda x,z: [y1(x),y2(x)],# y bounds
[x1,x2],# x bounds
[0,t]]) # z bounds
print(integral)
输出为 (-2.3258337548893193e-23,0.42437932862557626)
。
我有充分的理由相信值 -2.3258337548893193e-23
是正确的(至少,数量级符合预期)。但是,我担心绝对误差估计的 0.42437932862557626
的极大值,因为它不仅比积分大,而且超过了积分 22 个数量级!
我尝试使用 epsabs
和 epsrel
可选参数的几种方法来改进错误估计(例如,按照 this article 的建议),但它们都没有正常工作:我能得到的最好的结果是 2.7568352154734943e-06
的误差估计值与积分值相同。
我的最后一次尝试是更改常量参数的值,这恰好是物理长度。因此,我没有以米为单位设置它们的值,而是将它们设置为 nano 米:
# UNIT
a = 2.0 # unitless
w = 100 # nanometer
t = 20 # nanometer
delta = 20 # nanometer (used to calculate x1,y2)
检查积分,此更改应仅按 (1e9)**3=1e27
的系数缩放其值(然后,如果我想要其值以米表示的长度,我只需将其除以 1e27
)。而且,实际上,输出是 (-23258.337548893207,1.0567318670726077e-06)
,积分值按预期重新调整。更妙的是,该误差并没有随着时间的推移而扩大,并且现在具有可接受的值!
虽然我显然绕过了我的问题,但我仍然想知道我是否可以使用以米为单位表示的参数获得好的结果。或者,至少理解为什么在第一种情况下错误如此之大,而在第二种情况下却没有。看来得跟积分的特定值有关:“大”时误差较小,“小”时误差较大。
解决方法
暂无找到可以解决该程序问题的有效方法,小编努力寻找整理中!
如果你已经找到好的解决方法,欢迎将解决方案带上本链接一起发送给小编。
小编邮箱:dio#foxmail.com (将#修改为@)