三重嵌套和 Numpy

Question

我正在尝试对一个物理问题进行积分，而我编写的代码给我的结果大约大了 10 倍。我想知道是否有人可以为我指出正确的方向，到底是我可怕的三重 for 循环还是其他问题。

我正在尝试进行此计算。 (It's from this paper about calculating the ground state energy of a Lithium atom using Hylleraas co-ordinates if you're interested !!)

这是论文的相关部分，我将在下面解释我是如何分解它的。

长话短说，为了得到积分的值，无穷和 (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_integral。我相信它会为输入的任何值计算正确的值：我不太确定正确的值会进入它（因此打印语句）！

这是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)。当我一直在调查输出时，第一个值似乎不正确。

这是因为我的范围有误吗？

我知道这是一个非常重要的问题，但如果能提供任何帮助，我将不胜感激！

Answer 1

您是否尝试过使用一些数值方法，例如欧拉的前向和后向方法？

欧拉倒推法推导：

让我们考虑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 倍的导数。