偏微分方程(PDE)是研究各种自然现象的重要工具，被广泛用于解释各种物理规律。此外，许多工程和技术问题可以用偏微分方程进行建模和分析，如尾流湍流、光纤通信、大气污染物扩散等。因此，偏微分方程的研究对于航空航天、数值天气预报等许多领域都具有重要意义。目前，偏微分方程与其他许多学科的联系日益紧密，相互促进。因此，对偏微分方程的研究具有重要意义。然而，偏微分方程研究中的一个主要困难是往往不可能获得解析解。为此，相关研究人员提出了各种求解偏微分方程的数值方法，如有限差分法、有限元法、有限体积法等。数值方法极大地促进了偏微分方程的研究。如今，这些方法已经得到了广泛的应用，并且还在不断完善。与此同时，研究人员也在努力开发解决偏微分方程的新方法和新工具。
如今，越来越多的研究人员使用深度学习方法来研究偏微分方程。在众多的研究中，Raissi1等人提出的物理信息神经网络(PINN) 是一项不可忽视的重要工作。该神经网络模型考虑了偏微分方程所包含的物理规律，并将其作为正则项编码到神经网络中，提高了神经网络模型的性能。
微分方程求解介绍
传统的前向神经网络完全基于数据驱动的方法，没有考虑数据中包含的物理规律。因此，为了得到合理的模型，往往需要大量的数据对神经网络进行训练。相反，物理信息神经网络通过强迫网络输出满足相应的偏微分方程，将物理信息引入网络。具体地说，通过在损失函数中加入偏微分方程的正则化，使得模型在训练过程中考虑了物理规律。 这种处理使得训练过程需要较少的数据，并且加快了训练过程。物理信息神经网络不仅可以用来解决正向问题，即求偏微分方程的近似解，还可以用来求解反问题，即从训练数据中获得偏微分方程的参数。.
它由输入层、一个或多个隐藏层和一个输出层组成，每个层包含一个或多个人工神经元
迭代法求解微分方程
对于任意一个二阶微分方程，我们都可以用这个方程表示出
G
(
x
?
,
ψ
(
x
?
)
,
?
ψ
(
x
?
)
,
?
2
ψ
(
x
?
)
)
=
0
,
x
?
∈
D
G(\vec{x},\psi (\vec{x}), \nabla\psi (\vec{x}), \nabla^2\psi (\vec{x}))=0,\vec{x}\in D
G(x,ψ(x),?ψ(x),?2ψ(x))=0,x∈D
求解目的就是找出这样的一个方程： ψ ( x ? ) \psi (\vec{x}) ψ(x)，能够满足以上的 G ( x ? ) G(\vec{x}) G(x)函数。
求解微分方程的数值解，常常使用计算机辅助求解。
对于计算机求解，第一步要将其离散化处理：
G
(
x
i
?
,
ψ
(
x
i
?
)
,
?
ψ
(
x
i
?
)
,
?
2
ψ
(
x
i
?
)
)
=
0
,
?
x
i
?
∈
D
^
G(\vec{x_i},\psi (\vec{x_i}), \nabla\psi (\vec{x_i}), \nabla^2\psi (\vec{x_i}))=0, \forall\vec{x_i}\in \hat{D}
G(xi??,ψ(xi??),?ψ(xi??),?2ψ(xi??))=0,?xi??∈D^
离散化后，微分项就变成差分项：
?
ψ
?
x
≈
ψ
k
?
ψ
k
?
1
Δ
x
,
?
2
ψ
?
x
2
≈
ψ
k
+
1
?
2
ψ
k
+
ψ
k
?
1
Δ
x
2
\frac{\partial\psi}{\partial x} \approx\frac{\psi^k-\psi^{k-1}}{\Delta x}, \frac{\partial^2\psi}{\partial x^2} \approx\frac{\psi^{k+1}-2\psi^{k}+\psi^{k-1}}{\Delta x^2}
?x?ψ?≈Δxψk?ψk?1?,?x2?2ψ?≈Δx2ψk+1?2ψk+ψk?1?
然后使用计算机按照步长
Δ
x
{\Delta x}
Δx迭代求解即可。
PINN法求解微分方程
人工神经网络若要求解该方程，那就设
ψ
(
x
?
)
\psi (\vec{x})
ψ(x)为如下形式：
ψ
(
x
?
)
≡
N
m
(
x
?
,
{
w
,
b
?
}
)
\psi (\vec{x})\equiv N_m(\vec{x},\begin{Bmatrix} w ,\vec{b}\end{Bmatrix})
ψ(x)≡Nm?(x,{w,b?})
上式中，
N
m
N_m
Nm?表示神经网络的第
m
m
m层的输出，在这里也就表示整个前向神经网络的输出。
{
w
,
b
?
}
\begin{Bmatrix} w ,\vec{b}\end{Bmatrix}
{w,b?}表示全连接神经网络的权重与偏置。
微分项借助神经网络的连接规则，采用 反向传播(Back Propagation, BP)算法 进行计算。
传播过程中满足如下公式：
?
(
n
)
ψ
?
ψ
(
j
)
=
∑
i
,
j
∈
P
a
(
ψ
(
i
)
)
?
ψ
(
n
)
?
ψ
(
i
)
?
ψ
(
i
)
?
ψ
(
j
)
\frac{\partial^{(n)}\psi}{\partial \psi^{(j)}}= \sum_{i,j\in Pa(\psi^{(i)})}\frac{\partial\psi^{(n)}}{\partial \psi^{(i)}}\frac{\partial\psi^{(i)}}{\partial \psi^{(j)}}
?ψ(j)?(n)ψ?=i,j∈Pa(ψ(i))∑??ψ(i)?ψ(n)??ψ(j)?ψ(i)?
具体的，对于整个网络来说，输出对输入参数的微分即为：
?
k
ψ
?
x
k
=
∑
i
=
1
m
v
i
w
i
,
j
k
σ
i
(
k
)
\frac{\partial^{k}\psi}{\partial x^{k}}= \sum_{i=1}^mv_iw_{i,j}^k\sigma_i^{(k)}
?xk?kψ?=i=1∑m?vi?wi,jk?σi(k)?
其中
σ
i
\sigma_i
σi?表示第
i
i
i层网络的激活函数。
PINN方法是将物理公式编码到神经网络的正则化项中，其损失函数包含物理背景。
L
(
{
w
,
b
?
}
)
=
1
i
m
a
x
∑
i
,
m
∣
∣
F
^
(
x
i
?
,
ψ
^
m
(
x
?
i
)
,
…
?
2
ψ
^
(
x
i
?
)
)
∣
∣
+
∑
B
.
C
.
∣
∣
?
p
ψ
^
m
(
x
?
b
)
?
K
(
x
b
?
)
∣
∣
L(\begin{Bmatrix} w ,\vec{b}\end{Bmatrix})=\frac{1}{i_{max}}\sum_{i,m}\left|\right|\hat F(\vec{x^i},\hat \psi_m(\vec x^i), \dots \nabla^2\hat\psi(\vec{x^i}))\left|\right|+\sum_{B.C.}\left|\right|\nabla^p\hat\psi_m(\vec x^b)-K(\vec{x_b})\left|\right|
L({w,b?})=imax?1?i,m∑?∣∣F^(xi,ψ^?m?(xi),…?2ψ^?(xi))∣∣+B.C.∑?∣∣?pψ^?m?(xb)?K(xb??)∣∣
在这个损失函数中，第一项表示由物理公式产生的损失，第二项表示由初始值(Initial conditions)\边界值(Boundary conditions)带来的误差，将损失函数
L
(
{
w
,
b
?
}
)
L(\begin{Bmatrix} w ,\vec{b}\end{Bmatrix})
L({w,b?})的损失值降到最小，即表示神经网络的输出值能满足边界条件的情况下，也能满足物理公式！ 也就是建立起输出到输入的前向神经网络，该网络能够满足求解需要。
方法验证
伯格斯方程验证
伯格斯方程(Burgers equation) 是一个模拟冲击波的传播和反射的非线性偏微分方程。伯格斯方程也称为粘性伯格斯方程，只适用于一个空间维度。通过求解该方程验证PINN方法在一阶偏微分下的有效性。
伯格斯方程的表达形式：
ψ
t
+
ψ
ψ
x
?
(
0.01
/
π
)
ψ
x
x
=
0
,
x
∈
[
?
1
,
1
]
,
t
∈
[
0
,
1
]
\psi_t+\psi\psi_x-(0.01/\pi)\psi_{xx}=0,x\in [-1,1],t\in[0,1]
ψt?+ψψx??(0.01/π)ψxx?=0,x∈[?1,1],t∈[0,1]
边界条件如下：
B
.
C
.
{
ψ
(
0
,
x
)
=
?
s
i
n
(
π
x
)
,
ψ
(
t
,
?
1
)
=
u
(
t
,
1
)
=
0
B.C. \begin{cases} \psi(0,x)=-sin(\pi x), & \\ \psi(t,-1)=u(t,1)=0 \end{cases}
B.C.{ψ(0,x)=?sin(πx),ψ(t,?1)=u(t,1)=0?
下图为网络结构图，绿色部分是前向全连接神经网络，红色部分是偏微分方程构成的内嵌网络部分。全连接神经网络的输入即为伯格斯方程的两个自变量，输出为方程的解。灰色部分为内嵌物理公式与边界条件构成的损失值部分。
二阶偏微分方程验证
使用PINN网络求解二阶偏微分方程。偏微分方程的方程式由下面的方程给出。
?
2
ψ
(
x
,
y
)
=
e
?
x
(
x
?
2
+
y
3
+
6
y
)
,
x
,
y
∈
[
0
,
1
]
\nabla^2\psi(x,y)=e^{-x}(x-2+y^3+6y),x,y\in [0,1]
?2ψ(x,y)=e?x(x?2+y3+6y),x,y∈[0,1]
边界条件如下：
B
.
C
.
{
ψ
(
0
,
y
)
=
y
3
,
ψ
(
i
,
y
)
=
(
1
+
y
3
)
e
?
1
,
ψ
(
x
,
0
)
=
x
e
?
x
,
ψ
(
x
,
1
)
=
e
?
x
(
x
+
1
)
.
B.C. \begin{cases} \psi(0,y)=y^3, & \\ \psi(i,y)=(1+y^3)e^{-1},\\ \psi(x,0)=xe^{-x},\\ \psi(x,1)=e^{-x}(x+1). & \end{cases}
B.C.??????????ψ(0,y)=y3,ψ(i,y)=(1+y3)e?1,ψ(x,0)=xe?x,ψ(x,1)=e?x(x+1).??
这个方程有明确的解析解，解析解的表达式为：
ψ
a
(
x
,
y
)
=
e
?
x
(
x
+
y
3
)
\psi_a(x,y)=e^{-x}(x+y^3)
ψa?(x,y)=e?x(x+y3)
根据表达式绘制出解析解的图像：
二阶偏微分方程解析解图像。x1和x2分别表示两个输入。
使用神经网络求解的偏微分方程图像如图：
解析解与神经网络求解之间各个点的绝对误差使用下图展示。图中可以分析出PINN网络对求解二阶偏微分方程的准确度较高，可用于二阶偏微分方程的求解。
代码展示
下面展示求解PINN偏微分方程使用的代码，该代码求解的方程为一维形式的Allen-Cahn方程。更改 loss_ped() 方法代码可方便求解其他偏微分方程。
"""
Author: Minglang Yin
PINNs (physical-informed neural network) for solving time-dependent Allen-Cahn equation (1D).
u_t - 0.0001*u_xx + 5*u^3 - 5*u = 0, x in [-1, 1], t in [0, 1]
with
u(0, x) = x^2*cos(pi*x)
u(t, -1) = u(t, 1)
u_x(t, -1) = u_x(t, 1)
Input: [t, x]
Output: [u]
"""
import sys
# sys.path.insert(0,'../../Utils')
import torch
import torch.nn as nn
import numpy as np
import time
import random
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
import scipy.io
# from plotting import newfig, savefig
torch.manual_seed(1234)
np.random.seed(1234)
class Model(nn.Module):
def __init__(self):
super(Model, self).__init__()
self.net = nn.Sequential()#使用这个函数快速搭建神经网络
self.net.add_module('linear_layer_1', nn.Linear(2, 20))
self.net.add_module('tanh_layer_1', nn.Tanh())
for num in range(2, 5):
self.net.add_module('linear_layer_%d' % (num), nn.Linear(20, 20))
self.net.add_module('tanh_layer_%d' % (num), nn.Tanh())
self.net.add_module('linear_layer_50', nn.Linear(20, 1))
def forward(self, x):
return self.net(x)
def loss_pde(self, x):
u = self.net(x)
u_g = gradients(u, x)[0]
u_t, u_x = u_g[:, :1], u_g[:, 1:]
u_gg = gradients(u_g, x)[0]
u_xx = u_gg[:, 1:]
# loss = u_t - 0.0001*u_xx + 5.0*u**3 - 5.0*u
loss = u_t - 0.0001 * u_xx + u ** 3 - u
return (loss ** 2).mean()
def loss_f(self, x, u_f_train):
u_f_pred = self.net(x)
return ((u_f_pred - u_f_train) ** 2).mean()
def loss_bc(self, x_b_l_train, x_b_r_train):
u_b_l_pred = self.net(x_b_l_train)
u_b_r_pred = self.net(x_b_r_train)
u_b_l_pred_x = gradients(u_b_l_pred, x_b_l_train)[0][:, 1]
u_b_r_pred_x = gradients(u_b_r_pred, x_b_r_train)[0][:, 1]
return ((u_b_l_pred - u_b_r_pred) ** 2).mean() + ((u_b_l_pred_x - u_b_r_pred_x) ** 2).mean()
def loss_ic(self, x_i_train, u_i_train):
u_i_pred = self.net(x_i_train)
return ((u_i_pred - u_i_train) ** 2).mean()
def gradients(outputs, inputs):
return torch.autograd.grad(outputs, inputs, grad_outputs=torch.ones_like(outputs), create_graph=True)
def to_numpy(input):
if isinstance(input, torch.Tensor):
return input.detach().cpu().numpy()
elif isinstance(input, np.ndarray):
return input
else:
raise TypeError('Unknown type of input, expected torch.Tensor or ' \
'np.ndarray, but got {}'.format(type(input)))
def init_cond(x):
return np.sin(np.pi * x)
# return 1+np.cos(np.pi*x)
def main():
## parameters
cuda = 0
device = torch.device(f"cuda:{cuda}" if torch.cuda.is_available() else "cpu")
print(torch.cuda.is_available())
epochs = 100000
num_i_train = 200
num_b_train = 100
num_f_train = 10000
lr = 0.001
## pre-processing
data = scipy.io.loadmat('./AC.mat')
## x: array, x_grid: grid data, X: flatten data
t = data['tt'].flatten()[:, None]
x = data['x'].flatten()[:, None]
t_grid, x_grid = np.meshgrid(t, x)
exact_grid = np.real(data['uu'])
T = t_grid.flatten()[:, None]
X = x_grid.flatten()[:, None]
Exact = exact_grid.flatten()[:, None]
## Initial&Boundary data
id_i = np.random.choice(x.shape[0], num_i_train, replace=False)
id_b = np.random.choice(t.shape[0], num_b_train, replace=False)
id_f = np.random.choice(Exact.shape[0], num_f_train, replace=False)
x_i = x_grid[id_i, 0][:, None]
t_i = t_grid[id_i, 0][:, None]
x_i_train = np.hstack((t_i, x_i))
# u_i_train = init_cond(x_i)
u_i_train = exact_grid[id_i, 0][:, None]
x_b_l = x_grid[0, id_b][:, None]
x_b_r = x_grid[-1, id_b][:, None]
t_b_l = t_grid[0, id_b][:, None]
t_b_r = t_grid[-1, id_b][:, None]
x_b_l_train = np.hstack((t_b_l, x_b_l))
x_b_r_train = np.hstack((t_b_r, x_b_r))
x_f = X[id_f, 0][:, None]
t_f = T[id_f, 0][:, None]
x_f_train = np.hstack((t_f, x_f))
u_f_train = Exact[id_f, 0][:, None]
x_test = np.hstack((T, X))
## Form data tensor and send
x_i_train = torch.tensor(x_i_train, dtype=torch.float32).to(device)
x_b_l_train = torch.tensor(x_b_l_train, requires_grad=True, dtype=torch.float32).to(device)
x_b_r_train = torch.tensor(x_b_r_train, requires_grad=True, dtype=torch.float32).to(device)
x_f_train = torch.tensor(x_f_train, requires_grad=True, dtype=torch.float32).to(device)
x_test = torch.tensor(x_test, dtype=torch.float32).to(device)
u_i_train = torch.tensor(u_i_train, dtype=torch.float32).to(device)
u_f_train = torch.tensor(u_f_train, dtype=torch.float32).to(device)
## instantiate model
model = Model().to(device)
# Loss and optimizer
optimizer = torch.optim.Adam(model.parameters(), lr=lr)
# training
def train(epoch):
model.train()
def closure():
optimizer.zero_grad()
loss_pde = model.loss_pde(x_f_train)
loss_bc = model.loss_bc(x_b_l_train, x_b_r_train)
loss_ic = model.loss_ic(x_i_train, u_i_train)
loss = loss_pde + loss_bc + loss_ic
loss.backward()
return loss
loss = optimizer.step(closure)
loss_value = loss.item() if not isinstance(loss, float) else loss
print(f'epoch {epoch}: loss {loss_value:.6f}')
print('start training...')
tic = time.time()
for epoch in range(1, epochs + 1):
train(epoch)
toc = time.time()
print(f'total training time: {toc - tic}')
if __name__ == '__main__':
main()
结语
PINN方法是一种融合物理信息的神经网络求解偏微分方程的方法。该方法是对以往物理信息神经网络的改进。在该方法中，将偏微分方程中包含的物理定律作为正则化引入到神经网络中。这一改进激励神经网络更好地从有限数量的观测中学习偏微分方程的解。该方法对于物理模型明确但是数据难以测量的模型建模有重要意义。 物理信息神经网络最重要的任务是引入物理信息的合理规则化。物理信息的使用使神经网络能够更好地从较少的观测中学习偏微分方程的解。