蒙特卡洛方法求解偏微分方程Python代码

蒙特卡洛方法是一种基于随机抽样的数值计算方法，尤其适用于解决那些难以解析求解的复杂问题，包括部分偏微分方程（PDEs）。在Python中，我们可以使用像`scipy`、`numpy`和`matplotlib`这样的库来实现这个过程。下面是一个简单的例子，展示如何用Python的蒙特卡洛模拟求解二维扩散方程：

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import multivariate_normal

# 定义扩散系数和时间步长
D = 0.1
dt = 0.01
n_steps = 1000
grid_size = 100
grid = np.linspace(0, 1, grid_size)

# 创建一个二维网格
X, Y = np.meshgrid(grid, grid)

# 初始化随机噪声，模拟起始分布
random_noise = multivariate_normal(mean=[0.5, 0.5], cov=2 * D * dt)

# 模拟步骤
for _ in range(n_steps):
    # 从当前分布中生成新的随机点并应用扩散
    new_points = (X + D * dt * np.random.randn(*X.shape),
                  Y + D * dt * np.random.randn(*Y.shape))
    random_noise = random_noise + new_points - np.roll(np.roll(random_noise, shift=-1, axis=0), shift=-1, axis=1)

# 计算最终分布
final_distribution = random_noise / (2 * D * dt * n_steps)

# 可视化结果
plt.imshow(final_distribution, cmap='viridis', origin='lower')
plt.colorbar()
plt.xlabel('x')
plt.ylabel('y')
plt.title('Monte Carlo Solution of 2D Diffusion Equation')
plt.show()