Triple Nested Sum Numpy

Question

I'm trying to perform an integration for a physics problem and the code I've written is giving me results ~10 times too big. I wondered whether someone could point me in the right direction as to whether it is my hideous triple for loop or something else that is going wrong.

I'm trying to do this calculation. (It's from this paper about calculating the ground state energy of a Lithium atom using Hylleraas co-ordinates if you're interested !!)

This is the relevant part of the paper and I will explain how I have broken it down below.

Long story short, to get the value for the integral, the infinite sum (5) is truncated after 10 values of q. The formula T (q) is given in (6), and is a triple nested sum itself over three values: k_1_2, k_2_3 & k_3_1.

This is my code for 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

(Please excuse the print functions, I was using these to try to ensure the correct values of each iteration were firing - they were as far as I could see)

Each of these has a constant value whose formula is given by (4) which I calculate using this python function:


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

It relies on a capital pi product, binomial coefficient and product but is rather simple and I can't spot any errors.

It also relies on a mass of W_integrals added together. I'm confident that it calculates the correct value for whatever values entered into it: I'm less certain the right values are going into it (hence the print statements)!

This is the W code

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

Each T(q) value is then summed to give the final result in this function:

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

The problem is that currently, when compared to this table, my value of I_integral(0,0,0,-1,-1,1,1,1,1) is about 10 times off the values given by this table :

The completed T(q) sum is times by (64*pi^3) (~2000) at the end. When I've been investigating the output, the incorrect value seems to the very first one.

Is this because I am getting my ranges wrong?

I realise this is quite a meaty question but I would be most grateful for any help!

user_12713451 · Answer

Have you tried using some numerical approach, such as Euler's Forward and Backward Methods?

Derivation with Euler's Backwards Method:

Let's consider Ts our sampling time. Then, an approximation for the derivative is:

dX(t) / dt = [x(t) - x(t-Ts)] / Ts

If we map the continuous time plane into the discrete-time z-plane (i.e. s = e^{(s * Ts)}), we obtain:

dX[k] = [x[k] - x[k-1]] / Ts, where k is a discrete-time instant.

Let's consider the following signal as an example:

X = np.linspace(0,100,100000)
y = np.sin(X*0.1)

Then, in python, we could build a function like this:

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))]

For Ts = 0.001 (high frequency), we get the following output for X (euler_backard_method(X=y, Ts=1)):

Derivative Example

We could build an integration in the same fashion!

Integration

forward method: s <- (z-1) / Ts
backward method: s <- (z-1) / Ts * z
trapezoidal method: s <- 2 * (z-1) / Ts * (z+1)

The backward method would become: u[k] = u[k-1] + Ts * x[k], where u[k] is the integrated output of X. The corresponding function would be:

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

Notice: the results are highly influence by the sampling time, where you could get a derivative ~10 times too big, as you pointed out.

Triple Nested Sum Numpy

1 Answers