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)

My Result

The actual result

  • $\begingroup$ 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. $\endgroup$ Apr 6, 2020 at 16:13
  • $\begingroup$ Did my best, hopefully it is a little more understandable. $\endgroup$ Apr 6, 2020 at 16:28
  • $\begingroup$ 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... $\endgroup$ Apr 6, 2020 at 16:31
  • $\begingroup$ 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) $\endgroup$ Apr 6, 2020 at 16:42
  • 1
    $\begingroup$ 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 $\endgroup$ Apr 6, 2020 at 16:57

1 Answer 1


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.

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
  • 3
    $\begingroup$ 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". $\endgroup$
    – davidhigh
    Jan 3, 2021 at 20:20
  • $\begingroup$ Or search for "Gauss–Laguerre quadrature". $\endgroup$
    – cos_theta
    May 4, 2021 at 9:06
  • 2
    $\begingroup$ @cos_theta: but Gauss-Laguerre has a different weight function as the one in the OP, i.e. $e^{-r}$ instead of $e^{-r^2/2}$. $\endgroup$
    – davidhigh
    May 4, 2021 at 21:55

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.