计算科学 - 一个分量的高斯-厄米特双积分 - 吾爱随笔录

我正在尝试执行以下积分

\int_{0}^{2 π} \int_{0}^{+ \infty} \frac{r^{'} (e^{- r^{' 2} / 2 σ^{2}}) (r - r^{'} \cos (θ - θ^{'}))}{r^{2} + r^{' 2} - 2 r r^{'} \cos (θ - θ^{'})} d r^{'} d θ^{'}

$\int_{0}^{2\pi}\int_{0}^{+\infty} \frac{r'\left(e^{-r'^2/2\sigma^2}\right)\left(r-r'\cos(\theta-\theta')\right)}{r^2+r'^2-2rr'\cos(\theta-\theta')}dr'dθ'$

使用 Gauss-Hermite，使用Simpson 1/3 规则，但没有成功。我找不到我的错误，但输出应该如图 2 所示。这是我的代码（抱歉我的格式错误，这是我第一次在这里上传）。 $r$ $\theta$

$\sigma$ 应假定为 1。

import numpy as np
import matplotlib.pyplot as plt
import scipy.special as ss

def rt(d, r, theta ,sig):
    return r*(d-r*np.cos(theta))*np.exp(-r**2/(2*sig**2))/(d**2+r**2-2*d*r*np.cos(theta))
def intheta1(d, r, b, sig, N):
    h = b/N
    I = rt(d,r,0,sig) + rt(d,r,b,sig)
    for i in range(1, N, 2):
        I += 4*rt(d, r, i*h, sig)
    for j in range(2, N, 2):
        I += 2*rt(d, r, j*h, sig)
    return I*h/3

def intr1(d, b, sig, N, M):
    x, w = ss.roots_hermitenorm(N)
    s = 0
    for k in range(N):
        s += intheta1(d, x[k], b, sig, M)*w[k]
    return s/2

ps = np.linspace(0, 5, 1000)
qs = intr1(xs, 2*np.pi, 1, 1000, 90)

plt.plot(ps, qs)

python def integral_theta(r, rline, theta, sigma): rline = np.abs(rline) return np.exp((-(rline)**2*(1-sigma**2))/2/sigma**2) * rline * (r - rline*np.cos(theta))/(r**2 + rline**2 - 2*r*rline*np.cos(theta)) def i_theta(r, rline, sigma): a, b = 0, 2*np.pi N = 100 h = (b-a)/N s_odd = 0 for k in range(1,N,2): s_odd += integral_theta(r, rline, a+k*h, sigma) s_even = 0 for j in range(2, N-1,2): s_even += integral_theta(r, rline, a+j*h, sigma) return h/3*(integral_theta(r, rline, a, sigma) + integral_theta(r, rline, b, sigma) + 4*s_odd + 2*s_even) def i_r(r, sigma): M = 1000 x, w = ss.roots_hermitenorm(M) s = 0 for h in range(M): s += i_theta(r, x[h], sigma)*w[h] return s/2/np.pi/sigma**2