当前位置：首页 > news >正文

PINN学习（三）—— 发现方程问题的解决

news 2026/5/12 19:15:46

所用部分图片来源于PINN学习第一部分的视频资料。

一、背景

前两节讲解了PINN对正、逆问题的求解，以一维热传导问题为例，并附上了相关的基础PINN代码帮助理解，接下来是对PINN的最后一个功能——发现方程进行一维热传导问题的解决。

个人认为PINN用于发现方程其实是基于对逆问题的求解进行的，我们的逆问题是使用PINN求解方程的一个系数，可以想到，如果将方程的项进行提前设定为u_t、u、u_x、u_xx、u_xxx...等项的组合，然后额外设定一个列表，列表包含各个系数，对这些系数进行PINN逆问题求解更新，如果合适，应该会得到一个贴近实际PDE的一个系数。

可以按这个图片理解，可以看出这样的PINN仍然是一个半监督学习。并且，因为一般认为PDE不会包含很多的微分项，即各个系数应该有很多的值为“0”，所以我们额外加入了L2范数作为系数的Loss，强制大多数系数为0。

保持系数矩阵的L2范数为0是十分重要的，这一点在我训练PINN时体现的淋漓尽致，如果对系数的约束不足，往往会导致PINN发现方程变成了类似多项式插值一样，虽然数据拟合效果好，但是方程完全不正确（方程包含了各阶微分项）。

二、代码实例

import torch import torch.nn as nn import numpy as np import matplotlib.pyplot as plt # ---------------------- 1. 解析解（生成真实数据） ---------------------- def exact_solution_heat(a_sq, x, t): """一维热传导方程解析解：u = e^(-απ²t) * sin(πx)""" if isinstance(x, torch.Tensor): u = torch.exp(-a_sq * (np.pi ** 2) * t) * torch.sin(np.pi * x) else: u = np.exp(-a_sq * (np.pi ** 2) * t) * np.sin(np.pi * x) return u # ---------------------- 2. 简单全连接网络 ---------------------- class SimpleFCN(nn.Module): def __init__(self): super().__init__() self.net = nn.Sequential( nn.Linear(2, 20), nn.Tanh(), nn.Linear(20, 20), nn.Tanh(), nn.Linear(20, 1) ) def forward(self, x): return self.net(x) # ---------------------- 3. 数据准备 ---------------------- torch.manual_seed(123) np.random.seed(123) # 真实热传导方程：u_t = 0.01*u_xx（仅含u_xx项） true_a_sq = 0.01 # 候选微分项库（预设可能的项，模拟"未知方程结构"） term_names = ['u', 'u_x', 'u_xx', 'u_xxx', 'u*t', 'sin(x)'] # 候选项 n_terms = len(term_names) # 候选项数量：6 # 1. 观测数据 n_obs = 50 x_obs = torch.rand(n_obs, 1) t_obs = torch.rand(n_obs, 1) X_obs = torch.cat([x_obs, t_obs], dim=1) u_obs = exact_solution_heat(true_a_sq, x_obs, t_obs) + 0.001 * torch.randn_like(x_obs) # 2. 物理配点（用于构建微分项库） n_f = 500 x_f = torch.rand(n_f, 1).requires_grad_(True) t_f = torch.rand(n_f, 1).requires_grad_(True) X_f = torch.cat([x_f, t_f], dim=1) # ---------------------- 4. 初始化：网络+微分项系数 ---------------------- model = SimpleFCN() # 可学习系数：每个候选微分项对应一个系数（初始值随机） term_coeffs = nn.Parameter(torch.randn(n_terms, 1), requires_grad=True) optimizer = torch.optim.Adam(list(model.parameters()) + [term_coeffs], lr=1e-3) epochs = 15000 # 训练轮数 # ---------------------- 5. 训练：发现微分项+优化系数 ---------------------- print("开始训练（目标：识别微分方程结构+求解系数）...") for epoch in range(epochs): optimizer.zero_grad() # 1. 数据损失：拟合观测数据 u_pred_obs = model(X_obs) loss_data = torch.mean((u_pred_obs - u_obs) ** 2) # 2. 构建候选微分项库（计算所有预设微分项） u_pred_f = model(X_f) # 计算各阶导数 grads1 = torch.autograd.grad(u_pred_f, X_f, grad_outputs=torch.ones_like(u_pred_f), create_graph=True)[0] u_t = grads1[:, 1:2] # 方程左侧：u_t u_x = grads1[:, 0:1] # 一阶x导数 u_xx = torch.autograd.grad(u_x, X_f, grad_outputs=torch.ones_like(u_x), create_graph=True)[0][:, 0:1] # 二阶x导数 u_xxx = torch.autograd.grad(u_xx, X_f, grad_outputs=torch.ones_like(u_xx), create_graph=True)[0][:, 0:1] # 三阶x导数 # 构建候选项矩阵 (n_f, n_terms) terms = torch.cat([ u_pred_f, # u u_x, # u_x u_xx, # u_xx u_xxx, # u_xxx u_pred_f * t_f, # u*t torch.sin(np.pi * x_f) # sin(x) ], dim=1) # 3. 物理损失：约束 u_t = Σ(系数 * 候选项) pred_u_t = torch.matmul(terms, term_coeffs) # 预测的u_t loss_physics = torch.mean((u_t - pred_u_t) ** 2) # 4. L1正则：促进系数稀疏（让冗余项系数趋近于0，识别关键项） l1_reg = 1e-2 * torch.sum(torch.abs(term_coeffs)) # 稀疏正则 # 总损失 loss = loss_data + 10 * loss_physics + l1_reg loss.backward() optimizer.step() # 每500轮打印进度 if (epoch + 1) % 500 == 0: print(f"Epoch {epoch+1:4d} | Loss: {loss.item():.6f} | Physics Loss: {loss_physics.item():.6f}") # ---------------------- 6. 结果解析：识别微分项+系数 ---------------------- # 提取系数并筛选关键项（阈值：系数绝对值>1e-4） coeffs_np = term_coeffs.detach().numpy() threshold = 1e-4 key_indices = np.where(np.abs(coeffs_np) > threshold)[0] key_terms = [term_names[i] for i in key_indices] key_coeffs = coeffs_np[key_indices] # 构建发现的微分方程 eq_terms = [f"{coeff[0]:.6f}*{term}" for coeff, term in zip(key_coeffs, key_terms)] discovered_eq = f"u_t = {' + '.join(eq_terms)}" if eq_terms else "u_t = 0" # 输出结果 print("\n======== 微分方程发现结果 ========") print(f"真实方程：u_t = {true_a_sq:.6f}*u_xx") print(f"发现方程：{discovered_eq}") print("\n所有候选项系数：") for name, coeff in zip(term_names, coeffs_np): print(f"{name:6s} : {coeff[0]:.6f}") # ---------------------- 7. 可视化预测结果 ---------------------- # 测试数据 x_test = torch.linspace(0, 1, 50).view(-1, 1) t_test = torch.linspace(0, 1, 50).view(-1, 1) X_test = torch.cartesian_prod(x_test.squeeze(), t_test.squeeze()).view(-1, 2) u_exact = exact_solution_heat(true_a_sq, X_test[:, 0:1], X_test[:, 1:2]) u_pred = model(X_test).detach() plt.figure(figsize=(10, 4)) plt.subplot(1, 2, 1) plt.title("PINN Prediction") plt.tricontourf(X_test[:, 0], X_test[:, 1], u_pred.squeeze(), cmap="jet", levels=30) plt.colorbar(label="u") plt.xlabel("x") plt.ylabel("t") plt.subplot(1, 2, 2) plt.title("Exact Solution") plt.tricontourf(X_test[:, 0], X_test[:, 1], u_exact.squeeze(), cmap="jet", levels=30) plt.colorbar(label="u") plt.xlabel("x") plt.ylabel("t") plt.tight_layout() plt.show()

这是对候选项的设定。

terms = torch.cat([ u_pred_f, # u u_x, # u_x u_xx, # u_xx u_xxx, # u_xxx u_pred_f * t_f, # u*t torch.sin(np.pi * x_f) # sin(x) ], dim=1)

下面是对候选项的系数的设定。

# 可学习系数：每个候选微分项对应一个系数（初始值随机） term_coeffs = nn.Parameter(torch.randn(n_terms, 1), requires_grad=True) optimizer = torch.optim.Adam(list(model.parameters()) + [term_coeffs], lr=1e-3)

接下来给出训练15000次最后的损失以及方程的发现结果。

Epoch 15000 | Loss: 0.000126 | Physics Loss: 0.000002

======== 微分方程发现结果 ========
真实方程：u_t = 0.010000*u_xx
发现方程：u_t = 0.008118*u_xx

所有候选项系数：
u : 0.000020
u_x : 0.000078
u_xx : 0.008118
u_xxx : 0.000071
u*t : -0.000076
sin(x) : 0.000005

这是我不断调整各超参数最后得到的一个较为合适的解（很多时候得到的方程包含各阶微分项），其实这样也暴露了简单PINN来发现方程的缺点：

1、各损失项的超参数的设定问题

2、对系数矩阵进行稀疏性约束的方法问题