# montecarlo - Monte Carlo integration in case of coordinate transformation with Python

299 views

### montecarlo - Monte Carlo integration in case of coordinate transformation with Python

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
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)
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.

by (71.8m points)

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)
``````