Расщепление интегрированной плотности вероятности в двух пространственных областей

Question

У меня есть функция плотности вероятности:

T = 10000 
tmin = 0 
tmax = 10**20 
t = np.linspace(tmin, tmax, T) 
time = np.asarray(t)     #this line may be redundant 
for j in range(T): 
    timedep_PD[j]= probdensity_func(x,time[j],initial_state) 

Я хочу, чтобы интегрировать его в течение двух различных областей х. Я попытался следующим разделить timedep_PD массив в двух пространственных областях, а затем приступил к интеграции:

step = abs(xmin - xmax)/T 
l1 = int(np.floor((abs(ab - xmin)* T)/abs(xmin - xmax))) 
l2 = int(np.floor((abs(bd - ab)* T)/abs(xmin - xmax))) 

#For spatial region 1 
R1 = np.empty([l1]) 
R1 = x[:l1] 
for i in range(T): 
    Pd1[i] = Pd[i][:l1] 

#For spatial region 2 
Pd2 = np.empty([T,l2]) 
R2 = np.empty([l2]) 
R2 = x[l1:l1+l2] 
for i in range(T): 
    Pd2[i] = Pd[i][l1:l1+l2] 

#Integrating over each spatial region 
for i in range(T): 
    P[0][i] = np.trapz(Pd1[i],R1) 
    P[1][i] = np.trapz(Pd2[i],R2) 

Есть ли проще/более понятный способ идти о раскалывается функция плотности вероятности в двух пространственных областях, а затем интегрируя в каждой пространственной области на каждом временном шаге?

Answer 1

Петли могут быть исключены с помощью векторизованных операций. Неясно, является ли Pd двумерным массивом NumPy; это что-то еще (например, список списков), он должен быть преобразован в массив 2D NumPy с np.array(...). После этого вы можете сделать это:

Pd1 = Pd[:, :l1] 
Pd2 = Pd[:, l1:l1+l2] 

Не нужно зацикливаться на индексе времени; нарезка происходит на все времена сразу (с : вместо индекса означает «все действительные индексы»).

Аналогично, np.trapz может интегрировать все ломтиков время сразу:

P1 = np.trapz(Pd1, R1, axis=1) 
P2 = np.trapz(Pd2, R2, axis=1) 

Каждый P1 и P2 теперь временной ряд интегралов. Параметр axis определяет, по какой оси Pd1 интегрируется - это вторая ось, т. Е. Пробел.