用神经网络近似椭圆PDE时损失小但最大误差大的原因及优化

我尝试用神经网络近似椭圆偏微分方程的解，代码如下。训练后损失达到1e-4，但最大误差高达1.47，近似效果不佳。调整超参数后问题仍存在，请问原因是什么？如何降低最大误差？

#PDE
import tensorflow as tf
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D


class Net(tf.keras.Model):
    def __init__(self):
        super(Net, self).__init__()
        self.hidden1 = tf.keras.layers.Dense(100, activation='relu')
        self.hidden2 = tf.keras.layers.Dense(200, activation='relu')
        self.output_layer = tf.keras.layers.Dense(1)

    def call(self, x, y):
        inputs = tf.concat([x, y], axis=1)
        x = self.hidden1(inputs)
        x = self.hidden2(x)
        x = self.output_layer(x)
        return x

def f(x, y):
    return tf.exp(-x) * (x - 2 + y**3 + 6 * y)

def trial_func(x, y, net):
    A = (1 - x) * y**3 + x * (1 + y**3) * np.exp(-1) + (1 - y) * x * (tf.exp(-x) - np.exp(-1))     + y * [(1 + x) * tf.exp(-x) -(1 - x - 2 * x * np.exp(-1))]
    u = A + x * (1 - x) * y * (1 - y) * net(x, y)
    return u

def loss_func(d2udx2, d2udy2, x, y, net):
    u = trial_func(x, y, net)
    return tf.reduce_mean(tf.square(d2udx2 + d2udy2 - f(x, y)))

def exact_solution(x, y):
    return tf.exp(-x)*(x+y**3)

x_train = np.linspace(0, 1, 100, dtype=np.float32)
y_train = np.linspace(0, 1, 100, dtype=np.float32)
x_train, y_train = np.meshgrid(x_train, y_train)

x_train = tf.convert_to_tensor(x_train, dtype=tf.float32)
y_train = tf.convert_to_tensor(y_train, dtype=tf.float32)


u_train = tf.zeros_like(x_train)

model = Net()
optimizer = tf.keras.optimizers.Adam(learning_rate=0.01)

@tf.function
def compute_gradients(x, y, u):
    with tf.GradientTape(persistent=True) as tape:
        tape.watch([x, y])
        u_pred = trial_func(x, y, model)
        d2udx2_pred = tape.gradient(tape.gradient(u_pred, x), x)
        d2udy2_pred = tape.gradient(tape.gradient(u_pred, y), y)
        loss = loss_func(d2udx2_pred, d2udy2_pred, x, y, model)

    gradients = tape.gradient(loss, model.trainable_variables)
    optimizer.apply_gradients(zip(gradients, model.trainable_variables))
    return loss, gradients

for epoch in range(1000):
    loss, gradients = compute_gradients(x_train, y_train, u_train)
    
    if epoch % 100 == 0:
        print(f"Epoch: {epoch}, Loss: {loss.numpy()}")

# Reshape the input data for plotting
x_test = np.linspace(0, 1, 100)
y_test = np.linspace(0, 1, 100)
x_test, y_test = np.meshgrid(x_test, y_test)
x_test = tf.constant(x_test, dtype=tf.float32)
y_test = tf.constant(y_test, dtype=tf.float32)

# Evaluate the predicted values
u_pred = trial_func(x_test, y_test, model)

#-------------------------------------------------------------------------------

# Reshape the predicted values
u_pred = np.reshape(u_pred.numpy(), x_test.shape)

# Plotting the actual and predicted graph
fig = plt.figure(figsize=(10, 6))
ax = fig.add_subplot(111, projection='3d')

# Actual function values
ax.plot_surface(x_test, y_test, exact_solution(x_test, y_test), cmap='viridis', alpha=0.5, label='Actual')

# Predicted function values
ax.plot_surface(x_test, y_test, u_pred, cmap='plasma', alpha=0.8, label='Predicted')

ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('u')
ax.set_title('Actual vs. Predicted Solution')

# Show the plot
plt.show()

# Compute the error
error = np.abs(exact_solution(x_test, y_test) - u_pred)

# Maximum error
max_error = np.max(error)

# Plotting the error graph
fig = plt.figure(figsize=(8, 6))
ax = fig.add_subplot(111, projection='3d')

# Error values
ax.plot_surface(x_test, y_test, error, cmap='Reds', alpha=0.8, label='Error')

ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('Error')
ax.set_title('Error')

# Show the plot
plt.show()

print(f"Maximum error: {max_error}")

核心原因分析 #

1. 试函数边界条件不满足（最大误差的直接来源） #

你的最大误差正好对应y=1边界处的偏差，这是因为试函数A的构造错误，导致y=1时不满足PDE的边界条件：

• 精确解在y=1时为exact_solution(x,1) = exp(-x)(x+1)
• 但你的试函数A在y=1时计算结果为4x*exp(-1) + (1+x)*exp(-x)，比精确解多了4x*exp(-1)项，当x=1时这个额外项的值为4*exp(-1)≈1.475，与你观测到的最大误差完全一致。

2. 损失函数的局限性 #

当前使用的**均方误差（MSE）**是所有样本误差的平均值，它会优先优化大多数区域的拟合效果，但对少数极端误差（比如边界处的系统性偏差）不敏感。即使边界处误差极大，只要内部区域误差足够小，整体MSE仍能降到很低。

修复方案 #

1. 修正试函数的边界条件 #

重新构造试函数A，确保它在所有边界（x=0,x=1,y=0,y=1）上严格等于精确解。可以用双线性插值消除边界重叠冲突：

def trial_func(x, y, net):
    # 构造满足所有边界条件的A
    A = (1 - x) * y**3 \
        + x * (1 + y**3) * np.exp(-1) \
        + (1 - y) * x * tf.exp(-x) \
        + y * (x + 1) * tf.exp(-x) \
        - (1 - x)*(1 - y)*0 \
        - x*(1 - y)*np.exp(-1) \
        - (1 - x)*y*1 \
        - x*y*(1+1)*np.exp(-1)
    # 加入可训练的内部修正项
    u = A + x * (1 - x) * y * (1 - y) * net(x, y)
    return u

2. 优化损失函数与训练策略 #

• 加入边界误差惩罚：强制模型满足边界条件，给边界误差更高权重：

  def loss_func(d2udx2, d2udy2, x, y, net):
      # PDE残差损失
      pde_loss = tf.reduce_mean(tf.square(d2udx2 + d2udy2 - f(x, y)))
      # 边界解误差损失
      boundary_mask = (x == 0) | (x == 1) | (y == 0) | (y == 1)
      u_pred = trial_func(x, y, net)
      u_exact = exact_solution(x, y)
      boundary_loss = tf.reduce_mean(tf.square(tf.boolean_mask(u_pred, boundary_mask) - tf.boolean_mask(u_exact, boundary_mask)))
      # 组合损失
      return pde_loss + 10 * boundary_loss
  

• 改用鲁棒损失函数：用Huber损失替代MSE，降低极端误差的影响：

  def loss_func(d2udx2, d2udy2, x, y, net):
      residual = d2udx2 + d2udy2 - f(x, y)
      pde_loss = tf.reduce_mean(tf.keras.losses.huber(residual, tf.zeros_like(residual)))
      return pde_loss
  

• 调整训练参数：
  • 用学习率衰减：tf.keras.optimizers.schedules.ExponentialDecay(initial_learning_rate=0.01, decay_steps=1000, decay_rate=0.9)
  • 增加训练轮数至5000-10000轮
  • 优化网络结构：增加隐藏层数量（如3层）或神经元数，改用tanh激活函数避免边界处梯度消失

3. 验证梯度计算正确性 #

二阶导数的数值计算可能不稳定，建议用有限差分法验证GradientTape的结果：

# 有限差分计算x方向二阶导数
def finite_difference_second_x(u, x, y, model, h=1e-4):
    u_plus = trial_func(x+h, y, model)
    u_minus = trial_func(x-h, y, model)
    return (u_plus - 2*u + u_minus)/(h**2)

内容的提问来源于stack exchange，提问作者Tan