python - 三重嵌套和 Numpy
问题描述
我正在尝试对一个物理问题进行集成,而我编写的代码给我的结果大约是 10 倍大。我想知道是否有人可以为我指出正确的方向,即这是我可怕的三重 for 循环还是其他问题。
我正在尝试进行此计算。(如果你有兴趣,它来自这篇关于使用 Hylleraas 坐标计算锂原子的基态能量的论文!!)
这是本文的相关部分,我将在下面解释我是如何分解它的。
长话短说,为了得到积分的值,无穷和 (5) 在 q 的 10 个值之后被截断。公式 T (q) 在 (6) 中给出,它本身是三个值的三重嵌套和:k_1_2、k_2_3 和 k_3_1。
这是我的 T(q) 代码:
def T_q(j1, j2, j3, j_1_2, j_2_3, j_3_1, alpha, beta, gamma,q):
'''
T_q formula for I integral summation
'''
#print("q", q )
L_1_2 = 1/2 * (j_1_2 +1) #sets adjusted values of j12 etc.
#print(L_1_2, "L12")
L_2_3 = 1/2 * (j_2_3 +1)
#print(L_2_3, "l23")
L_3_1 = 1/2 * (j_3_1 +1)
#print(L_3_1, "L31")
j_1 = j1 +2
j_2 = j2 +2
j_3 = j3 +2
t_q = 0
for k_1_2 in np.arange(L_1_2 + 1): # Triple for loop for the triple sum
#print("k_1_2", k_1_2)
for k_2_3 in np.arange(L_2_3 + 1):
# print("k_2_3", k_2_3)
for k_3_1 in np.arange(L_3_1 + 1):
# print("new loop")
#print("k_3_1", k_3_1)
W_mess = (W_integral((j_1 + 2*q + 2*k_1_2 + 2*k_3_1), (j_2 + j_1_2 - 2*k_1_2 + 2*k_2_3), (j_3 + j_2_3 -2*q -2*k_2_3 + j_3_1 - 2*k_3_1),alpha, beta, gamma) +
W_integral((j_1 + 2*q + 2*k_1_2 + 2*k_3_1), (j_3 + j_3_1 - 2*k_3_1 + 2*k_2_3), (j_2 + j_1_2 -2*q -2*k_1_2 + j_2_3 - 2*k_2_3),alpha, gamma, beta) +
W_integral((j_2 + 2*q + 2*k_1_2 + 2*k_2_3), (j_1 + j_1_2 - 2*k_1_2 + 2*k_3_1), (j_3 + j_2_3 -2*q -2*k_2_3 + j_3_1 - 2*k_3_1),beta, alpha, gamma) +
W_integral((j_2 + 2*q + 2*k_1_2 + 2*k_2_3), (j_3 + j_2_3 - 2*k_2_3 + 2*k_3_1), (j_1 + j_1_2 -2*q -2*k_1_2 + j_3_1 - 2*k_3_1),beta, gamma, alpha) +
W_integral((j_3 + 2*q + 2*k_2_3 + 2*k_3_1), (j_1 + j_3_1 - 2*k_3_1 + 2*k_1_2), (j_2 + j_1_2 -2*q -2*k_1_2 + j_2_3 - 2*k_2_3),gamma, alpha, beta) +
W_integral((j_3 + 2*q + 2*k_2_3 + 2*k_3_1), (j_2 + j_2_3 - 2*k_2_3 + 2*k_1_2), (j_1 + j_1_2 -2*q -2*k_1_2 + j_3_1 - 2*k_3_1),gamma, beta, alpha))
t_q += (1/((2*q+1)**2)) * C_constant(j_1_2,q,k_1_2) * C_constant(j_2_3,q,k_2_3) * C_constant(j_3_1,q,k_3_1) * W_mess
#print("t_q, ",t_q)
#print("t_q final",t_q)
return t_q
(请原谅打印功能,我正在使用这些功能来确保每次迭代的正确值被触发 - 它们是我所看到的)
其中每一个都有一个常数值,其公式由(4)给出,我使用这个 python 函数计算:
def C_constant(j,q,k):
'''
Calculates C constant
'''
S_q_j = np.minimum( (q-1), (j+1)/2 ) # takes minimum
constant_term = (2*q+1)/(j+2)
binomial_term = sc.binom(j+2,(2*k+1))
product = mp.nprod(lambda t: ((2*k + 2*t -j )/(2*k +2*q - 2*t +1)), [0,S_q_j] )
numpy_product = np.double(product)
C = constant_term * binomial_term * numpy_product
return C
它依赖于大写 pi 乘积、二项式系数和乘积,但相当简单,我无法发现任何错误。
它还依赖于加在一起的大量 W_integrals。我相信它会为输入的任何值计算正确的值:我不太确定正确的值会进入它(因此打印语句)!
这是W代码
def W_integral(l,m,n,alpha,beta,gamma):
'''
W integral taken from this paper https://journals.aps.org/pra/abstract/10.1103/PhysRevA.52.3681
Asks for l m n values + alpha beta gamma and returns equation (7), in said paper
Checked against Matlab code
'''
constant = np.math.factorial(l)/((alpha +beta +gamma)**(l+m+n+3))
W_sum = mp.nsum(lambda p: ((np.math.factorial(l+m+n+p+2))/((l+m+2+p)*np.math.factorial(l+1+p)) * ((alpha/(alpha +beta +gamma))**p)) * constant * mp.hyp2f1(1,l+m+n+p+3,l+m+p+3,(alpha+beta)/(alpha +beta +gamma)) ,[0,mp.inf])
numpy_W= np.double(W_sum)
return numpy_W
然后对每个 T(q) 值求和,以在此函数中给出最终结果:
def I_integral(j1, j2, j3, j_1_2, j_2_3, j_3_1, alpha, beta, gamma):
'''
Takes values for power of electron co-ordinates and returns the value of I "
'''
N = 10
I_0_N = ((4*np.pi)**3) * np.array([T_q(j1, j2, j3, j_1_2, j_2_3, j_3_1, alpha, beta, gamma,q) for q in np.arange(N+1)]).sum()
#print("I_0_N before constant")
return I_0_N
问题是,目前,与该表相比,我的值I_integral(0,0,0,-1,-1,1,1,1,1)
大约是该表给出的值的 10 倍:
完成的 T(q) 总和最后乘以 (64*pi^3) (~2000)。当我一直在调查输出时,第一个值似乎不正确。
这是因为我的范围错了吗?
我意识到这是一个非常重要的问题,但我将不胜感激任何帮助!
解决方案
您是否尝试过使用一些数值方法,例如欧拉的前向和后向方法?
用欧拉倒推法推导:
让我们考虑一下Ts
我们的采样时间。那么,导数的近似值是:
dX(t) / dt = [x(t) - x(t-Ts)] / Ts
如果我们将连续时间平面映射到离散时间 z 平面(即s = e (s * Ts)),我们得到:
dX[k] = [x[k] - x[k-1]] / Ts
,其中 k 是离散时间瞬间。
让我们以以下信号为例:
X = np.linspace(0,100,100000)
y = np.sin(X*0.1)
然后,在 python 中,我们可以构建一个这样的函数:
def euler_backard_method(X, Ts):
"""
Computes the Euler's Backwards Method Numerical Derivative.
Arguments:
X: an input array
Ts: the sampling time
Output:
the derivative of X
"""
return [(X[idx]-X[idx-1])/Ts for idx in range(1,len(X))]
对于 Ts = 0.001(高频),我们得到 X ( euler_backard_method(X=y, Ts=1)
) 的以下输出:
导数示例
我们可以以同样的方式构建集成!
一体化
- 转发方法:
s <- (z-1) / Ts
- 后退法:
s <- (z-1) / Ts * z
- 梯形法:
s <- 2 * (z-1) / Ts * (z+1)
后向方法将变为:u[k] = u[k-1] + Ts * x[k]
,其中u[k]
是 的积分输出X
。相应的功能将是:
def backward_integration(X, Ts):
"""
Numerical Integration using backward method.
Arguments:
X: the input data
Ts: the sampling period
Output:
The derivative of X
"""
U = []
u = 0
for idx in range(len(X)):
u+=Ts*X[idx]
U.append(u)
return U
注意:结果受采样时间的影响很大,正如您所指出的,您可能会得到大约 10 倍的导数。
推荐阅读
- python - Tiled / pytmx - 在不同地图中为相同图像获取相同 gid 的问题
- hyperlink - 从 youtube 页面的侧边栏中提取 youtube 链接
- python - pipenv CLI 可以使用 --keep-outdated 选项安装到开发包中吗?
- sql - 在 SQL Server 中将 15 分钟转换为 10 分钟时间序列
- swift - Mac OS 桌面音频主机应用程序无法加载某些 AVAudioUnit 插件
- git - 有没有办法指定一个自定义工具在进行合并时始终由 git 使用?
- python-3.x - 无论如何要格式化由 map 函数制成的字符串?
- java - Java 中的 If Else
- python - 关于从标准输入随机化的 Python 问题
- c# - 使用 EF 防止公开更改收集内容