【数值计算方法】蒙特卡洛方法积分的Python实现

原文：https://www.cnblogs.com/aksoam/p/18378332

原理不做赘述,参见【数值计算方法】数值积分&微分-python实现 - FE-有限元鹰 - 博客园,直接上代码,只实现1d,2d积分,N维积分的蒙特卡洛方法也类似.

• 代码

from typing import Callable,Union,List
def MonteCarloInt2d(f:Callable,x:Union[float,List[float]],
                    y:Union[float,List[float]],n:int=10000)->float:
    """蒙特卡洛2d积分"""
    
    import random
    rseed=random.randint(0,10000)
    np.random.seed(rseed)
    n=int(n)
    low_x,high_x=x
    low_y,high_y=y
    
    S=(high_x-low_x)*(high_y-low_y)
    # 生成随机浮点数在[low_x,high_x]和[low_y,high_y]之间的坐标
    xys=np.random.uniform(low=[low_x,low_y], high=[high_x,high_y], size=(n,2))
    return (S/n)*sum([f(xys[ind,0],xys[ind,1]) for ind in range(xys.shape[0])])

def MonteCarloInt2d(f:Callable,x:Union[float,List[float]],
                    y:Union[float,List[float]],n:int=10000)->float:
    """蒙特卡洛2d积分"""
    
    import random
    rseed=random.randint(0,10000)
    np.random.seed(rseed)
    n=int(n)
    low_x,high_x=x
    low_y,high_y=y
    
    S=(high_x-low_x)*(high_y-low_y)
    # 生成随机浮点数在[low_x,high_x]和[low_y,high_y]之间的坐标
    xys=np.random.uniform(low=[low_x,low_y], high=[high_x,high_y], size=(n,2))
    return (S/n)*sum([f(xys[ind,0],xys[ind,1]) for ind in range(xys.shape[0])])

• 测试: 1d积分

import scipy
def f(x):
    return x**2 + 2*x + 1

a,b=-1,2
v=MonteCarloInt1d(f,a,b,1e6)
print(f"integral of f(-1,2) is {v:.4f}")

scipy_result,error = scipy.integrate.quad(f, a, b)
print(f"SciPy积分结果: {scipy_result},error: {error}")

# <!-- 输出结果 -->
# integral of f(-1,2) is 8.9909
# SciPy积分结果: 9.0,error: 9.992007221626409e-14

• 测试: 2d积分

from scipy.integrate import dblquad
def f(x, y):
    return x**2 + y**2

# 使用SciPy计算2d积分
def scipy_integrate_2d(f, a, b, c, d):
    result, error = dblquad(f, a, b, lambda x: c, lambda x: d)
    return result

a, b = 0, 1
c, d = 0, 1

# 使用蒙特卡洛方法计算2d积分
monte_carlo_result = MonteCarloInt2d(f, [a,b],[ c, d])
print(f"蒙特卡洛2d积分结果: {monte_carlo_result}")

# 使用SciPy计算2d积分
scipy_result = scipy_integrate_2d(f, a, b, c, d)
print(f"SciPy 2d积分结果: {scipy_result}")

# <!-- 输出结果 -->
# 蒙特卡洛2d积分结果: 0.6676998082037032
# SciPy 2d积分结果: 0.6666666666666669