基于物理信息神经网络(PINN)求解三维声波方程的MATLAB实战指南-开发者社区

简介：本资源是一套基于物理信息神经网络（PINN）求解三维声波波动方程的MATLAB实现方案，面向计算物理、声学仿真与AI驱动科学计算领域的研究生及科研工程师，解决传统数值方法在高维复杂边界下建模成本高、泛化性弱的问题。压缩包共3个文件（2个核心MATLAB脚本+1个演示动画），总大小466KB，其中main.m负责问题定义、网络构建与训练主流程，modelLoss.m封装PDE残差、初始/边界条件联合损失函数，MP4动画直观展示三维波场随时间演化的传播过程。已有204人学习下载，代码结构高度模块化，完整覆盖数据生成、网络设计、损失构造、Adam优化训练及多维度可视化（损失曲线、空间切片、动态传播），无需额外依赖库，开箱即用，适合快速复现PINN在偏微分方程求解中的前沿应用。

1. 项目概述：当神经网络学会“物理定律”

如果你正在处理三维声波传播、地震波模拟或者复杂介质中的波动问题，并且对传统数值方法（如有限差分、有限元）的网格划分、稳定性条件和计算成本感到头疼，那么你找对地方了。今天要聊的，是一个将前沿人工智能技术与经典物理方程求解相结合的实战项目：基于物理信息神经网络（PINN）求解三维声波波动方程。

简单来说，PINN 是一种“聪明”的神经网络。它不像普通的深度学习模型那样，需要海量的、已经标注好的“输入-输出”数据对来学习。相反，它通过学习物理定律本身——在这里就是三维声波波动方程——来直接求解问题。你只需要告诉网络方程是什么、边界条件在哪里、初始状态如何，它就能自己“推导”出整个时空域内的声波场分布。这尤其适合那些数据稀缺、但物理规律明确的科学计算场景。

这个项目的核心价值在于，它提供了一套完整的、可复现的 MATLAB 代码框架。你拿到手的不再是零散的概念代码，而是一个从理论到实现、从数据生成到结果可视化的完整工具包。无论你是计算物理、地球物理勘探、无损检测还是声学设计领域的研究者或工程师，这套代码都能帮你快速搭建起自己的 PINN 求解器，绕过繁琐的网格生成和迭代计算，用一种更“优雅”的方式洞察波动现象。

2. 核心思路：PINN如何“捆绑”物理方程？

在深入代码之前，我们必须搞清楚 PINN 的“灵魂”所在。它之所以能求解偏微分方程，核心在于将物理定律作为“软约束”直接嵌入到神经网络的训练过程中。这与我们熟知的“数据驱动”模式有本质区别。

2.1 从数据驱动到物理驱动

传统的监督学习神经网络，可以看作一个强大的函数逼近器。我们给它看很多组(x, y, z, t) -> p的样本（即空间坐标和时间对应声压值），它学习一个复杂的映射关系f_theta(x, y, z, t) ≈ p。这里的theta是网络权重。这种方法严重依赖高质量、高密度的训练数据，而对于一个三维时空问题，获取这样的数据成本极高，甚至不现实。

PINN 换了一种思路。我们不再提供p的“标准答案”，而是提供一个衡量答案是否“合理”的标尺——物理方程。对于三维声波波动方程：

∂²p/∂t² = c² (∂²p/∂x² + ∂²p/∂y² + ∂²p/∂z²)

其中p(x, y, z, t)是声压，c(x, y, z)是声速。PINN 构建一个神经网络N(x, y, z, t; theta)来直接预测p。那么，如何确保N的输出满足上面的方程呢？

我们让神经网络自动计算其输出对输入（x, y, z, t）的偏导数。利用自动微分技术，可以精确地得到∂²N/∂t²、∂²N/∂x²等项。然后，我们构造一个物理残差（Physics Residual）：

r_pde = ∂²N/∂t² - c² (∂²N/∂x² + ∂²N/∂y² + ∂²N/∂z²)

如果N是波动方程的精确解，那么在任何点(x, y, z, t)上，r_pde都应该等于零。因此，我们的训练目标就从“拟合数据”变成了“最小化物理残差”。

2.2 损失函数的精心设计

单一的物理残差最小化还不足以定解，我们必须引入边界条件和初始条件。这是 PINN 设计中的关键一步，也是与实际问题对接的桥梁。

假设我们的求解域是一个立方体，时间从0到T。我们需要定义以下损失项：

PDE损失（L_pde）：在求解域内部随机采样一大批“配置点（Collocation Points）”，计算这些点上物理残差r_pde的均方误差（MSE）。这迫使网络学习满足控制方程。
边界条件损失（L_bc）：在立方体的六个边界面上采样点，计算网络预测值N与指定的边界条件（如 Dirichlet 边界：p=0；或 Neumann 边界：∂p/∂n=0）之间的 MSE。
初始条件损失（L_ic）：在时间t=0的整个空间截面上采样点，计算网络预测值N与给定的初始声压分布p0(x,y,z)及其时间导数（如果需要）之间的 MSE。

最终的总损失函数是这些项的加权和：

L_total = λ_pde * L_pde + λ_bc * L_bc + λ_ic * L_ic

这里的λ是权重系数，它们的平衡至关重要。如果λ_ic太小，网络可能无法捕捉正确的初始状态；如果λ_pde太大，可能会为了满足方程而轻微牺牲边界条件。在实际操作中，常常采用自适应权重策略，或者简单地将它们都设为1，通过调整采样点的数量来间接控制各损失的贡献度。

实操心得：损失平衡是门艺术初期训练时，经常发现边界或初始条件拟合得很好，但内部点的物理残差仍然很大。这是因为PDE损失相对于边界/初始条件损失是一个更高阶的约束（涉及二阶导数），优化起来更困难。一个有效的技巧是，在训练初期，适当增大λ_ic和λ_bc（例如设为10或100），让网络先快速“记住”边界和初始状态。训练几百轮后，再逐渐将它们调回1，让网络专注于降低内部物理残差。这类似于一个“课程学习”的过程。

3. MATLAB实现全流程拆解

下面，我们进入实战环节，一步步拆解这个三维PINN求解器的MATLAB实现。我将按照代码模块的顺序，解释每个部分的设计意图和关键细节。

3.1 环境与问题定义

首先，我们需要明确求解的问题域和物理参数。假设我们在一个[0, Lx] x [0, Ly] x [0, Lz]的三维空间和[0, T]的时间域内求解。声速c可以是常数，也可以是空间的函数（例如模拟分层介质）。

% 定义计算域和参数 Lx = 1.0; Ly = 1.0; Lz = 1.0; % 空间范围 T = 1.0; % 时间范围 c = 1.0; % 声速 (假设为常数) % 定义初始条件函数 (例如，一个高斯脉冲) p0 = @(x,y,z) exp(-100*((x-0.5).^2 + (y-0.5).^2 + (z-0.5).^2)); % 定义初始时间导数 (假设初始静止) dp0_dt = @(x,y,z) zeros(size(x)); % 定义边界条件 (假设为全吸收边界，简化处理为Dirichlet零边界) bc_func = @(x,y,z,t) zeros(size(x));

这里，初始条件我们设置了一个位于模型中心的高斯声压脉冲。边界条件简化为零，这对应着完全吸收边界（在无穷远处声压为零）。对于更复杂的吸收边界条件（如PML），需要在PDE残差或边界损失中引入更复杂的表达式。

3.2 神经网络架构搭建

在MATLAB中，我们可以使用Deep Learning Toolbox来构建神经网络。对于PINN，一个全连接的前馈网络（多层感知机，MLP）通常就足够了。关键是如何处理输入（x,y,z,t）和输出(p)。

function net = createPINN(layers) % layers: 一个数组，例如 [4, 20, 20, 20, 20, 1]，表示输入层4维，4个隐藏层每层20个神经元，输出层1维。 net = feedforwardnet(layers(2:end-1)); net = configure(net, zeros(layers(1),1), zeros(layers(end),1)); % 配置输入输出尺寸 % 非常重要：修改隐藏层的激活函数。默认的‘tansig’在二阶导数时表现可能不佳。 for i = 1:length(net.layers)-1 net.layers{i}.transferFcn = 'tanh'; % 或 ‘sin’， ‘swish’等。tanh是PINN中常见且稳定的选择。 end % 输出层通常使用线性激活函数，因为声压值范围没有限制。 net.layers{end}.transferFcn = 'purelin'; % 初始化权重。Xavier/Glorot初始化对深度网络更友好。 net = init(net); end

注意事项：激活函数的选择tanh函数是PINN中的“常青树”，因为它处处光滑可微，且导数有界，这对自动微分计算高阶导数（如波动方程中的二阶导数）非常友好。避免使用ReLU这类分段线性函数，因为它在零点不可微，会导致二阶导数为零，无法正确传播波动方程的物理信息。近年来，sin(SIREN网络) 或自适应激活函数也被证明能提升PINN的性能，但对于入门，tanh是最稳妥的选择。

3.3 数据采样策略

采样点的分布直接影响训练效率和最终精度。我们需要为PDE损失、边界损失和初始损失分别采样。

% 1. 内部配置点采样 (用于PDE损失) N_pde = 50000; % 内部点数量 x_pde = Lx * rand(N_pde, 1); y_pde = Ly * rand(N_pde, 1); z_pde = Lz * rand(N_pde, 1); t_pde = T * rand(N_pde, 1); points_pde = [x_pde, y_pde, z_pde, t_pde]; % 2. 边界点采样 (每个面采样) N_bc_per_face = 2000; bc_points = []; faces = {'x0', 'xL', 'y0', 'yL', 'z0', 'zL'}; % 对应六个面 for face = faces [x_bc, y_bc, z_bc] = sampleOnFace(face{1}, Lx, Ly, Lz, N_bc_per_face); t_bc = T * rand(N_bc_per_face, 1); bc_points = [bc_points; [x_bc, y_bc, z_bc, t_bc]]; end % 3. 初始点采样 (t=0时刻的整个空间) N_ic = 10000; x_ic = Lx * rand(N_ic, 1); y_ic = Ly * rand(N_ic, 1); z_ic = Lz * rand(N_ic, 1); t_ic = zeros(N_ic, 1); % 时间固定为0 points_ic = [x_ic, y_ic, z_ic, t_ic];

sampleOnFace是一个辅助函数，用于在立方体的特定面上生成均匀随机点。边界点和初始点的数量可以少于内部点，因为它们约束的是低维流形（边界是三维时空中的二维面，初始时刻是三维空间）。一个常见的比例是 PDE点 : BC点 : IC点 = 10 : 1 : 1。

3.4 自动微分与损失计算

这是PINN实现中最核心、也最容易出错的部分。我们需要计算网络输出对输入的二阶偏导。MATLAB的dlgradient和dlfeval函数是实现自动微分的利器。

function [loss, gradients] = computeLoss(net, params, points_pde, points_bc, points_ic) % params: 包含网络权重、损失权重等的结构体 % 解构参数 weights = params.weights; lambda_pde = params.lambda_pde; lambda_bc = params.lambda_bc; lambda_ic = params.lambda_ic; c = params.c; % 将数据转换为 dlarray，以启用自动微分跟踪 dl_pde = dlarray(points_pde', 'CB'); % 格式： [特征维， 批大小] dl_bc = dlarray(points_bc', 'CB'); dl_ic = dlarray(points_ic', 'CB'); % 计算PDE损失 [loss_pde, ~] = dlfeval(@pdeResidual, net, weights, dl_pde, c); loss_pde = mean(loss_pde.^2); % 计算边界损失 pred_bc = forward(net, weights, dl_bc); true_bc = dlarray(zeros(1, size(points_bc,1)), 'CB'); % 假设零边界 loss_bc = mean((pred_bc - true_bc).^2); % 计算初始损失 (包括声压和声压时间导数) [loss_ic, loss_ic_dt] = dlfeval(@icResidual, net, weights, dl_ic); loss_ic_total = mean(loss_ic.^2) + mean(loss_ic_dt.^2); % 总损失 loss = lambda_pde * loss_pde + lambda_bc * loss_bc + lambda_ic * loss_ic_total; % 计算梯度（用于优化器） gradients = dlgradient(loss, weights); end function [residual, p] = pdeResidual(net, weights, dl_input, c) % 计算PDE残差： r = p_tt - c^2*(p_xx + p_yy + p_zz) [p, dp_dx, dp_dy, dp_dz, dp_dt] = computeDerivatives(net, weights, dl_input); % 计算二阶导数。注意：dlgradient需要标量输出，因此我们分别计算。 d2p_dx2 = dlgradient(sum(dp_dx, 'all'), dl_input, 'RetainData', true); d2p_dx2 = d2p_dx2(1,:); % 取对x的梯度 d2p_dy2 = dlgradient(sum(dp_dy, 'all'), dl_input, 'RetainData', true); d2p_dy2 = d2p_dy2(2,:); d2p_dz2 = dlgradient(sum(dp_dz, 'all'), dl_input, 'RetainData', true); d2p_dz2 = d2p_dz2(3,:); d2p_dt2 = dlgradient(sum(dp_dt, 'all'), dl_input); d2p_dt2 = d2p_dt2(4,:); laplacian_p = d2p_dx2 + d2p_dy2 + d2p_dz2; residual = d2p_dt2 - c^2 * laplacian_p; end function [p, dp_dx, dp_dy, dp_dz, dp_dt] = computeDerivatives(net, weights, dl_input) % 前向传播并计算一阶导数 p = forward(net, weights, dl_input); % 计算一阶导数 dp_dx = dlgradient(sum(p, 'all'), dl_input, 'RetainData', true); dp_dx = dp_dx(1,:); % 对x的偏导 dp_dy = dlgradient(sum(p, 'all'), dl_input, 'RetainData', true); dp_dy = dp_dy(2,:); % 对y的偏导 dp_dz = dlgradient(sum(p, 'all'), dl_input, 'RetainData', true); dp_dz = dp_dz(3,:); % 对z的偏导 dp_dt = dlgradient(sum(p, 'all'), dl_input); dp_dt = dp_dt(4,:); % 对t的偏导 end

这段代码有几个关键点：

dlarray与dlgradient：这是MATLAB进行自动微分的基础。dlarray封装数据并记录计算图。dlgradient计算梯度，‘RetainData’参数在计算高阶导时必须小心使用，以保留中间计算图。
高阶导数计算：计算波动方程的二阶导数需要调用两次dlgradient。第一次计算一阶导dp_dx，然后对dp_dx再求关于x的梯度，得到d2p_dx2。dlgradient(sum(dp_dx, ‘all’), dl_input, …)中的sum(…, ‘all’)是为了得到一个标量输出，这是dlgradient的要求。
计算图管理：频繁的高阶微分会导致计算图膨胀，消耗大量内存。代码中通过有选择地使用‘RetainData’来平衡。对于不用于后续梯度的中间变量，应避免保留其数据。

3.5 训练循环与优化器配置

有了损失函数，我们就可以用优化器来训练网络了。ADAM优化器是PINN训练的首选，因为它能自适应调整学习率，对非凸的损失地形有较好的适应性。

% 创建网络和参数 layers = [4, 50, 50, 50, 50, 1]; % 4输入，4层隐藏层每层50神经元，1输出 net = createPINN(layers); weights = getwb(net); % 获取初始权重向量 weights = dlarray(weights); % 转换为 dlarray % 训练参数 numEpochs = 20000; learningRate = 1e-3; lambda_pde = 1.0; lambda_bc = 1.0; lambda_ic = 10.0; % 初始条件权重稍大 % 使用ADAM优化器 averageGrad = []; averageSqGrad = []; decayRate = 0.9; % 一阶矩衰减率 sqDecayRate = 0.999; % 二阶矩衰减率 epsilon = 1e-8; lossHistory = []; for epoch = 1:numEpochs % 前向传播并计算损失和梯度 [loss, gradients] = computeLoss(net, struct('weights', weights, 'lambda_pde', lambda_pde, 'lambda_bc', lambda_bc, 'lambda_ic', lambda_ic, 'c', c), points_pde, points_bc, points_ic); % 记录损失 lossHistory = [lossHistory, extractdata(loss)]; % ADAM更新规则 [weights, averageGrad, averageSqGrad] = adamupdate(weights, gradients, averageGrad, averageSqGrad, epoch, learningRate, decayRate, sqDecayRate, epsilon); % 每1000轮打印一次损失 if mod(epoch, 1000) == 0 fprintf('Epoch %d, Loss: %.6e\n', epoch, extractdata(loss)); % 可以在这里动态调整学习率或损失权重 % if epoch > 5000 % lambda_ic = 1.0; % 降低初始条件权重 % end end end % 将训练好的权重写回网络对象 net = setwb(net, extractdata(weights));

训练PINN需要耐心。损失曲线通常不会像监督学习那样平滑下降，可能会在某个阶段停滞很长时间，然后突然下降。20000轮迭代是一个常见的起点，对于复杂问题可能需要更多。

3.6 结果验证与可视化

训练完成后，我们需要验证网络是否真的学会了波动方程的物理规律。最直接的方法是在一组规则网格点上进行预测，并与参考解（如果有的话，如有限差分法的解）进行对比，或者直观地检查波动的物理合理性。

% 在规则网格上生成测试点 [x_grid, y_grid, z_grid, t_grid] = ndgrid(linspace(0, Lx, 30), linspace(0, Ly, 30), linspace(0, Lz, 30), linspace(0, T, 10)); x_test = x_grid(:); y_test = y_grid(:); z_test = z_grid(:); t_test = t_grid(:); points_test = [x_test, y_test, z_test, t_test]; % 使用训练好的网络进行预测 p_pred = sim(net, points_test'); % 注意输入转置为 [特征维， 样本数] % 重塑预测结果以便可视化 p_pred_4d = reshape(p_pred, size(x_grid)); % 可视化：选择一个时间切片和深度切片 time_idx = 5; % 查看第5个时间步 depth_idx = 15; % 查看z方向的中间层 slice_xy = squeeze(p_pred_4d(:,:,depth_idx, time_idx)); figure; imagesc(linspace(0,Lx,30), linspace(0,Ly,30), slice_xy'); xlabel('X'); ylabel('Y'); title(sprintf('声压分布 (Z=%.2f, T=%.2f)', z_grid(1,1,depth_idx,1), t_grid(1,1,1,time_idx))); colorbar; axis image; colormap jet; % 可视化：空间某一点的声压随时间变化 point_x = 0.5; point_y = 0.5; point_z = 0.5; % 找到最近的网格索引 [~, idx_x] = min(abs(linspace(0,Lx,30)-point_x)); [~, idx_y] = min(abs(linspace(0,Ly,30)-point_y)); [~, idx_z] = min(abs(linspace(0,Lz,30)-point_z)); p_time_series = squeeze(p_pred_4d(idx_x, idx_y, idx_z, :)); t_series = squeeze(t_grid(1,1,1,:)); figure; plot(t_series, p_time_series, 'b-o', 'LineWidth', 1.5); xlabel('时间 (t)'); ylabel('声压 (p)'); title('中心点声压随时间变化'); grid on;

通过二维切片图，我们可以观察波前的传播是否呈圆形（对于均匀介质），以及边界反射是否被抑制（取决于边界条件设置）。通过时间序列图，可以检查波动是否符合预期的频率和衰减特性。这些都是定性判断解是否物理合理的重要手段。

4. 性能调优与高级技巧

基础框架跑通后，你可能会遇到精度不足、训练缓慢或难以收敛的问题。以下是一些经过实战检验的调优技巧。

4.1 网络架构与初始化

深度与宽度：对于三维波动方程，4-8个隐藏层，每层50-200个神经元是一个合理的起点。太小的网络容量不足，太大的网络则更难训练且容易过拟合。一个经验法则是，增加宽度比增加深度更能有效提升PINN的表达能力。
权重初始化：使用initnw(Nguyen-Widrow) 或init(默认) 初始化都可以。有研究表明，使用sin激活函数时，特定的初始化策略（如SIREN论文中的方法）能带来显著提升。对于tanh，保持默认初始化通常可行。
输入归一化：将输入坐标(x,y,z,t)归一化到[-1, 1]或[0, 1]区间，可以加速训练并提高稳定性。尤其是时间t和空间坐标尺度差异较大时，归一化至关重要。

4.2 损失函数与采样策略进阶

自适应权重：手动调整λ_pde，λ_bc，λ_ic很繁琐。可以采用“学习权重”的方法，将这些权重也作为可训练参数，让网络在训练中自动平衡各项损失。或者，使用基于损失值大小动态调整权重的算法，如 “Learning Rate Annealing for Physics-Informed Neural Networks” 中提出的方法。
残差自适应采样（RAR）：初始的随机采样可能无法捕捉到解变化剧烈的区域（如波前）。RAR策略在训练过程中，定期在物理残差r_pde较大的区域额外增加采样点，从而更高效地分配计算资源。实现方法是：每隔一定轮次，用当前网络评估一批新随机点的残差，选择残差最大的前N%的点，加入到训练点集中。
小批量训练：对于超大规模采样点（>10^6），一次性计算所有点的损失和梯度可能导致内存溢出。可以采用小批量随机梯度下降（SGD）或其变体。每次迭代从PDE点、BC点、IC点中分别随机抽取一个小批次进行计算。这引入了噪声，但能处理更大规模的问题。

4.3 针对波动方程的特殊处理

时间域分解：求解长时间演化问题时，单个PINN可能难以捕捉所有细节。可以采用“时间分段”策略，用多个PINN分别负责不同时间窗口[T_i, T_{i+1}]的求解，并在时间接口处施加连续性条件作为额外的损失项。
频域求解：对于时间谐波问题（即假设解为p(x,y,z,t)=P(x,y,z)*exp(iωt)），可以将波动方程转化为亥姆霍兹方程（Helmholtz Equation），从而消去时间变量，简化成三维空间问题。PINN同样可以求解亥姆霍兹方程。
硬边界条件：对于简单的Dirichlet零边界，除了通过损失函数约束，还可以通过构造特殊的网络结构来“硬编码”满足边界条件。例如，对于边界x=0和x=Lx上p=0，可以将网络输出设计为N(x,y,z,t) = x*(Lx-x)*N_raw(x,y,z,t)，这样无论N_raw输出什么，在边界上N自动为零。这能显著降低优化难度。

5. 常见问题与排查实录

在实际运行代码时，你几乎一定会遇到下面这些问题。这里是我的排查笔记。

5.1 损失不下降或震荡剧烈

症状：训练了几千轮，总损失居高不下，或者在某个值附近剧烈震荡。
排查清单：
1. 学习率过大：这是最常见的原因。将学习率从1e-3降到1e-4或5e-5试试。ADAM优化器对学习率相对鲁棒，但过大依然会导致震荡。
2. 损失权重失衡：检查L_pde，L_bc，L_ic各个分量的值。如果某一个比其他大几个数量级（例如L_pde是1e-3，而L_ic是1e+1），那么优化器会主要优化大的那一项。尝试调整λ权重，使各项损失在训练初期处于同一数量级（例如都在1e-1到1e+1之间）。
3. 网络表达能力不足：尝试增加网络层数或每层神经元数量。同时，检查激活函数是否正确（输出层应为‘purelin’）。
4. 梯度爆炸/消失：检查梯度值gradients。如果出现NaN或Inf，可能是计算图太深或激活函数选择不当。使用tanh而非sigmoid可以缓解梯度消失。也可以尝试梯度裁剪。
5. 采样点不足或分布不合理：增加N_pde的数量。确保边界点和初始点覆盖了所有边界和整个初始空间。

5.2 预测结果完全错误（如全零或常数）

症状：训练损失看起来下降了，但用网络预测任何输入，输出都是一个接近常数的值，没有波动现象。
排查清单：
1. 初始条件损失权重过低：网络可能找到了一个简单的解——满足波动方程和边界条件的常数解（例如p=0）。为了“打破对称性”，必须确保初始条件损失L_ic有足够大的权重，迫使网络在t=0时匹配非平凡的初始状态。这就是为什么我在示例代码中设置lambda_ic = 10.0。
2. 自动微分错误：这是最隐蔽的bug。务必单独验证自动微分计算的正确性。选择一个简单的测试函数，如p = sin(x)*cos(y)*exp(z)*t^2，手动计算其p_tt和laplacian(p)，然后与你pdeResidual函数计算的结果对比。在computeDerivatives和pdeResidual函数中设置断点，逐步检查每个一阶和二阶导数的值。
3. 输入数据未归一化：如果x,y,z的范围是[0,1]，而t的范围是[0, 100]，巨大的尺度差异会导致优化陷入病态。务必将所有输入特征归一化到相近的范围。

5.3 训练速度太慢

症状：每一轮迭代都耗时极长，尤其是当采样点很多时。
排查清单：
1. 向量化操作：确保你的代码充分利用了MATLAB的矩阵运算。避免在循环内对单个点调用dlfeval。computeLoss函数应一次性处理所有采样点（或一个批次）。
2. 减少高阶导数计算：波动方程需要二阶导数，计算开销大。在保证精度的前提下，可以尝试减少N_pde。或者，先使用较少的点训练一个粗糙解，再逐步增加点进行精细化训练（一种简单的延续法）。
3. 使用GPU：如果安装了Parallel Computing Toolbox且有兼容的GPU，将数据和网络转换为gpuArray可以带来数十倍的加速。将dlarray创建在GPU上：dl_pde = dlarray(gpuArray(points_pde’), ‘CB’);。
4. 检查计算图保留：不必要的‘RetainData’, true会积累计算图，占用大量内存并减慢速度。确保只在计算链式导数时保留必要的数据。

5.4 结果有物理不合理现象（如负能量、发散）

症状：波在传播过程中能量异常增长（发散）或出现非物理的振荡。
排查清单：
1. PDE残差未充分最小化：损失函数中的L_pde可能仍然较大。这意味着网络输出并没有很好地满足波动方程。继续训练，或检查PDE残差项的计算是否正确。
2. 边界条件处理不当：如果使用的是简单的零Dirichlet边界，在有限计算域内，波会在边界发生强反射，与“开放边界”的物理假设不符。你看到的发散可能是反射波的叠加。对于模拟无限域，需要引入吸收边界条件。在PINN中，这可以通过在边界损失项中使用特殊的吸收边界条件公式来实现（如 Clayton-Engquist 边界条件），或者更简单地，在物理域外设置一个“海绵层”，并在该层内增加阻尼项到波动方程中。
3. 数值色散：虽然PINN是基于网格无关的方法，但神经网络本身作为一种近似，在解高频分量时也可能引入误差。尝试增加网络容量，或在训练数据中增加对高频区域的采样（RAR策略对此有帮助）。

最后，我想分享一个最深的体会：PINN不是一个“即插即用”的黑箱。它更像一个需要精心调试的物理模拟器。成功的诀窍不在于追求最复杂的网络结构，而在于对物理问题的深刻理解（如何设置边界和初始条件）、对优化过程的细致观察（分析各项损失的变化）以及耐心的迭代调试。当你看到神经网络输出的波场动画与物理直觉完美吻合时，那种成就感是传统编程方法难以比拟的。这套MATLAB代码为你提供了一个坚实的起点，剩下的探索和创新，就交给你了。

本文还有配套的精品资源，点击获取