2. 二维拉普拉斯方程边值问题

任务目标

在矩形域 \([0,1]\times[0,1]\) 上求解二维拉普拉斯方程 \(\nabla^2 u = 0\)，边界条件为：

• \(u(0,y)=0,\quad u(1,y)=0\)
• \(u(x,0)=0,\quad u(x,1)=\sin(\pi x)\)

使用物理信息神经网络（PINN）逼近势场 \(u(x,y)\)，并将预测结果与有限差分法（FDM）数值解对比，验证PINN求解边值问题的有效性。

实现要点

PINN将偏微分方程残差与边界条件残差共同作为损失函数。网络输入 \((x,y)\)，输出 \(\hat{u}(x,y)\)。总损失为：

\[ \mathcal{L} = \lambda_{PDE}\cdot \frac{1}{N_f}\sum_{i=1}^{N_f} |\nabla^2 \hat{u}(x_i,y_i)|^2 + \lambda_{BC}\cdot \frac{1}{N_b}\sum_{j=1}^{N_b} |\hat{u}(x_j,y_j) - g(x_j,y_j)|^2 \]

其中 \(g\) 为Dirichlet边界值。使用自动微分计算二阶导数，Adam优化器训练。有限差分解作为参考基准，采用五点差分格式在相同网格上求解。

完整代码

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

# ---------- 1. PINN 网络定义 ----------
class PINN(nn.Module):
    def __init__(self):
        super().__init__()
        self.net = nn.Sequential(
            nn.Linear(2, 50), nn.Tanh(),
            nn.Linear(50, 50), nn.Tanh(),
            nn.Linear(50, 50), nn.Tanh(),
            nn.Linear(50, 1)
        )
    def forward(self, x, y):
        return self.net(torch.cat([x, y], dim=1))

# ---------- 2. 拉普拉斯算子（自动微分） ----------
def laplacian(model, x, y):
    x.requires_grad_(True)
    y.requires_grad_(True)
    u = model(x, y)
    u_x = torch.autograd.grad(u, x, grad_outputs=torch.ones_like(u), create_graph=True)[0]
    u_y = torch.autograd.grad(u, y, grad_outputs=torch.ones_like(u), create_graph=True)[0]
    u_xx = torch.autograd.grad(u_x, x, grad_outputs=torch.ones_like(u_x), create_graph=True)[0]
    u_yy = torch.autograd.grad(u_y, y, grad_outputs=torch.ones_like(u_y), create_graph=True)[0]
    return u_xx + u_yy

# ---------- 3. 训练数据生成 ----------
N_f = 2000   # 内部配点数
N_b = 400    # 边界配点数
# 内部点
x_f = torch.rand(N_f, 1) * 1.0
y_f = torch.rand(N_f, 1) * 1.0
# 边界点（四条边均匀采样）
t = torch.linspace(0, 1, N_b//4).view(-1, 1)
x_b = torch.cat([torch.zeros_like(t), torch.ones_like(t), t, t])
y_b = torch.cat([t, t, torch.zeros_like(t), torch.ones_like(t)])
# 边界真值
g_b = torch.zeros_like(x_b)
g_b[(x_b==1.0) & (y_b==1.0)] = 0  # 角点修正
# 上边界 y=1: u = sin(pi*x)
mask_top = (y_b == 1.0)
g_b[mask_top] = torch.sin(np.pi * x_b[mask_top])

# ---------- 4. 训练 ----------
model = PINN()
optimizer = torch.optim.Adam(model.parameters(), lr=1e-3)
loss_fn = nn.MSELoss()
lambda_pde, lambda_bc = 1.0, 1.0

for epoch in range(3000):
    optimizer.zero_grad()
    # PDE残差
    lap = laplacian(model, x_f, y_f)
    loss_pde = loss_fn(lap, torch.zeros_like(lap))
    # 边界残差
    u_b = model(x_b, y_b)
    loss_bc = loss_fn(u_b, g_b)
    # 总损失
    loss = lambda_pde * loss_pde + lambda_bc * loss_bc
    loss.backward()
    optimizer.step()
    if epoch % 500 == 0:
        print(f'Epoch {epoch:4d} | PDE Loss: {loss_pde.item():.2e} | BC Loss: {loss_bc.item():.2e}')

# ---------- 5. 有限差分参考解 ----------
N = 101
x = np.linspace(0, 1, N)
y = np.linspace(0, 1, N)
dx = dy = 1/(N-1)
U = np.zeros((N, N))
U[-1, :] = np.sin(np.pi * x)  # y=1 边界
# 迭代求解（Gauss-Seidel）
for _ in range(10000):
    U_old = U.copy()
    for i in range(1, N-1):
        for j in range(1, N-1):
            U[i, j] = 0.25 * (U[i+1, j] + U[i-1, j] + U[i, j+1] + U[i, j-1])
    if np.max(np.abs(U - U_old)) < 1e-6:
        break

# ---------- 6. 对比与可视化 ----------
# PINN预测
x_grid, y_grid = np.meshgrid(x, y)
x_flat = torch.tensor(x_grid.flatten(), dtype=torch.float32).view(-1, 1)
y_flat = torch.tensor(y_grid.flatten(), dtype=torch.float32).view(-1, 1)
with torch.no_grad():
    u_pinn = model(x_flat, y_flat).numpy().reshape(N, N)

# 误差统计
err = np.abs(u_pinn - U)
print(f'\n最大误差: {err.max():.4e}')
print(f'平均误差: {err.mean():.4e}')

# 绘制对比
fig, axes = plt.subplots(1, 3, figsize=(15, 4))
im0 = axes[0].contourf(x_grid, y_grid, U, levels=50, cmap='viridis')
axes[0].set_title('FDM Reference')
plt.colorbar(im0, ax=axes[0])
im1 = axes[1].contourf(x_grid, y_grid, u_pinn, levels=50, cmap='viridis')
axes[1].set_title('PINN Prediction')
plt.colorbar(im1, ax=axes[1])
im2 = axes[2].contourf(x_grid, y_grid, err, levels=50, cmap='Reds')
axes[2].set_title('Absolute Error')
plt.colorbar(im2, ax=axes[2])
plt.tight_layout()
plt.savefig('pinn_vs_fdm.png', dpi=150)
plt.show()

运行与验证

1. 依赖安装：需要 torch、numpy、matplotlib。执行 pip install torch numpy matplotlib。
2. 运行方式：将代码保存为 pinn_laplace2d.py，执行 python pinn_laplace2d.py。
3. 通过条件：训练结束后终端输出最大误差应小于 \(10^{-2}\) 量级（通常约 \(10^{-3}\) 左右）；生成的 pinn_vs_fdm.png 中，PINN与FDM云图分布一致，误差图颜色整体偏淡（误差值小）。

常见问题

• 训练不收敛：增大训练轮数至5000以上，或调整学习率为 \(5\times10^{-4}\)。
• 边界角点误差大：角点处Dirichlet条件不连续（\(u(1,1)=0\) 与 \(u(1,1)=\sin(\pi)=0\) 一致，但 \(u(0,1)=0\) 与 \(\sin(0)=0\) 也一致，本案例角点无冲突）。若遇不一致边界，需对角点单独加权处理。
• FDM迭代慢：可改用向量化更新或使用 scipy.sparse.linalg.spsolve 直接求解线性方程组。