基于PINN求解薄膜四阶抛物方程自相似解的完整流程
本文结合前序PINN框架核心逻辑与薄膜破裂四阶抛物方程特性,设计一套从坐标变换到稳定解验证的全流程方案,精准求解方程的自相似解,具体步骤如下:
一、问题前置:明确薄膜四阶抛物方程及自相似性核心
1. 目标方程(无量纲化)
薄膜在范德华力与表面张力作用下的演化方程为四阶非线性抛物方程(轴对称/二维通用形式): - 轴对称: $$ \frac{\partial h}{\partial t}+\frac{1}{r} \frac{\partial}{\partial r}\left[r h^{3} \frac{\partial}{\partial r}\left(\frac{1}{r} \frac{\partial}{\partial r}\left(r \frac{\partial h}{\partial r}\right)\right)+\frac{r}{h} \frac{\partial h}{\partial r}\right]=0 $$
-
二维: $$ \frac{\partial h}{\partial t}+\frac{\partial}{\partial x}\left[\left(h^{3} \frac{\partial^{3} h}{\partial x^{3}}+\frac{1}{h} \frac{\partial h}{\partial x}\right)\right]=0 $$
-
核心目标:找到满足\(h(r,t)=(t_R-t)^\lambda H(\eta)\)(\(\eta=\frac{r}{(t_R-t)^\beta}\))的自相似解,确定标度参数 \(\lambda\)、\(\beta\) 及自相似剖面 \(H(\eta)\)。
2. 自相似性先验认知
根据文献已知薄膜破裂的自相似标度 \(\lambda\)、\(\beta\),可以用于后期验证实验结果,本流程将通过PINN自主识别求解。
根据物理约束条件,H需要额外满足两个条件 1 原点正则性(\(H'(0)=H'''(0)=0\)) 2 远场渐近性(\(H(\eta)\sim A\eta^{1/2}\),\(A\) 为常数,也是需要求解的)。
接着,就是根据ode来确定约束。
3.所要求解的重要参数
\(\eta\),\(\beta\)作为尺度参数,于后文的识别部分来求解
A作为远场渐进性的参数,将被整合进函数的要求作为PINN求解的参数
二、步骤1:引入自相似坐标,转化稳态常微分方程
1. 构造自相似拟设(Ansatz)
基于四阶抛物方程的量纲分析与自相似特性,定义: - 自相似坐标:\(\eta=\frac{r}{(t_R-t)^\beta}\)(\(r\) 为径向/水平坐标,\(t_R\) 为破裂时刻,\(\beta\) 为坐标标度参数) - 薄膜厚度拟设:\(h(r,t)=(t_R-t)^\lambda H(\eta)\)(\(\lambda\) 为厚度标度参数,\(H(\eta)\) 为无量纲自相似剖面)
2. 代入方程消去时间项
利用链式法则计算时间导数 \(\partial_t h\) 和空间导数 \(\partial_r h,\partial_r^4 h\),代入原始四阶抛物方程: -时间导数:\(\partial_t h=-\lambda(t_R-t)^{\lambda-1}H(\eta)+\beta(t_R-t)^{\lambda-1}\eta H'(\eta)\) - 空间导数:\(\partial_r h=(t_R-t)^{\lambda-\beta}H'(\eta)\).高阶导数同理按链式法则展开 - 化简结果:消去 \((t_R-t)^{\lambda-1-\beta^4}\) 等标度因子后,得到仅关于 \(\eta\) 的**稳态常微分方程(ODE)**: - 轴对称场景:\(\frac{1}{5}(H-2\eta H')=\frac{1}{\eta}\left[\frac{\eta H'}{H}+\eta H^3\left(\frac{(\eta H')'}{\eta}\right)'\right]'\)(已代入文献标度 \(\lambda=1/5,\beta=2/5\),PINN可自主优化验证) - 二维对称场景:\(H\)
三、步骤2:嵌入数学结构,构造约束化解拟设
1. 坐标紧致化(处理无限域)
自相似坐标 \(\eta\in[0,+\infty)\) 为无限域,构造非线性变换映射至有限域 \(q\in(0,1]\): 由正则性约束,还需要满足关于\(\eta\)为偶函数的条件 $$ q=\frac{\eta^2}{\sqrt{1+\eta^4}} \quad \text{或} \quad q=1-e^{-\eta^2/5} $$ (适配远场 \(H(\eta)\sim\eta^{1/2}\) 的渐近衰减,避免边界数值失真)
2. 结构化解拟设(硬编码约束)
将原点正则性、符号约束等数学结构嵌入拟设,避免无效搜索: $$ H(\eta)=\exp\left(NN(q)\right) \cdot \sqrt{1+\eta^2} \cdot \frac{1}{(1+\eta^2)^{1/4}} $$
- 解析项作用: -- \(H_0\):原点厚度初始猜测(文献值 \(H(0)\approx0.768\),轴对称场景,在此处需要作为哦参数输入); -- \(\sqrt{1+\eta^2}\):适配远场 \(H(\eta)\sim\eta^{1/2}\) 的渐近行为; -- \(\frac{1}{(1+\eta^2)^{1/4}}\):压制高频振荡,保证解的光滑性; 这两项实际上就是合起来,在无限远处凑成\(\eta^{1/2}\),而且最好还要满足不会受到导数条件的限制
-- \(\exp(NN(q))\):神经网络输出项(隐含层用ReLU+Sigmoid组合,保证梯度平滑,而且最终有界,而且q这一项在\(\eta->\infty\)时是1,\(e^{NN(1)}\)也是一个有界量,用于拟合A这一项,而且不会影响无穷远处的阶数渐进
四、步骤3:构造二级神经网络(MSNN)
分阶段逼近解 采用两阶段神经网络架构,先粗拟合再精细化修正,适配四阶方程的高梯度特性:
1. 第一阶段网络(粗拟合核心剖面)
- 网络结构:输入层(\(q\in(0,1]\))→ 隐含层(3层,每层128神经元,激活函数ReLU)→ 输出层(1个神经元,输出 \(NN_1(q)\)),再带入原始式子,得到\(H_1(\eta)=\exp\left(NN(q)\right) \cdot \sqrt{1+\eta^2} \cdot \frac{1}{(1+\eta^2)^{1/4}}\)
- 训练目标:最小化自相似ODE的残差 \(L_1=\frac{1}{N}\sum_{i=1}^N \frac{|R(\eta_i; NN_1)|^2}{exp(NN_1)}\)(\(R\) 为ODE残差,\(N\) 为配点数);
- 各参数解释:分子部分残差,将\(\eta_i,NN_1\)代入\(H\),
- 输出结果:初步自相似剖面 \(H_1(\eta)\)、标度参数初始值 \(\lambda_1,\beta_1\)。
2. 第二阶段网络(修正高梯度区域)
- 误差拟设:定义误差项 \(\Delta H(\eta)=H(\eta)-H_1(\eta)\),代入自相似ODE线性化得到误差方程;
- 网络结构:输入层(\(q\in(0,1]\))→ 隐含层(4层,每层256神经元,激活函数GELU)→ 输出层(1个神经元,输出 \(\Delta H(NN_2(q))\));
- 训练目标:最小化误差方程残差 \(L_2=\frac{1}{N}\sum_{i=1}^N \left|\frac{R(\eta_i; H_1+\Delta H)}{1+|\nabla H_1(\eta_i)|}\right|^2\)(引入梯度归一化,压制高梯度区域残差主导);
- 最终输出:\(H(\eta)=H_1(\eta)+\Delta H(\eta)\),提升高梯度区域(如 \(\eta\approx3\sim5\))的拟合精度。
五、步骤4:配置二阶优化器,加速收敛
摒弃一阶优化器(如Adam),采用**全矩阵高斯-牛顿优化器**(基于kfac-jax实现),适配四阶方程的强非线性:
1. 优化器参数设置
- 学习率:自适应调整(初始值 \(1e-3\),残差下降停滞时衰减至 \(1e-5\));
- 动量系数:0.95(平衡全局搜索与局部收敛);
- 正则化:L2正则化(系数 \(1e-8\)),避免过拟合。
2. 训练策略
- 第一阶段:训练5000迭代,残差收敛至 \(O(10^{-6})\) 后停止;
- 第二阶段:加载第一阶段参数,继续训练8000迭代,目标残差 \(O(10^{-10})\)(接近双精度舍入误差); -- 优势:二阶优化器直接利用残差的Hessian信息,解决四阶方程残差曲面非凸导致的收敛停滞问题。
六、步骤5:混合配点操作,聚焦关键区域
采用“位置采样+自适应采样”混合策略,确保配点覆盖性与关键区域密度:
1. 位置采样(全域覆盖)
-采样区域:\(\eta\in[0,20]\)(覆盖远场渐近区域,\(\eta>20\) 时 \(H(\eta)\sim\eta^{1/2}\),误差可忽略); - 采样分布:\(\eta\in[0,5]\) 均匀采样(原点及高梯度区域),\(\eta\in[5,20]\) 对数采样(远场区域,适配渐近衰减); - 采样数量:基础配点500个。
2. 自适应采样(聚焦高残差区)
- 采样规则:每100迭代计算当前配点的残差 \(R(\eta_i)\),按残差平方加权采样,补充200个高残差配点(如 \(\eta\approx2\sim4\) 的梯度突变区域);
- 动态更新:每次补充配点后,删除残差小于 \(1e-11\) 的低价值配点,保持总配点数700个,提升训练效率。
七、步骤6:识别标度参数(\(\lambda,\beta\)),双策略验证
1. 解析推断(初步定位)
利用原点正则性条件(\(H'(0)=H'''(0)=0\)),将 \(\eta\to0\) 代入自相似ODE: - 泰勒展开 \(H(\eta)=H(0)+\frac{1}{2}H''(0)\eta^2+o(\eta^2)\),代入ODE后得到关于 \(\lambda,\beta\) 的代数方程: $$ \lambda H(0) = \beta H(0) + \text{低阶项} $$ - 结合量纲分析,初步推断 \(\lambda\approx1/5\)、\(\beta\approx2/5\),作为数值搜索初始值。
2. 迭代漏斗搜索(精准优化)
以残差最小化为目标,采用割线搜索优化 \(\lambda,\beta\): - 搜索范围:\(\lambda\in[0.15,0.25]\)、\(\beta\in[0.35,0.45]\)(围绕文献值小幅拓展); - 目标函数:\(L(\lambda,\beta)=\max\left(\frac{1}{N}\sum|R(\eta_i)|^2, \left|\frac{H(\eta\to\infty)}{\eta^{1/2}}-A\right|\right)\)(\(A\) 为远场系数,文献值 \(A\approx0.676\)); - 收敛判据:\(\Delta L<1e-12\),最终输出 \(\lambda^*,\beta^*\)(理论上应逼近 \(\lambda=1/5\)、\(\beta=2/5\))。
八、步骤7:求解稳定点,验证解的有效性
1. 稳定点判定标准 自相似解的稳定点需满足:
- 残差条件:全域ODE残差 \(\max|R(\eta)|<1e-10\);
- 远场匹配:\(|H(\eta=20)/20^{1/2}-A|<1e-4\);
- 正则性条件:\(|H'(0)|<1e-8\)、\(|H'''(0)|<1e-8\)。
2. 稳定性验证(适配薄膜破裂特性)
- 对比文献数据:轴对称场景下,\(H(0)\approx0.768\)、\(H''(0)\approx0.169\),与表1第一组解误差小于 \(1e-3\);
- 振荡特性:稳定解的 \(H'''(\eta)\) 应最少振荡(文献指出“最小振荡曲率剖面为物理可观测解”),需验证 \(H'''(\eta)\) 在 \(\eta\in[0,10]\) 内无额外零点。
3. 结果输出
- 最终自相似解:\(h(r,t)=(t_R-t)^{\lambda^*} H\left(\frac{r}{(t_R-t)^{\beta^*}}\right)\);
- 关键参数:\(\lambda^*,\beta^*,H(0),H''(0),A\)(远场系数);
- 验证指标:残差分布图、\(H(\eta)\) 与文献解的对比图、远场渐近拟合度。
核心流程总结表
| 步骤 | 核心操作 | 关键目标 |
|---|---|---|
| 自相似坐标引入 | 构造 \(\eta=\frac{r}{(t_R-t)^\beta}\)、\(h=(t_R-t)^\lambda H(\eta)\) | 转化为稳态ODE |
| 数学结构嵌入 | 坐标紧致化+结构化拟设(硬编码正则性、渐近性) | 缩小解空间,保证物理合理性 |
| 二级神经网络 | 粗拟合(ReLU)+ 精细化修正(GELU+梯度归一化) | 适配四阶方程高梯度特性 |
| 二阶优化器 | 高斯-牛顿优化(kfac-jax)+ 自适应学习率 | 突破非凸残差曲面的收敛停滞 |
| 混合配点 | 位置采样(全域覆盖)+ 自适应采样(高残差区聚焦) | 平衡覆盖性与计算效率 |
| \(\lambda/\beta\) 识别 | 解析推断(初步定位)+ 迭代漏斗搜索(精准优化) | 确定自相似标度参数 |
| 稳定点求解 | 残差+正则性+远场匹配验证,对比文献数据 | 输出物理可观测的稳定自相似解 |
| 该流程完全复用PINN求解自相似解的核心框架,同时针对薄膜四阶抛物方程的高非线性、无限域、正则性约束等特性定制化调整,可高效得到与文献一致的稳定自相似解。 |