Python中的概率分佈函數（PDF）加權的數值雙積分

Question

我在執行函數中的雙定積分時出現問題，該函數取決於2個變量（q,r），並且有一個額外的積分。 我想用一個高斯函數與重量的功能是：

F(q,r)=f(q,r)+int_{0,r}(h(q,r')dr')

而在必須重新整合與高斯加權：

I(q)=int_{0,inf}(F(q,r)^2*g(r)dr)

高斯g(r)在中心座標R。

您可以觀察到的主要問題是我正在將數組與標量混合。使用與高斯相同的方法（np.ogrid和求和軸）可能是一個解決方案，但我不知道如何實現它。

import numpy as np 
from scipy.integrate import quad 
import math as m 

R=53. 
R0=40. 
delta=50. 
c=2. 
qm, rm = np.ogrid[0.0005:2.0:0.0005, 20:100:500j] 

#normalized gauss function 
#g(r) 
def gauss_grid(r,Rmin,pd): 
    def gauss(r,Rmin,pd): 
     sigma=1.5 
     return (1/sigma)*np.exp(-((r-Rmin)**2)/(2*sigma**2)) 
    gauss_grid = gauss(r,Rmin,pd) 
    #normalization of gaussian 
    gauss_grid /= np.sum(gauss_grid) 
    return gauss_grid 

#spherical function 
#f(q,r) 
def form(q,R): 
    return (4/3)*m.pi*3*(np.sin(q*R)-q*R*np.cos(q*R))/(q**3) 

#FINAL function 
#I(q) 
def helfand(): 
    def F(q,R): 
     #integral (0,R) of h(q,r) 
     def integral(q,Rmax): 
      #h(q,r) 
      def integrand(r,q): 
       return np.sin(q*r)*(r**2)/(q*r*(1+np.exp(c*(R0-r)))) 
      return quad(integrand, 0, Rmax, args=(q))[0] 
     return (form(q,R)+delta*integral(q,R))**2 

    FF_hel=F(qm,rm) 
    FF_hel *= gauss_grid(rm,R,pd) 
    I=FF_hel.sum(axis=1) 
    return I,qm.ravel() 

helfand() 

* UPDATE * * **

我試圖與scipy.integrate庫（與quad），我不能讓它做。這就像它沒有將正確的參數（q）傳遞給下一個函數。這裏我想要一個非常簡化的版本：

import numpy as np 
from scipy.integrate import quad 
import matplotlib.pyplot as plt 

R=53. 
R0=41. 
pd=15. 
sigma=1.5 

def I(q): 
    #(function with integral inside) squared 
    def FF(q,r): 
     def integral_f(q,r): 
      def f(r1,q): 
       return np.sin(q*r1) 
      return quad(f,0,r,args=(q))[0] 

     def h(q,r): 
      return (r*np.cos(q*r)) 
     return (h(q,r)+integral_f(q,r))**2 

    #gaussian function normalized 
    def g(r,R0): 
     def gauss(r,R0): 
      return (1/sigma)*np.exp(-((r-R0)**2)/(2*sigma**2)) 
     return gauss(r,R0)/(quad(gauss,0,np.inf,args=(R0))[0]) 

    #main function to be integrated with gaussian 
    def function(r,q): 
     return FF(q,r)*g(r,R) 

    return quad(function,0,np.inf,args=(q))[0] 

q=np.arange(0.001,1.,0.001) 
plt.plot(q,I(q)) 

錯誤說：

提供的函數不返回一個有效的浮動。

Answer 1

我會創建一個簡單的2D矩形網格點，跨越積分點的極限。然後，我寧願在每個元素上進行高斯積分來評估積分。這意味着在每個積分點調用函數，加權或不加權，乘以正交權重和求和。

它與二維四邊形有限元相似，並通過數值積分來評估剛度矩陣。

在SciPy中有2D正交方法。我會在寫我自己的之前使用它。

http://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.quadrature.html

Answer 2

我認爲你可以用兩個單積分計算這一點。

如果你寫出來的雙重積分，你得到兩個部分：

int_ {0，INF}（F（Q，R）*克（R）DR）+ int_（0，INF）（int_ { 0，R}（H（q，R'）博士）* G（R），DR）

我們可以在第二交換整合以獲得

int_（0，INF）（int_ {r'，inf}（g（r）dr）* h（q，r'）dr'）

內積分可以用互補的 誤差函數表示。