使用蒙特卡洛方法求解定积分

发表于 更新于

基本原理

蒙特卡洛方法求解定积分的理论依据是积分中值定理和大数定理。即

\[\int_a^bf(x)dx=f(\xi)\cdot(b-a)\]

其中\(a\leq\xi\leq b\). 蒙特卡洛方法是随机在\((a,b)\)区间取点，然后将这个点当作\(\xi\)来计算面积，但是由于点\(\xi\)当前是随机的，这个积分值就不是准确的，但是当大量取点时，由于函数\(f(x)\)本身的作用，根据大数定量最终会稳定到准确值的附近，因为接近真实值的情况出现的概率更大。

一维函数的定积分

假设我们想要计算如下函数在区间 \([0, 1]\) 上的积分:

\[\int_0^1x^2dx\]

该积分的解析解是\(\frac{1}{3}\). 使用Python代码实现蒙特卡洛积分的代码为:

一维函数积分
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
import numpy as np

def monte_carlo_integral(f, a, b, N):
    """
    使用蒙特卡洛方法估计函数f在[a,b]上的积分。
    
    参数:
        f: 要积分的函数。
        a (float): 积分下限。
        b (float): 积分上限。
        N (int): 随机样本的数量。
    返回:
        float: 积分的蒙特卡洛估计值。
    """
    samples = np.random.uniform(a, b, N)  # 在[a, b]上生成N个均匀分布的随机数
    evaluations = f(samples)  # 计算每个随机样本点处函数的值
    integral_estimate = (b - a) * np.mean(evaluations)  # 积分估计
    return integral_estimate

# 定义要积分的函数
def f(x):
    return x**2

# 运行蒙特卡洛积分
a, b = 0, 1  # 积分范围
N = 100000   # 样本数量
integral = monte_carlo_integral(f, a, b, N)

print(f"蒙特卡洛积分估计值: {integral}")
print(f"真实值: {1/3}")

二维函数的定积分

二维函数的定积分同样也满足积分中值定理，所以蒙特卡洛方法是求解定积分的有效工具。本节我们求解定义在\([0,1]\times[0,1]\)上的二重积分:

\[\iint_{[0,1]\times[0,1]]}(x^2+y^2)dxdy\]

该积分的解可以直接解析求得为\(\frac{2}{3}\). 其 Python 代码实现为:

二维函数定积分
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
import numpy as np

def monte_carlo_double_integral(f, a1, b1, a2, b2, N):
    """
    使用蒙特卡洛方法估计函数f在[a1,b1] x [a2,b2]上的二重积分。
    
    参数:
        f: 要积分的函数，接受两个参数x和y。
        a1, b1 (float): 第一维的积分下限和上限。
        a2, b2 (float): 第二维的积分下限和上限。
        N (int): 随机样本的数量。
    返回:
        float: 积分的蒙特卡洛估计值。
    """
    # 在[a1, b1]和[a2, b2]上生成N个均匀分布的随机数对
    x_samples = np.random.uniform(a1, b1, N)
    y_samples = np.random.uniform(a2, b2, N)
    
    # 计算每个随机样本点处函数的值
    evaluations = f(x_samples, y_samples)
    
    # 积分估计
    area = (b1 - a1) * (b2 - a2)
    integral_estimate = area * np.mean(evaluations)
    
    return integral_estimate

# 定义要积分的函数
def integrand(x, y):
    return x**2 + y**2

# 运行蒙特卡洛二重积分
a1, b1 = 0, 1  # 第一维积分范围
a2, b2 = 0, 1  # 第二维积分范围
N = 100000     # 样本数量

integral = monte_carlo_double_integral(integrand, a1, b1, a2, b2, N)

print(f"蒙特卡洛积分估计值: {integral}")
print(f"真实值: {2/3}")

通过增加样本数量 N，可以提高积分估计的精度。尽管对于这个简单的二重积分问题，传统数值积分方法可能更高效且准确，但对于更高维度或复杂形状的积分区域，蒙特卡洛方法就显得尤为重要了。这种方法提供了一种通用的方法来解决那些难以通过常规手段求解的积分问题。