📚 课程目录第4章波动方程初值问题

4. 波动方程初值问题

任务目标

使用物理信息神经网络(PINN)求解一维波动方程初值问题，模拟弦的自由振动，验证网络对双曲型偏微分方程的求解能力，并与解析解进行定量对比。

问题定义

求解如下一维波动方程初值问题：

∂²u/∂t² = c² · ∂²u/∂x²,  x ∈ [0, 1], t ∈ [0, 1]

初始条件：u(x, 0) = sin(πx)，∂u/∂t(x, 0) = 0
边界条件：u(0, t) = u(1, t) = 0
波速：c = 1

该问题的解析解为 u(x, t) = sin(πx)cos(πt)，用于验证PINN精度。

实现要点

PINN损失函数包含三部分：方程残差（内部配点满足波动方程）、初始位移残差、初始速度残差（对t求偏导为0）。边界条件通过强制构造满足：令网络输出乘以 x(1-x) 因子，自动满足两端零边界。网络输入为 (x, t)，输出为位移 u。

完整代码

import torch
import torch.nn as nn
import numpy as np

# ---------- 网络定义 ----------
class PINN(nn.Module):
    def __init__(self):
        super().__init__()
        self.net = nn.Sequential(
            nn.Linear(2, 64), nn.Tanh(),
            nn.Linear(64, 64), nn.Tanh(),
            nn.Linear(64, 64), nn.Tanh(),
            nn.Linear(64, 1)
        )
    def forward(self, x, t):
        # 强制边界条件 u(0,t)=u(1,t)=0
        u = self.net(torch.cat([x, t], dim=1))
        return u * x * (1 - x)

# ---------- 损失函数 ----------
def loss_fn(model, x_eq, t_eq, x_ic, t_ic, x_ic_v, t_ic_v):
    x_eq.requires_grad_(True)
    t_eq.requires_grad_(True)
    
    # 方程残差: u_tt - u_xx = 0
    u = model(x_eq, t_eq)
    u_t = torch.autograd.grad(u, t_eq, grad_outputs=torch.ones_like(u), create_graph=True)[0]
    u_tt = torch.autograd.grad(u_t, t_eq, grad_outputs=torch.ones_like(u_t), create_graph=True)[0]
    u_x = torch.autograd.grad(u, x_eq, grad_outputs=torch.ones_like(u), create_graph=True)[0]
    u_xx = torch.autograd.grad(u_x, x_eq, grad_outputs=torch.ones_like(u_x), create_graph=True)[0]
    loss_eq = torch.mean((u_tt - u_xx) ** 2)
    
    # 初始位移: u(x,0) = sin(πx)
    u_ic = model(x_ic, t_ic)
    loss_ic = torch.mean((u_ic - torch.sin(np.pi * x_ic)) ** 2)
    
    # 初始速度: u_t(x,0) = 0
    u_v = model(x_ic_v, t_ic_v)
    u_v_t = torch.autograd.grad(u_v, t_ic_v, grad_outputs=torch.ones_like(u_v), create_graph=True)[0]
    loss_iv = torch.mean(u_v_t ** 2)
    
    return loss_eq + loss_ic + loss_iv

# ---------- 训练数据 ----------
torch.manual_seed(42)
N_eq = 2000   # 内部配点
N_ic = 200    # 初始条件点

# 方程配点: 拉丁超立方采样
x_eq = torch.rand(N_eq, 1)
t_eq = torch.rand(N_eq, 1) * 1.0

# 初始位移点 t=0
x_ic = torch.rand(N_ic, 1)
t_ic = torch.zeros(N_ic, 1)

# 初始速度点 t=0 (需要二阶导，单独一组)
x_ic_v = torch.rand(N_ic, 1)
t_ic_v = torch.zeros(N_ic, 1)
t_ic_v.requires_grad_(True)

# ---------- 训练 ----------
model = PINN()
optimizer = torch.optim.Adam(model.parameters(), lr=1e-3)

for epoch in range(5000):
    optimizer.zero_grad()
    loss = loss_fn(model, x_eq, t_eq, x_ic, t_ic, x_ic_v, t_ic_v)
    loss.backward()
    optimizer.step()
    if epoch % 500 == 0:
        print(f'Epoch {epoch:5d}, Loss: {loss.item():.2e}')

# ---------- 验证 ----------
# 在网格上对比解析解
x_test = torch.linspace(0, 1, 101).unsqueeze(1)
t_test = torch.linspace(0, 1, 101).unsqueeze(1)
X, T = torch.meshgrid(x_test.squeeze(), t_test.squeeze(), indexing='ij')
Xf = X.reshape(-1, 1)
Tf = T.reshape(-1, 1)

u_pred = model(Xf, Tf).detach().reshape(101, 101)
u_exact = torch.sin(np.pi * X) * torch.cos(np.pi * T)

# 相对L2误差
err = torch.norm(u_pred - u_exact) / torch.norm(u_exact)
print(f'相对L2误差: {err.item():.4f}')

# 保存预测结果供可视化
np.savez('wave_pred.npz', x=x_test.numpy(), t=t_test.numpy(), 
         u_pred=u_pred.numpy(), u_exact=u_exact.numpy())

运行方式

需要安装 torch 和 numpy。直接运行脚本即可，训练约5000轮（CPU约2～3分钟）。程序会打印每500轮的损失值和最终相对L2误差。

验证方法与判定标准

运行结束后观察输出：

损失收敛：训练末期损失应降至1e-4以下（理想情况可达1e-5～1e-6）；
相对L2误差：与解析解 u(x,t)=sin(πx)cos(πt) 对比，误差小于0.01即认为PINN成功求解该波动方程；
物理合理性：可绘制 t=0, 0.25, 0.5 时刻的波形，应分别对应 sin(πx)、0.707·sin(πx)、0，体现驻波特征。

常见问题

初始速度条件难以收敛：若初始速度损失下降缓慢，可适当增大其权重（如乘以10），或增加 t=0 附近的配点密度。
长时间域发散：本案例 t∈[0,1] 较短，若扩展时间域需增加配点数量或采用分段训练。
边界强制技巧：乘以 x(1-x) 因子使网络输出在 x=0 和 x=1 处恒为零，避免了边界损失项的调参，但会略微增加训练难度，可配合适当增大网络宽度缓解。