Monte Carlo坐标变换情况下的积分Python

Question

我正在尝试使用 Monte Carlo 方法计算积分，其中被积函数经历了从圆柱坐标到笛卡尔坐标的转换。被积函数本身非常简单，可以使用 scipy.integrate.quad 计算，但我需要它在笛卡尔坐标中用于稍后的特定目的。

所以这是被积函数：rho*k_r**2*(k_r**2*kn(0,k_r*rho)**2+k_n**2*kn(1,k_r*rho)**2) d rho

这里的kn(i,rho)是第二类修正贝塞尔函数

用quad求解得到以下结果：

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)

输出：Expected = 7.641648442007296

现在使用 Monte Carlo 方法（从 here 查找的命中或未命中方法）的相同积分给出几乎相同的结果：

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)

输出：MC_integral_polar = 7.637391399699502

自从 Monte Carlo 使用此示例后，我认为笛卡尔的情况也会顺利进行，但我无法得到正确的答案。

对于笛卡尔情况，与之前的情况类似，我采用了命中或未命中方法，其中 rho = np.sqrt(x**2+y**2) 和被积函数变为 k_r**2*(k_r**2*kn(0,k_r*rho)**2+k_n**2*kn(1,k_r*rho)**2) dx dy，其中域超过 x 和 y:

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

这是我的尝试：

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)

输出：MC_integral_cartesian = 48.83166430996952

如您所见，笛卡尔坐标中的 Monte Carlo 高估了积分。我不确定为什么会这样，但认为这可能与我应该集成该功能的不正确限制或域有关。

感谢任何帮助，因为我被困了几天没有任何进展。

Answer 1

正如我所说，问题出在 jacobian 上。在极地的情况下，你有超过

的积分

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

你对 dφ 进行解析积分（你的 f(ρ) 不依赖于 φ），得到 2π

如果是笛卡尔坐标则没有解析积分，所以它超过了dx*dy，没有因数 2π。说明它的代码，Python 3.9.1，Windows 10 x64，它产生了几乎相同的答案

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)

Monte Carlo坐标变换情况下的积分Python

Monte Carlo integration in case of coordinate transformation with Python

python

montecarlo