4. 反应-扩散方程的反应速率常数辨识

任务目标

已知一维反应-扩散系统中浓度场 \(u(x,t)\) 的稀疏观测数据，控制方程为：

\[ \frac{\partial u}{\partial t} = D \frac{\partial^2 u}{\partial x^2} + k u(1-u) \]

其中扩散系数 \(D=0.01\) 已知，待反演参数为反应速率常数 \(k\)。要求：

1. 构造合成观测数据（正问题）；
2. 使用物理信息神经网络（PINN）反演 \(k\)；
3. 对比不同观测稀疏程度下的反演精度。

实现方案

采用全连接神经网络近似 \(u(x,t)\)，将反应-扩散方程残差作为物理约束加入损失函数。网络输入为 \((x,t)\)，输出为 \(\hat{u}\)。总损失由三部分构成：

• 数据损失：网络输出与观测值的均方误差；
• 物理损失：在配点处计算方程残差；
• 边界/初始损失：满足 \(u(0,t)=u(1,t)=0\)，\(u(x,0)=\sin(\pi x)\)。

反演参数 \(k\) 作为可训练变量，与网络权重联合优化。通过控制训练数据的时间采样密度，考察反演精度变化。

关键完整代码

import numpy as np
import torch
import torch.nn as nn
import matplotlib.pyplot as plt

# ---------- 正问题：生成合成观测数据 ----------
def solve_reaction_diffusion(D, k, nx=101, nt=200, T=1.0):
    dx = 1.0 / (nx - 1)
    dt = T / nt
    x = np.linspace(0, 1, nx)
    u = np.sin(np.pi * x)
    u_history = [u.copy()]
    for _ in range(nt):
        u_new = u.copy()
        u_new[1:-1] += D * dt / dx**2 * (u[2:] - 2*u[1:-1] + u[:-2])
        u_new[1:-1] += dt * k * u[1:-1] * (1 - u[1:-1])
        u = u_new
        u_history.append(u.copy())
    return x, np.array(u_history)

# ---------- PINN 模型 ----------
class PINN(nn.Module):
    def __init__(self, n_input=2, n_hidden=50, n_output=1):
        super().__init__()
        self.net = nn.Sequential(
            nn.Linear(n_input, n_hidden), nn.Tanh(),
            nn.Linear(n_hidden, n_hidden), nn.Tanh(),
            nn.Linear(n_hidden, n_hidden), nn.Tanh(),
            nn.Linear(n_hidden, n_output)
        )
        self.k = nn.Parameter(torch.tensor(0.5, dtype=torch.float32))  # 待反演参数

    def forward(self, x, t):
        return self.net(torch.cat([x, t], dim=1))

# ---------- 训练配置 ----------
D_true = 0.01
k_true = 2.0
x, u_full = solve_reaction_diffusion(D_true, k_true)

# 构造训练数据（不同稀疏度）
def prepare_data(time_stride):
    idx_t = np.arange(0, 201, time_stride)
    idx_x = np.arange(0, 101, 5)
    X_mesh, T_mesh = np.meshgrid(x[idx_x], np.arange(0, 1.0+1e-8, 1.0/200)[idx_t])
    u_obs = u_full[idx_t][:, idx_x]
    x_obs = torch.tensor(X_mesh.reshape(-1, 1), dtype=torch.float32)
    t_obs = torch.tensor(T_mesh.reshape(-1, 1), dtype=torch.float32)
    u_obs_t = torch.tensor(u_obs.reshape(-1, 1), dtype=torch.float32)
    return x_obs, t_obs, u_obs_t

# ---------- 训练与反演 ----------
def train_pinn(time_stride, n_iter=3000, lr=1e-3):
    x_obs, t_obs, u_obs = prepare_data(time_stride)
    model = PINN()
    optimizer = torch.optim.Adam(model.parameters(), lr=lr)

    # 物理配点（拉丁超立方采样）
    n_colloc = 5000
    x_c = torch.rand(n_colloc, 1)
    t_c = torch.rand(n_colloc, 1) * 1.0

    # 边界/初始配点
    x_b0 = torch.zeros(500, 1)
    x_b1 = torch.ones(500, 1)
    t_b = torch.rand(500, 1) * 1.0
    x_init = torch.rand(500, 1)
    t_init = torch.zeros(500, 1)
    u_init = torch.sin(np.pi * x_init)

    for it in range(n_iter):
        optimizer.zero_grad()

        # 数据损失
        u_pred = model(x_obs, t_obs)
        loss_data = nn.MSELoss()(u_pred, u_obs)

        # 物理损失（自动微分）
        x_c.requires_grad_(True)
        t_c.requires_grad_(True)
        u_c = model(x_c, t_c)
        u_t = torch.autograd.grad(u_c, t_c, grad_outputs=torch.ones_like(u_c), create_graph=True)[0]
        u_x = torch.autograd.grad(u_c, x_c, grad_outputs=torch.ones_like(u_c), create_graph=True)[0]
        u_xx = torch.autograd.grad(u_x, x_c, grad_outputs=torch.ones_like(u_x), create_graph=True)[0]
        f_res = u_t - D_true * u_xx - model.k * u_c * (1 - u_c)
        loss_phys = torch.mean(f_res**2)

        # 边界/初始损失
        u_b0 = model(x_b0, t_b)
        u_b1 = model(x_b1, t_b)
        u_i = model(x_init, t_init)
        loss_bc = torch.mean(u_b0**2) + torch.mean(u_b1**2) + nn.MSELoss()(u_i, u_init)

        loss = loss_data + 0.1 * loss_phys + loss_bc
        loss.backward()
        optimizer.step()

        if it % 500 == 0:
            print(f"iter {it}, loss={loss.item():.2e}, k={model.k.item():.4f}")

    return model.k.item()

# ---------- 运行不同稀疏度实验 ----------
print("=== 时间步长=1（密集） ===")
k_est_dense = train_pinn(time_stride=1)
print(f"反演 k = {k_est_dense:.4f} (真值 2.0)")

print("\n=== 时间步长=10（稀疏） ===")
k_est_sparse = train_pinn(time_stride=10)
print(f"反演 k = {k_est_sparse:.4f} (真值 2.0)")

运行与验证

依赖：Python 3.8+，PyTorch 1.12+，NumPy，Matplotlib。

运行方式：直接执行上述脚本。正问题使用显式有限差分生成参考解，时间步长 \(\Delta t=0.005\)，空间步长 \(\Delta x=0.01\)，满足数值稳定性条件。

验证方法：

1. 观察训练日志中损失是否下降至 \(10^{-4}\) 量级以下；
2. 比较反演得到的 \(k\) 与真值 2.0 的相对误差；
3. 对比 time_stride=1（201个时间帧）与 time_stride=10（21个时间帧）的反演结果，稀疏数据下误差应明显增大。

预期输出：密集数据下 \(k\) 估计值接近 2.0（误差通常 <5%）；稀疏数据下误差增大至 10%～30% 范围，具体数值因网络初始化与训练随机性略有波动。

数据稀疏性影响讨论

时间采样间隔从 \(\Delta t=0.005\) 增大至 \(\Delta t=0.05\) 时，训练数据量减少为原来的 1/10。物理损失项在总损失中占比相对提升，但配点数量不变，导致网络更依赖方程约束而缺乏足够观测锚定，反演参数方差增大。实际应用中若观测时间分辨率不足，可考虑增加配点密度或采用时间分段训练策略缓解。

上图定性展示了时间采样间隔增大时反演误差的上升趋势。实际数值会因网络结构、训练迭代次数及随机种子不同而波动，但稀疏数据导致精度下降的规律是稳定的。