Monte Carlo integration in case of coordinate transformation with Python

Question

I am trying to evaluate the integral using the Monte Carlo method where the integrand underwent a transformation from cylindrical to cartesian coordinates. The integrand itself is quite simple and could be calculated using scipy.integrate.quad, but I need it to be in cartesian coordinates for a specific purpose later on.

So here is the integrand: rho*k_r**2*(k_r**2*kn(0,k_r*rho)**2+k_n**2*kn(1,k_r*rho)**2) d rho

Here the kn(i,rho) is modified Bessel function of 2nd kind.

Solving it with quad gives the following result:

from scipy.special import kn
from scipy.integrate import quad
import random
k_r        = 6.2e-2
k_n        = k_r/1.05
C_factor   = 2*np.pi*1e5
lmax,lmin  = 80,50

def integration_polar():
    def K_int(rho):
        return rho*k_r**2*(k_r**2*kn(0,k_r*rho)**2+k_n**2*kn(1,k_r*rho)**2)
   
    rho = np.linspace(lmin,lmax,200)
    I,_ = quad(K_int,lmin,lmax)
    Gamma = I*C_factor
    print("expected=",Gamma)

Output: Expected = 7.641648442007296

Now the same integral using Monte Carlo method (hit-or-miss method looked up from here) gives almost the same result:

def integration_polar_MC():    
    random.seed(1)
    n = 100000
    def K_int(rho):
        return rho*k_r**2*(k_r**2*kn(0,k_r*rho)**2+k_n**2*kn(1,k_r*rho)**2)
    def sampler():
        x = random.uniform(lmin,lmax)
        y = random.uniform(0,c_lim)
        return x,y

    c_lim = 2*K_int(50) #Upper limit of integrand
    sum_I = 0
    for i in range(n):
        x,y = sampler()
        func_Int = K_int(x)
        if y>func_Int:
            I = 0
        elif y<=func_Int:
            I = 1
        sum_I += I
     Gamma = C_factor*(lmax-lmin)*c_lim*sum_I/n
     print("MC_integral_polar:",Gamma)

Output: MC_integral_polar = 7.637391399699502

Since Monte Carlo worked with this example, I thought the cartesian case would go smoothly as well but I couldn't get the right answer.

For the cartesian case, similarly as in previous case I've employed the hit-or-miss method, with rho = np.sqrt(x**2+y**2) and integrand becoming k_r**2*(k_r**2*kn(0,k_r*rho)**2+k_n**2*kn(1,k_r*rho)**2) dx dy where domain over x and y:

-80 <= x <= 80
-80 <= y <= 80
 50 <= np.sqrt(x**2+y**2) <= 80

Here is my attempt:

def integration_cartesian_MCtry():          
    random.seed(1)
    lmin,lmax = -100,100
    n = 100000
    def K_int(x,y):  
        rho = np.sqrt(x**2+y**2)
        if rho>=50 and rho<=80:
            return k_r**2*(k_r**2*kn(0,k_r*rho)**2+k_n**2*kn(1,k_r*rho)**2)
        else:
            return 0
    def sampler():
        x = random.uniform(lmin,lmax)
        y = random.uniform(lmin,lmax)
        z = random.uniform(0,c_lim)
        return x,y,z

    c_lim = K_int(50,0)
    sum_I = 0
    for i in range(n):
        x,y,z = sampler()
        func_Int = K_int(x,y)
        if z>func_Int:
            I = 0
        elif z<=func_Int:
            I = 1
        sum_I += I

    Gamma  = C_factor*(lmax-lmin)**2*c_lim*sum_I/n
    print("MC_integral_cartesian:",Gamma)

Output: MC_integral_cartesian = 48.83166430996952

As you can see Monte Carlo in cartesian overestimates the integral. I am not sure why it is happening but think that it may be related to the incorrect limits or domain over which I should integrate the function.

Any help appreciated as I am stuck without any progress for a few days.

Answer 1

Problem, as I said, is with jacobian. In case of polar, you have integration over

f(ρ)*ρ*dρ*dφ

You integrate over dφ analytically (your f(ρ) doesn't depend on φ), and get 2π

In case of cartesian there are no analytical integration, so it is over dx*dy, no factor of 2π. Code to illustrate it, Python 3.9.1, Windows 10 x64, and it produced pretty much the same answer

import numpy as np
from scipy.special import kn

k_r        = 6.2e-2
k_n        = k_r/1.05
C_factor   = 2*np.pi*1e5

lmin       = 50
lmax       = 80

def integration_polar_MC(rng, n):

    def K_int(rho):
        if rho>=50 and rho<=80:
            return rho*k_r**2*(k_r**2*kn(0, k_r*rho)**2 + k_n**2*kn(1, k_r*rho)**2)
        return 0.0

    def sampler():
        x = rng.uniform(lmin, lmax)
        y = rng.uniform(0.0, c_lim)
        return x,y

    c_lim = 2*K_int(50) # Upper limit of integrand
    sum_I = 0
    for i in range(n):
        x,y = sampler()
        func_Int = K_int(x)
        I = 1
        if y>func_Int:
            I = 0

        sum_I += I

    Gamma = C_factor*(lmax-lmin)*c_lim*sum_I/n
    return Gamma


def integration_cartesian_MC(rng, n):

    def K_int(x,y):
        rho = np.hypot(x, y)
        if rho>=50 and rho<=80:
            return k_r**2*(k_r**2*kn(0,k_r*rho)**2+k_n**2*kn(1,k_r*rho)**2)

        return 0.0

    def sampler():
        x = rng.uniform(lmin,lmax)
        y = rng.uniform(lmin,lmax)
        z = rng.uniform(0,c_lim)
        return x,y,z

    lmin,lmax = -100,100
    c_lim = K_int(50, 0)
    sum_I = 0
    for i in range(n):
        x,y,z = sampler()
        func_Int = K_int(x,y)
        I = 1
        if z>func_Int:
            I = 0
        sum_I += I

    Gamma  = C_factor*(lmax-lmin)**2*c_lim*sum_I/n
    return Gamma/(2.0*np.pi) # to compensate for 2π in the constant

rng = np.random.default_rng()
q = integration_polar_MC(rng, 100000)

print("MC_integral_polar:", q)

q = integration_cartesian_MC(rng, 100000)
print("MC_integral_cart:", q)