Matlab代码深度学习模拟Helmholtz教程：硬约束边界条件PINNs物理信息神经网络



Matlab代码深度学习模拟Helmholtz教程：硬约束边界条件PINNs物理信息神经网络

数据驱动与力学

2025-06-24

1. 数学背景与原理

1.1 亥姆霍兹方程

亥姆霍兹方程是描述波动现象的重要偏微分方程，一维形式为：

d²u/dx² + k²u = 0

其中：

u(x) 是待求解的声压场函数
k = 2πf/c₀ 是波数
f 是频率 (Hz)
c₀ 是介质中的声速 (m/s)

1.2 边界条件

本程序考虑狄利克雷边界条件：

u(0) = 1 Pa    (左边界)
u(1) = -1 Pa   (右边界)

1.3 解析解

对于给定边界条件的一维亥姆霍兹方程，解析解为：

u(x) = cos(kx) - (csc(k) + cot(k))sin(kx)

其中：

csc(k) = 1/sin(k) (余割函数)
cot(k) = cos(k)/sin(k) (余切函数)

2. 物理信息神经网络(PINNs)方法

2.1 基本思想

PINNs将物理定律(微分方程)作为约束条件嵌入到神经网络的损失函数中，使得网络在学习数据的同时自动满足物理定律。

2.2 试探函数

为了自动满足边界条件，构造试探函数：

G(x) = φ₁(x)·u₀₁ + φ₂(x)·u₀₂ + φₑ(x)·N(x)

其中：

φ₁(x) = 1-x (左边界形状函数)
φ₂(x) = x (右边界形状函数)
φₑ(x) = φ₁(x)·φ₂(x) = x(1-x) (等效形状函数)
u₀₁ = 1, u₀₂ = -1 (边界条件值)
N(x) 是神经网络输出

2.3 损失函数

损失函数基于微分方程残差：

Loss = ||d²G/dx² + k²G||²

2.4 硬约束边界条件

根据代码分析，MATLAB程序使用了一种非常巧妙的试探解方法（Trial Solution Method）来实现边界条件，这是PINNs中的一种先进技术。下面详细解释这种实现方式：

1. 边界条件定义

x0BC1 = 0;                      % 左边界 (m)
x0BC2 = 1;                      % 右边界 (m)

u0BC1 = 1;                      % 左边界值 (Pa)
u0BC2 = -1;                     % 右边界值 (Pa)

2. 形状函数构造

在 modelLoss 函数中：

phi1 = 1-X;
phi2 = X;

phi_eqv = phi1.*phi2;

3. 试探解构造

% 试探神经网络
G = phi1*U0(1)+phi2*U0(2)+phi_eqv.*U;

边界条件满足机制

这种构造方式的巧妙之处在于：

在左边界 (x=0) 时：

phi1 = 1-0 = 1
phi2 = 0
phi_eqv = 1×0 = 0
因此：G = 1×U0(1) + 0×U0(2) + 0×U = U0(1) = 1

在右边界 (x=1) 时：

phi1 = 1-1 = 0
phi2 = 1
phi_eqv = 0×1 = 0
因此：G = 0×U0(1) + 1×U0(2) + 0×U = U0(2) = -1

方法优势

自动满足边界条件：无论神经网络输出什么值，试探解在边界处都严格等于指定的边界值
无需惩罚项：不需要在损失函数中添加边界条件的惩罚项，避免了权重平衡问题
数学保证：边界条件的满足是通过数学构造保证的，而不是通过优化近似获得的
计算效率：减少了损失函数的复杂性，提高了训练效率

这种试探解方法是现代PINNs中处理边界条件的一种优雅且有效的方式，特别适用于狄利克雷边界条件（给定边界值）的问题。

3. 代码详细解读

3.1 参数设置和数据生成

% 整合的1维亥姆霍兹方程求解程序 - 使用物理信息神经网络(PINNs)
% 包含所有必要的函数

clc;
clear all;

%% 生成训练数据
freq = 2000;                     % 频率 (Hz)
c0 = 340;                       % 空气中声速 (m/s)
k = 2*pi*freq/c0;               % 波数 (1/m)

解读：

freq = 2000：设置声波频率为2000Hz
c0 = 340：空气中标准声速
k = 2*pi*freq/c0：根据色散关系计算波数

x0BC1 = 0;                      % 左边界 (m)
x0BC2 = 1;                      % 右边界 (m)

u0BC1 = 1;                      % 左边界值 (Pa)
u0BC2 = -1;                     % 右边界值 (Pa)

X0 = [x0BC1 x0BC2];     
U0 = [u0BC1 u0BC2];

解读：

定义求解区域：x ∈ [0,1]
设置边界条件：u(0)=1, u(1)=-1
将边界位置和值存储为数组

numInternalCollocationPoints = 14000;

pointSet = sobolset(1);                   % 高度均匀填充空间的二进制数字序列
points = net(pointSet,numInternalCollocationPoints);    % 生成准随机点集

dataX = points;         % 在0和1之间创建随机x数据点

解读：

numInternalCollocationPoints = 14000：设置配点法的内部点数量
sobolset(1)：创建一维Sobol序列生成器，产生高度均匀分布的准随机数
net(pointSet,14000)：生成14000个在[0,1]区间均匀分布的配点
这些配点用于在物理方程约束下训练神经网络

3.2 神经网络架构定义

%% 定义深度学习模型
numLayers = 5;
numNeurons = 90;

parameters = buildNet(numLayers,numNeurons);

解读：

创建5层神经网络，每层90个神经元
调用buildNet函数构建网络参数

buildNet函数详解：

function parameters = buildNet(numLayers,numNeurons)
% 构建神经网络参数

parameters = struct;

% 输入层
sz = [numNeurons 1];
parameters.fc1_Weights = initializeHe(sz,1,"double");
parameters.fc1_Bias = initializeZeros([numNeurons 1],"double");

解读：

输入层：1个输入(x坐标) → 90个神经元
使用He初始化方法初始化权重
偏置初始化为零

% 隐藏层
for layerNumber=2:numLayers-1
    name = "fc"+layerNumber;

    sz = [numNeurons numNeurons];
    numIn = numNeurons;
    parameters.(name + "_Weights") = initializeHe(sz,numIn,"double");
    parameters.(name + "_Bias") = initializeZeros([numNeurons 1],"double");
end

解读：

隐藏层：90个神经元 → 90个神经元
循环创建第2到第4层(共3个隐藏层)
每层都使用He初始化

% 输出层
sz = [1 numNeurons];
numIn = numNeurons;
parameters.("fc" + numLayers + "_Weights") = initializeHe(sz,numIn,"double");
parameters.("fc" + numLayers + "_Bias") = initializeZeros([1 1],"double");

解读：

输出层：90个神经元 → 1个输出
网络最终输出一个标量值N(x)

3.3 前向传播模型

function U = model(parameters,X)
% 神经网络模型前向传播

numLayers = numel(fieldnames(parameters))/2;

% 第一个全连接操作
weights = parameters.fc1_Weights;
bias = parameters.fc1_Bias;
U = fullyconnect(X,weights,bias);

% 对剩余层进行sin激活和全连接操作
fori=2:numLayers
    name = "fc" + i;

    U = sin(U);

    weights = parameters.(name + "_Weights");
    bias = parameters.(name + "_Bias");
    U = fullyconnect(U, weights, bias);
end

解读：

numLayers = numel(fieldnames(parameters))/2：计算网络层数(每层有权重和偏置两个参数)
第一层：线性变换但不激活
后续层：先sin激活，再线性变换
重要：使用sin激活函数有利于学习周期性解

3.4 损失函数和梯度计算

function [loss,gradients] = modelLoss(parameters,X,U0,freq,c0)
% 计算模型损失和梯度
% 使用初始条件进行预测
k = 2*pi*freq/c0;
U = model(parameters,X);
phi1 = 1-X;
phi2 = X;

phi_eqv = phi1.*phi2;

% 试探神经网络
G = phi1*U0(1)+phi2*U0(2)+phi_eqv.*U;

解读：

U = model(parameters,X)：神经网络输出N(x)
phi1 = 1-X：左边界形状函数
phi2 = X：右边界形状函数
phi_eqv = phi1.*phi2：等效形状函数 = x(1-x)
G = ...：构造试探函数，自动满足边界条件

% 计算对X的一阶导数
Gx = dlgradient(sum(G,'all'),X,'EnableHigherDerivatives',true);

% 计算对X的二阶导数
Gxx = dlgradient(sum(Gx,'all'),X,'EnableHigherDerivatives',true);

% 计算损失
f = Gxx+k^2*G;
zeroTarget = zeros(size(f),"like",f);
loss = l2loss(f, zeroTarget);

% 计算相对于可学习参数的梯度
gradients = dlgradient(loss,parameters);

解读：

dlgradient：自动微分计算梯度
EnableHigherDerivatives=true：允许计算高阶导数
Gx：dG/dx (一阶导数)
Gxx：d²G/dx² (二阶导数)
f = Gxx+k^2*G：亥姆霍兹方程左边的表达式
loss = l2loss(f, zeroTarget)：L2损失，要求f≈0
gradients：损失函数对网络参数的梯度

3.5 优化设置和训练

%% 指定优化选项
options = optimoptions("fmincon", ...
    HessianApproximation="lbfgs", ...
    MaxIterations=14000, ...
    MaxFunctionEvaluations=14000, ...
    OptimalityTolerance=1e-3, ...
    SpecifyObjectiveGradient=true, ...
    Display='iter');

解读：

使用fmincon进行约束优化
lbfgs：有限内存BFGS算法，适合大规模优化
MaxIterations=14000：最大迭代次数
SpecifyObjectiveGradient=true：使用解析梯度加速收敛

%% 训练网络
start = tic;

[parametersV,parameterNames,parameterSizes] = parameterStructToVector(parameters);
parametersV = extractdata(parametersV);

%% 将变量转换为深度学习变量
X = dlarray(dataX,"BC");
X0 = dlarray(X0,"CB");
U0 = dlarray(U0,"CB");

objFun = @(parameters) objectiveFunction(parameters,X,U0,parameterNames,parameterSizes,freq,c0);

parametersV = fmincon(objFun,parametersV,[],[],[],[],[],[],[],options);

parameters = parameterVectorToStruct(parametersV,parameterNames,parameterSizes);

toc(start)

解读：

parameterStructToVector：将结构体参数转换为向量形式
dlarray：创建深度学习数组，支持自动微分
"BC"和"CB"：指定数组维度标签(B=batch, C=channel)
objFun：创建目标函数句柄
fmincon：执行优化求解
训练完成后将参数转换回结构体形式

3.6 模型评估和可视化

%% 评估模型精度
numPredictions = 500;
XTest = linspace(x0BC1,x0BC2,numPredictions);

dlXTest = dlarray(XTest,'CB');
U = model(parameters,dlXTest);

% 空间函数
phi1_test = 1-dlXTest;
phi2_test = dlXTest;
phi_eqv_test = phi1_test.*phi2_test;

% 预测解
dlUPred = phi1_test*U0(1)+phi2_test*U0(2)+phi_eqv_test.*U;

% 真实解
UTest = solve1DWaveEqn(XTest,k);

% 计算误差
err = norm(extractdata(dlUPred) - UTest) / norm(UTest);

解读：

在500个测试点上评估模型
使用训练好的参数计算预测解
调用解析解函数计算真实解
计算相对L2误差：||u_pred - u_true||/||u_true||

f1 = figure;
% 绘制预测解
plot(XTest,extractdata(dlUPred),'-','LineWidth',2);

% 绘制真实解
hold on
plot(XTest, UTest, '--','LineWidth',2)
hold off

xlabel('x (m)','FontSize',14,'FontWeight','bold')
ylabel('声压 (Pa)','FontSize',14,'FontWeight','bold')
title("频率 = " + freq + " Hz;" + " 误差 = " + gather(err));

legend('预测解','真实解')

解读：

可视化预测解和真实解的对比
显示频率和相对误差信息

3.7 解析解函数

function U = solve1DWaveEqn(X,k)
% 1维波动方程的解析解

% 初始化解
U = zeros(size(X));

% 循环遍历x值
for i = 1:numel(X)
    x = X(i);
    U(i) = cos(k*x)-(csc(k)+cot(k))*sin(k*x);
end