使用蒙特卡洛方法求解一阶微分方程

发表于 更新于

使用蒙特卡洛方法求解一阶线性常微分方程

在这个例子中，我们将展示如何使用蒙特卡洛方法来近似求解一个简单的一阶线性常微分方程，并与解析解进行比较。

问题描述

考虑如下形式的初值问题：

\[\frac{dy}{dx}=-2xy,\quad y(0)= 1\]

其解析解为 \(y(x) = e^{-x^2}\)。我们将使用蒙特卡洛方法来近似这个解。

Python代码实现

以下是基于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
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
import numpy as np
import matplotlib.pyplot as plt

def monte_carlo_ode_solver(f, x0, y0, x_end, num_samples=1000, dx=0.05):
    """
    使用蒙特卡洛方法求解一阶常微分方程。
    
    参数:
        f: 函数，定义了dy/dx = f(x, y)
        x0 (float): 初始x值。
        y0 (float): 初始y值。
        x_end (float): x的最大值。
        num_samples (int): 每一步中使用的样本数量。
        dx (float): x方向上的步长。
    返回:
        x_points (numpy.ndarray): x坐标点。
        y_estimated (numpy.ndarray): 对应每个x点估计的y值。
    """
    x_points = np.arange(x0, x_end+dx, dx)
    y_estimated = np.zeros_like(x_points)
    y_estimated[0] = y0  # 设置初始条件
    
    for i in range(1, len(x_points)):
        current_x = x_points[i-1]
        current_y = y_estimated[i-1]
        
        dy_samples = []
        for _ in range(num_samples):
            # 随机扰动dx
            delta_x = dx * (np.random.rand() - 0.5)
            # 根据当前点的斜率计算dy
            delta_y = f(current_x, current_y) * (dx + delta_x)
            dy_samples.append(delta_y)
        
        # 计算平均dy并更新y值
        avg_dy = np.mean(dy_samples)
        y_estimated[i] = current_y + avg_dy
    
    return x_points, y_estimated

# 定义导数函数
def f(x, y):
    return -2 * x * y

# 运行求解器
x, y_mc = monte_carlo_ode_solver(f, x0=0, y0=1, x_end=2, num_samples=1000)

# 计算真实解用于比较
y_true = np.exp(-x**2)

# 绘制结果
plt.plot(x, y_mc, label='Monte Carlo Approximation')
plt.plot(x, y_true, label='True Solution', linestyle='--')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.show()

分析

在使用蒙特卡洛方法求解微分方程时，最关键的一步是 delta_x = dx * (np.random.rand() - 0.5), 这表示对于步长\(dx\)在其中心附近扰动，并且保存每次扰动的所造成的\(dy\)值，当完成大量的数据采样时，根据大数定理，\(dy\)值会稳定到真实的\(dy\)的情况。进一步将此值追加到\(y\)，同时\(dx\)追加到\(x\)后，再进行下一个点的判断。所以蒙特卡洛方法基于大数定理是一个有效的方法。