Double Integral with Gauss- Hermite for one component

I am trying to perform the following integral

$$\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θ'$$

Using Gauss-Hermite for $$r$$ and Simpson 1/3 rule for $$\theta$$ with no success. I can't find my mistake but the output should look like Fig. 2. This was my code (sorry for my bad formatting, this is my first time uploading here).

$$\sigma$$ should be assumed as 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)


• I cannot match rt function in your Python script with your formula written above. Are you sure the formula is correctly implemented? Also, please remove the image and write your formula by using our LaTeX here. – Alone Programmer Apr 6 '20 at 16:13
• Did my best, hopefully it is a little more understandable. – Tomás Lopes Apr 6 '20 at 16:28
• Still I can't match it with your Python implementation. For example: in your formula you have: $\exp{(-\frac{(r^{'})^{2}}{2 \sigma^{2}})}$ but in your code you have: np.exp(-r**2/(2*sig**2)) and as far as I understand you use d in your code for showing $r^{'}$, but it's clearly in conflict with your formula. So, something is wrong here for sure... – Alone Programmer Apr 6 '20 at 16:31
• I am using r as r' and d as r....My goal is to make the graph of the magnetic field (the result of the Integral) as a function of r (represented by d in python) – Tomás Lopes Apr 6 '20 at 16:42
• Sorry but I don't think you can easily map a region of $[0,5]$ to $[-\infty,\infty]$ easily without a complex changing variable. More convenient way to do this is to use Gauss quadrature and map $[0,5]$ region to $[-1,1]$: en.wikipedia.org/wiki/Gaussian_quadrature – Alone Programmer Apr 6 '20 at 16:57

I have been studying this type of numerical integration and I believe I understood my mistake. First of all I am using Gauss-Hermite which work with limits $${-\infty}$$ to $${\infty}$$ so using the fact that this function is even makes it so that to integrate from $$0$$ to $${\infty}$$ I have to use np.abs() of my integration variable. Also, using Gauss-Hermite makes it so that I have to remove the exponential function. In this case I am using roots_hermitenorm() so I had to find a way to remove $$\exp(-r^2/2)$$ from the expression.
I got to these answer which is currently working flawlessly.

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

• No need to do any tricky stuff: There is a corresponding Gaussian quadrature with this weight function on the interval $[0,\infty]$. Just search for "half-open Hermite quadrature". – davidhigh Jan 3 at 20:20