Как мне численно интегрировать функцию в Python? - PullRequest
Новая
       36

Как мне численно интегрировать функцию в Python?

3 голосов
/ 24 мая 2019

Я пытаюсь написать код Python (на основе записной книжки Mathematica), который включает в себя некоторые интеграции, но мне пока не повезло.
Интеграл, который я пытаюсь оценить:

Integral

, где

Ls matrices

* Это для Bjorken-Mtingwa Intra-расчеты рассеяния луча.

В Python я пытаюсь: 

import numpy as np
import math
from sympy import Symbol
from scipy import interpolate
import scipy.integrate as integrate
from scipy.integrate import quad

def Aint(Np, r, c , beta_rel, gamma_rel, emit_x, emit_y, Sigma_s, Sigma_M):
    return (Np * r**2 *c) / ( 64 * np.pi**2 * beta_rel**3 * gamma_rel**4 * emit_x * emit_y * Sigma_s * Sigma_M)

def MatrixSum(M1 ,M2 ,M3):
    return [[M1[i,j] + M2[i,j] + M3[i,j] for j in range (M1.shape[0])] for i in range(M1.shape[1])]

#Then I'm initializing all my parameters by loading a file into a dataframe and doing some calculations. I will not include most of that part to keep it short. I have no complex numbers.
gamma_rel = 1.00451006711
beta_rel  = 0.09465451392
beta_x = 7.890105185
beta_y = 13.61578059
Phi_x = -1.957881913
Phi_y = 0.0
emx = 2.95814809e-06
emy = 2.95814809e-06
H_x = 32.68714662287
H_y = 0.0
bl = 4.256474951
r0 = 2.16775224067e-17
Sigma_M = 0.00118124786
II = np.identity(3)

ABM = Aint(Np, r0, c, beta_rel, gamma_rel, emx, emy, bl, Sigma_M)
Clog = 13.53496204 #Evaluating the Coulomb logarithm with a function but seems correct so I will not include the calculations.

Lp = (gamma_rel**2 / Sigma_M) * np.matrix( [ [0,0,0], [0,1,0], [0,0,0] ] )
Lx = (beta_x / emx) * np.matrix( [ [1, - gamma_rel * Phi_x, 0], [- gamma_rel * Phi_x, gamma_rel**2 * H_x / beta_x , 0], [0, 0, 0] ] )
Ly = (beta_y / emy) * np.matrix( [ [0, 0, 0], [0, gamma_rel**2 * H_y / beta_y, - gamma_rel * Phi_y], [0, - gamma_rel * Phi_y, 1] ] )
L  = np.matrix(MatrixSum(Lx, Ly, Lp))

Ix = integrate.quad(lambda x: (x**(1/2.0) * (np.trace(Lx) * np.trace(np.linalg.inv(L + x * II)) - 3 * np.trace(np.matmul(Lx, np.linalg.inv(L + x * II)))) / np.linalg.det(L + x * II)**(1/2.0)) , 0, np.inf)
Ixx = 4 * np.pi * ABP * Clog * Ix[0]

#Similarly for the other 2 integrals. In reality, all 3 integrals are evaluated in a double loop.

, но я получаю разные результаты от mathematica.Я также попробовал scipy.integrate.simps, но это тоже не помогло.

В Mathematica я просто интегрирую его с: 

Ix = NIntegrate[Intx, {x, 0, inf}, MaxRecursion -> 100];

с Intx, являющимся неотъемлемой частью фотографии, и той же процедурой, что и ранее.

Есть ли рекомендации по эффективной интеграции этой функции?Что-то не так с моим методом?

...