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)
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. $\endgroup$np.exp(-r**2/(2*sig**2))
and as far as I understand you used
in your code for showing $r^{'}$, but it's clearly in conflict with your formula. So, something is wrong here for sure... $\endgroup$