4. Burgers方程粘性解
任务目标
利用物理信息神经网络(PINN)求解一维Burgers方程，在给定初始条件下捕捉激波形成过程，并对比不同粘性系数对解的影响。通过本案例掌握PINN处理非线性对流占优问题的基本流程。
问题定义
求解如下Burgers方程：
u_t + u*u_x = ν*u_xx, x∈[-1,1], t∈[0,1]初始条件：u(0,x) = -sin(πx)
边界条件：u(t,-1) = u(t,1) = 0
粘性系数ν分别取0.01（弱粘性，形成近似激波）和0.1（强粘性，解光滑化），对比两种情况下解的演化差异。
环境依赖
pip install tensorflow numpy matplotlib关键完整代码
import numpy as np
import tensorflow as tf
import matplotlib.pyplot as plt
# 设置随机种子保证可复现
tf.random.set_seed(42)
np.random.seed(42)
class BurgersPINN:
def __init__(self, nu):
self.nu = nu
# 神经网络：输入(x,t)，输出u
self.model = tf.keras.Sequential([
tf.keras.layers.Dense(20, activation='tanh', input_shape=(2,)),
tf.keras.layers.Dense(20, activation='tanh'),
tf.keras.layers.Dense(20, activation='tanh'),
tf.keras.layers.Dense(1)
])
self.optimizer = tf.keras.optimizers.Adam(learning_rate=1e-3)
def loss(self, x, t):
with tf.GradientTape() as tape2:
tape2.watch([x, t])
with tf.GradientTape() as tape1:
tape1.watch([x, t])
u = self.model(tf.concat([x, t], axis=1))
grads = tape1.gradient(u, [x, t])
u_x = grads[0]
u_t = grads[1]
u_xx = tape2.gradient(u_x, x)
# PDE残差: u_t + u*u_x - nu*u_xx = 0
pde_res = u_t + u * u_x - self.nu * u_xx
# 初始条件残差: u(0,x) + sin(πx) = 0
ic_res = self.model(tf.concat([tf.zeros_like(x), x], axis=1)) + tf.sin(np.pi * x)
# 边界条件残差: u(t,-1) = u(t,1) = 0
bc_left = self.model(tf.concat([t, tf.ones_like(t) * -1.0], axis=1))
bc_right = self.model(tf.concat([t, tf.ones_like(t) * 1.0], axis=1))
return tf.reduce_mean(tf.square(pde_res)) + \
tf.reduce_mean(tf.square(ic_res)) + \
tf.reduce_mean(tf.square(bc_left)) + \
tf.reduce_mean(tf.square(bc_right))
def train(self, epochs=5000):
# 生成配点
x = tf.random.uniform((2000, 1), -1.0, 1.0)
t = tf.random.uniform((2000, 1), 0.0, 1.0)
x_ic = tf.random.uniform((500, 1), -1.0, 1.0)
t_bc = tf.random.uniform((500, 1), 0.0, 1.0)
for epoch in range(epochs):
with tf.GradientTape() as tape:
loss_val = self.loss(x, t)
grads = tape.gradient(loss_val, self.model.trainable_variables)
self.optimizer.apply_gradients(zip(grads, self.model.trainable_variables))
if epoch % 1000 == 0:
print(f'Epoch {epoch}, Loss: {loss_val.numpy():.6f}')
def predict(self, x_star, t_star):
return self.model(tf.concat([t_star, x_star], axis=1)).numpy()
# 训练两个不同粘性系数的模型
print("=== 训练 ν=0.01 模型 ===")
model_low_nu = BurgersPINN(nu=0.01)
model_low_nu.train()
print("\n=== 训练 ν=0.1 模型 ===")
model_high_nu = BurgersPINN(nu=0.1)
model_high_nu.train()
# 在网格上预测并可视化
x_grid = np.linspace(-1, 1, 200)
t_grid = np.linspace(0, 1, 100)
X, T = np.meshgrid(x_grid, t_grid)
x_flat = X.flatten().reshape(-1, 1).astype(np.float32)
t_flat = T.flatten().reshape(-1, 1).astype(np.float32)
u_low = model_low_nu.predict(x_flat, t_flat).reshape(100, 200)
u_high = model_high_nu.predict(x_flat, t_flat).reshape(100, 200)
# 绘制不同时刻的解对比
fig, axes = plt.subplots(2, 3, figsize=(15, 8))
time_indices = [0, 25, 50] # t=0, 0.25, 0.5
for i, idx in enumerate(time_indices):
# ν=0.01
axes[0, i].plot(x_grid, u_low[idx], 'b-', label=f'ν=0.01, t={t_grid[idx]:.2f}')
axes[0, i].set_title(f'ν=0.01, t={t_grid[idx]:.2f}')
axes[0, i].set_xlim([-1, 1])
axes[0, i].set_ylim([-1.5, 1.5])
axes[0, i].grid(True)
# ν=0.1
axes[1, i].plot(x_grid, u_high[idx], 'r-', label=f'ν=0.1, t={t_grid[idx]:.2f}')
axes[1, i].set_title(f'ν=0.1, t={t_grid[idx]:.2f}')
axes[1, i].set_xlim([-1, 1])
axes[1, i].set_ylim([-1.5, 1.5])
axes[1, i].grid(True)
plt.tight_layout()
plt.savefig('burgers_comparison.png', dpi=150)
plt.show()
运行方式
python burgers_pinn.py验证方法与判定标准
训练过程中观察loss值应持续下降至1e-3以下。训练完成后，程序生成burgers_comparison.png对比图，验证标准如下：
- ν=0.01：在t≈0.25后，x=0附近应出现陡峭的梯度过渡带（近似激波），解从u≈+1快速过渡到u≈-1，过渡带宽度随t增大而略微变宽；
- ν=0.1：解保持光滑，无激波形成，梯度过渡带明显宽于ν=0.01情形；
- 两种情况下t=0时刻的解都应接近
-sin(πx)的初始波形。
常见问题
- loss不下降：尝试增大配点数量至5000以上，或增加训练轮数至10000；
- 激波处振荡：ν=0.01时激波附近可能出现Gibbs振荡，可增加网络宽度（如50个神经元）或使用分段激活函数缓解；
- 训练时间过长：本案例在CPU上约需5-10分钟，可适当减少配点数至1000加速验证流程。