原文：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 randomrseed=random.randint(0,10000)np.random.seed(rseed)n=int(n)low_x,high_x=xlow_y,high_y=yS=(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 randomrseed=random.randint(0,10000)np.random.seed(rseed)n=int(n)low_x,high_x=xlow_y,high_y=yS=(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 + 1a,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 resulta, 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
