如何迭代森林火灾模拟函数并整合代码实现密度计算？

以下是修正并整合后的完整代码，解决迭代逻辑、状态转换与数据可视化问题：

完整可运行代码 #

import numpy as np
import matplotlib.pyplot as plt

# 初始参数设置
p, f = 0.5, 0.3
nx, ny = 100, 100
# 初始化森林：70%空地块，30%树木，随机设置1个初始火源
X = np.random.choice([0, 1], size=(ny, nx), p=[0.7, 0.3])
X[np.random.randint(ny), np.random.randint(nx)] = 2

adjacent = ((-1,0), (0,-1), (0, 1), (1,0))
E, T, F = 0, 1, 2
xvalues = [0]
yvalues = [0]

# 修正后的状态迭代函数
def iterate(X):
    Xnew = np.copy(X)
    for ix in range(1, nx-1):
        for iy in range(1, ny-1):
            current = X[iy, ix]
            if current == E:
                # 空地块随机长出树木
                if np.random.random() <= p:
                    Xnew[iy, ix] = T
            elif current == T:
                # 检查是否被相邻火源点燃
                has_fire_neighbor = False
                for dx, dy in adjacent:
                    if X[iy+dy, ix+dx] == F:
                        has_fire_neighbor = True
                        break
                if has_fire_neighbor:
                    Xnew[iy, ix] = F
                else:
                    # 随机自然着火
                    if np.random.random() <= f:
                        Xnew[iy, ix] = F
            elif current == F:
                # 燃烧后变为空地块
                Xnew[iy, ix] = E
    return Xnew

# 执行最多1000次迭代
max_iterations = 1000
for i in range(1, max_iterations + 1):
    X = iterate(X)
    # 计算各地块数量
    num_empty = (X == E).sum()
    num_tree = (X == T).sum()
    num_fire = (X == F).sum()
    # 计算密度（避免分母为0）
    denominator = num_fire + num_empty
    density = num_tree / denominator if denominator != 0 else 0
    # 存储数据
    xvalues.append(i)
    yvalues.append(density)
    # 打印当前迭代密度
    print(f"Iteration {i}: Density = {density:.4f}")

# 绘制树木密度变化图
plt.figure(figsize=(10, 6))
plt.plot(xvalues, yvalues, label='树木密度')
plt.xlabel('迭代次数')
plt.ylabel('树木密度（树木数/(空地块+燃烧地块)）')
plt.title('森林火灾模拟：树木密度随迭代变化')
plt.legend()
plt.grid(True)
plt.show()

关键修正与说明 #

• 修复迭代函数逻辑：
  • 原代码中else与for循环绑定的逻辑错误，改为显式检查相邻火源
  • 补充燃烧地块（F）的状态转换：燃烧后变为空地块（E），符合森林火灾模型常规逻辑
  • 使用np.copy(X)初始化新状态数组，避免原数组被提前修改导致计算错误
• 初始条件优化：原代码全为空地块无法触发火灾，添加了随机树木与初始火源生成
• 迭代循环实现：通过for循环执行最多1000次迭代，每次更新状态并计算密度，存储数据到列表
• 异常处理：添加分母为0的判断，避免极端情况下的除以0错误

内容的提问来源于stack exchange，提问作者tho