Я пытаюсь выполнить интеграцию для задачи физики, и код, который я написал, дает мне результаты - в 10 раз больше. Я задавался вопросом, может ли кто-то указать мне правильное направление относительно того, является ли это моя отвратительная тройка для l oop или что-то еще, что идет не так.
Я пытаюсь сделать этот расчет. ( Это из этой статьи о расчете энергии основного состояния атома лития с использованием координат Hylleraas, если вам интересно !! )
Это соответствующая часть статьи, и я объясню, как я разбил его ниже.
Короче говоря, чтобы получить значение для интеграла, бесконечного сумма (5) усекается после 10 значений q. Формула 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
Оно основано на продукте с пи-полными числами, биномиальном коэффициенте и продукте, но довольно просто, и я не могу определить любые ошибки.
Также полагается на массу 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) при конец. Когда я изучал вывод, неправильное значение кажется самым первым.
Это потому, что я неправильно понимаю свои диапазоны?
Я понимаю, что это довольно сложный вопрос, но я бы больше всего благодарен за любую помощь!