neurodiffeq

神经差分

引文

@article{chen2020neurodiffeq,
  title={NeuroDiffEq: A Python package for solving differential equations with neural networks},
  author={Chen, Feiyu and Sondak, David and Protopapas, Pavlos and Mattheakis, Marios and Liu, Shuheng and Agarwal, Devansh and Di Giovanni, Marco},
  journal={Journal of Open Source Software},
  volume={5},
  number={46},
  pages={1931},
  year={2020}
}

---

您是否知道 Neurodiffeq 支持解决方案包并可用于解决逆向问题？看这里！

？已经熟悉 Neurodiffeq 了？ ？跳转至常见问题解答。

---

介绍

neurodiffeq是一个用神经网络求解微分方程的包。微分方程是将某些函数与其导数联系起来的方程。它们出现在各个科学和工程领域。传统上这些问题可以通过数值方法（例如有限差分、有限元）来解决。虽然这些方法是有效且充分的，但它们的可表达性受到其函数表示的限制。如果我们能够计算连续且可微的微分方程的解，那将会很有趣。

作为通用函数逼近器，人工神经网络已被证明具有求解具有某些初始/边界条件的常微分方程 (ODE) 和偏微分方程 (PDE) 的潜力。 neurodiffeq的目标是实现使用 ANN 求解微分方程的现有技术，使软件足够灵活，能够处理各种用户定义的问题。

安装

使用点

与大多数标准库一样， neurodiffeq托管在 PyPI 上。要安装最新的稳定版本，

 pip install -U Neurodiffeq # '-U' 表示更新到最新版本

手动

或者，您可以手动安装该库以尽早访问我们的新功能。对于想要为库做出贡献的开发人员，这是推荐的方式。

 git clone https://github.com/NeuroDiffGym/neurodiffeq.gitcd Neurodiffeq && pip install -r 要求
点安装。  # 要对库进行更改，请使用 `pip install -e .`pytesttests/ # 运行测试。选修的。

入门

我们很乐意帮助您解决任何问题。同时，您可以查看常见问题解答。

要查看neurodiffeq的完整教程和文档，请查看官方文档。

除了文档之外，我们最近还制作了带有幻灯片的快速演练演示视频。

用法示例

进口

从neurodiffeq导入diff从neurodiffeq.solvers导入Solver1D，Solver2D从neurodiffeq.conditions导入IVP，DirichletBVP2D从neurodiffeq.networks导入FCNN，SinActv

常微分方程系统示例

在这里，我们求解两个 ODE 的非线性系统，称为 Lotka-Volterra 方程。有两个未知函数（ u和v ）和一个自变量（ t ）。

 def ode_system(u, v, t): 返回 [diff(u,t)-(uu*v), diff(v,t)-(u*vv)]条件 = [IVP(t_0=0.0, u_0=1.5) ), IVP(t_0=0.0, u_0=1.0)]nets = [FCNN(actv=SinActv), FCNN(actv=SinActv)]solver = Solver1D(ode_system, 条件, t_min=0.1, t_max=12.0, nets=nets)solver.fit(max_epochs=3000)solution =solver.get_solution()

solution是一个可调用对象，您可以将 numpy 数组或火炬张量传递给它，例如

u, v = Solution(t, to_numpy=True) # t 可以是 np.ndarray 或 torch.Tensor

将u和v相对于它们的解析解绘制出来，结果如下：

偏微分方程系统示例

在这里，我们在矩形上求解具有狄利克雷边界条件的拉普拉斯方程。请注意，我们选择拉普拉斯方程是因为它计算解析解的简单性。在实践中，只要您足够好地调整求解器，您就可以尝试任何非线性、混沌偏微分方程。

求解 2-D PDE 系统与求解 ODE 非常相似，不同之处在于有两个变量x和y （用于边值问题）或x和t （用于初始边值问题），这两个变量均受支持。

 def pde_system(u, x, y):返回 [diff(u, x, order=2) + diff(u, y, order=2)]条件 = [DirichletBVP2D(x_min=0, x_min_val=lambda y: torch. sin(np.pi*y),x_max=1, x_max_val=lambda y: 0, y_min=0, y_min_val=lambda x: 0, y_max=1, y_max_val=lambda x: 0,
                       ）
]nets = [FCNN(n_input_units=2, n_output_units=1,hidden_​​units=(512,))]solver = Solver2D(pde_system, 条件, xy_min=(0, 0), xy_max=(1, 1), nets=nets) solver.fit(max_epochs=2000)solution =solver.get_solution()

二维偏微分方程的solution的签名略有不同。同样，它接受 numpy 数组或 torch 张量。

 u = 解(x, y, to_numpy=True)

在[0,1] × [0,1]上计算 u 会得到以下图

基于人工神经网络的解决方案	偏微分方程的残差

使用显示器

监视器是一种用于可视化 PDE/ODE 解决方案以及训练期间损失历史和自定义指标的工具。 Jupyter Notebooks 用户需要运行%matplotlib notebook魔法。对于 Jupyter Lab 用户，请尝试%matplotlib widget 。

从 Neurodiffeq.monitors 导入 Monitor1D...monitor = Monitor1D(t_min=0.0, t_max=12.0, check_every=100)solver.fit(..., 回调=[monitor.to_callback()])

您应该看到图每 100 个 epoch以及最后一个 epoch更新一次，显示两个图 - 一个用于区间[0,12]上的解决方案可视化，另一个用于损失历史记录（训练和验证）。

定制网络

为了方便起见，我们实现了FCNN全连接神经网络，其隐藏单元和激活函数可以定制。

 from Neurodiffeq.networks import FCNN# 默认：n_input_units=1, n_output_units=1,hidden_​​units=[32, 32],activation=torch.nn.Tanhnet1 = FCNN(n_input_units=..., n_output_units=...,hidden_​​units=[ ..., ..., ...]，激活=...)
 ...网络 = [net1, net2, ...]

FCNN通常是一个很好的起点。对于高级用户，求解器与任何自定义torch.nn.Module兼容。唯一的限制是：

1. 这些模块接受形状为(None, n_coords)的张量，输出形状为(None, 1)的张量。

2. nets中必须总共有n_funcs模块才能传递给solver = Solver(..., nets=nets) 。

实际上， neurodiffeq有一个single_net功能，它不遵守上述规则，这里不再介绍。

阅读有关构建您自己的网络（也称为模块）架构的 PyTorch 教程。

迁移学习

通过将old_solver.nets （火炬模块列表）序列化到磁盘，然后加载它们并传递给新的求解器，可以轻松完成迁移学习：

 old_solver.fit(max_epochs=...)# ... 将 `old_solver.nets` 转储到磁盘# ... 从磁盘加载网络，将它们存储在某个 `loaded_nets` 变量中 new_solver = Solver(..., nets=loaded_nets )new_solver.fit(max_epochs=...)

我们目前正在研究包装函数来保存/加载网络和求解器的其他内部变量。同时，您可以阅读有关保存和加载网络的 PyTorch 教程。

抽样策略

在 Neurodiffeq 中，通过最小化在域中一组点上评估的损失（ODE/PDE 残差）来训练网络。每次都会对点进行随机重新采样。要控制采样点的数量、分布和边界域，您可以指定自己的训练/验证generator 。

 from Neurodiffeq.generators import Generator1D# 默认 t_min=0.0, t_max=1.0, method='uniform', Noise_std=Noneg1 = Generator1D(size=..., t_min=..., t_max=..., method=.. .,noise_std=...)g2 = Generator1D(大小=...,t_min=...,t_max=...,方法=...,噪声_std=...)求解器 = Solver1D(..., train_generator=g1, valid_generator=g2)

以下是Generator1D的一些示例分布。

Generator1D(8192, 0.0, 1.0, method='uniform')	Generator1D(8192, -1.0, 0.0, method='log-spaced-noisy', noise_std=1e-3)

请注意，当同时指定train_generator和valid_generator时， Solver1D(...)中可以省略t_min和t_max 。事实上，即使你同时传递t_min 、 t_max 、 train_generator 、 valid_generator ， t_min和t_max仍然会被忽略。

组合发电机

生成器的另一个不错的功能是您可以连接它们，例如

g1 = Generator2D((16, 16), xy_min=(0, 0), xy_max=(1, 1))g2 = Generator2D((16, 16), xy_min=(1, 1), xy_max=(2, 2) ))g = g1 + g2

这里， g将是一个生成器，输出g1和g2的组合样本

g1	g2	g1 + g2

更高维度的采样

您可以使用Generator2D 、 Generator3D等对更高维度的采样点。但还有另一种方法

g1 = Generator1D(1024, 2.0, 3.0, method='uniform')g2 = Generator1D(1024, 0.1, 1.0, method='log-spaced-noisy',noise_std=0.001)g = g1 * g2

这里， g将是一个生成器，每次在 2-D 矩形(2,3) × (0.1,1)中生成 1024 个点。它们的x坐标是使用策略uniform从(2,3)绘制的，y坐标是使用策略log-spaced-noisy从(0.1,1)绘制的。

g1	g2	g1 * g2

解决捆绑和逆向问题

有时，一次求解一组方程是很有趣的。例如，您可能想要在初始条件u(0) = U0下求解du/dt + λu = 0形式的微分方程。您可能希望通过将所有λ和U0视为神经网络的输入来一次性解决此问题。

其中一种应用是化学反应，其中反应速率未知。不同的反应速率对应不同的解决方案，并且只有一种解决方案与观察到的数据点匹配。您可能有兴趣首先求解一组解，然后确定最佳反应速率（也称为方程参数）。第二步称为逆问题。

以下是如何使用neurodiffeq执行此操作的示例：

1. 假设我们有一个方程du/dt + λu = 0和初始条件u(0) = U0 ，其中λ和U0是未知常数。我们还有一组观察值t_obs和u_obs 。我们首先导入BundleSolver和BundleIVP这是获取解决方案包所必需的：

   从neurodiffeq.conditions导入BundleIVP从neurodiffeq.solvers导入BundleSolver1D导入matplotlib.pyplot作为plt导入numpy作为np导入torch从neurodiffeq导入diff

2. 我们确定输入t的域，以及参数λ和U0的域。我们还需要决定参数的顺序。即，哪个应该是第一个参数，哪个应该是第二个参数。出于本演示的目的，我们选择λ作为第一个参数（索引 0），选择U0作为第二个参数（索引 1）。跟踪参数的索引非常重要。

    T_MIN, T_MAX = 0, 1LAMBDA_MIN, LAMBDA_MAX = 3, 5 # 第一个参数，索引 = 0U0_MIN, U0_MAX = 0.2, 0.6 # 第二个参数，索引 = 1

3. 然后，我们像往常一样定义conditions和solver ，只是我们使用BundleIVP和BundleSolver1D而不是IVP和Solver1D 。这两个的接口与IVP和Solver1D非常相似。您可以在 API 参考中了解更多信息。

    # 方程参数位于输入（通常是时间和空间坐标）之后 diff_eq = lambda u, t, lmd: [diff(u, t) + lmd * u]# 关键字参数必须在 BundleIVP 中命名为“u_0”。如果你使用其他的东西，例如`y0`，`u0`等，它将不起作用。conditions = [BundleIVP(t_0=0, u_0=None, bundle_param_lookup={'u_0': 1}) # u_0 has索引 1]solver = BundleSolver1D(ode_system=diff_eq,条件=条件,t_min=T_MIN, t_max=T_MAX, theta_min=[LAMBDA_MIN, U0_MIN], # u_0 的索引为 1theta_max=[LAMBDA_MAX, U0_MAX], # u_0 的索引为 1eq_param_index=(0,), # λ 是唯一的方程参数，其索引为0n_batches_valid=1，
   ）

   由于λ是方程中的参数， U0是初始条件中的参数，因此我们必须在diff_eq中包含λ ，在条件中包含U0 。如果某个参数同时出现在方程和条件中，则必须将其包含在两个位置中。传递给BundleSovler1D的所有conditions元素都必须是Bundle*条件，即使它们没有参数。

4. 现在，我们可以像平常一样训练它并获得解决方案。

    solver.fit(max_epochs=1000)solution =solver.get_solution(best=True)

   该解决方案需要三个输入 - t 、 λ和U0 。所有输入必须具有相同的形状。例如，如果您有兴趣固定λ=4和U0=0.4并绘制解u与t ∈ [0,1]关系，您可以执行以下操作

   t = np.linspace(0, 1)lmd = 4 * np.ones_like(t)u0 = 0.4 * np.ones_like(t)u = 解决方案(t, lmd, u0, to_numpy=True)导入 matplotlib.pyplot 作为 pltplt .plot(t, u)

5. 一旦有了捆绑solution ，您就可以找到一组与观察到的数据点(t_i, u_i)最匹配的参数(λ, U0) 。这是使用简单的梯度下降来实现的。在下面的玩具示例中，我们假设只有三个数据点u(0.2) = 0.273 、 u(0.5)=0.129和u(0.8) = 0.0609 。以下是经典的 PyTorch 工作流程。

    # 观测数据点st_obs = torch.tensor([0.2, 0.5, 0.8]).reshape(-1, 1)u_obs = torch.tensor([0.273, 0.129, 0.0609]).reshape(-1, 1)#随机初始化λ 和 U0；跟踪它们的梯度lmd_tensor = torch.rand(1) * (LAMBDA_MAX - LAMBDA_MIN) + LAMBDA_MINu0_tensor = torch.rand(1) * (U0_MAX - U0_MIN) + U0_MINadam = torch.optim.Adam([lmd_tensor.requires_grad_(True), u0_tensor.requires_grad_(True)], lr=1e-2)# 运行梯度下降 10000 epochsfor _ in range(10000):output = Solution(t_obs, lmd_tensor * torch.ones_like(t_obs), u0_tensor * torch.ones_like(t_obs))loss = ((output - u_obs) ** 2).mean()loss.backward()adam.step()adam.zero_grad()
      print(f"λ = {lmd_tensor.item()}, U0={u0_tensor.item()}, 损失 = {loss.item()}")

常问问题

Q：如何使用GPU进行训练？

简单的。导入 Neurodiffeq 时，库会自动检测 CUDA 在您的计算机上是否可用。由于该库基于 PyTorch，因此如果找到兼容的 GPU 设备，它将默认张量类型设置为torch.cuda.DoubleTensor 。

问：如何使用预训练网络？

请参阅自定义网络和迁移学习部分。

问：如何改变学习率？

标准 PyTorch 方式。

1. 按照自定义网络中的说明构建网络： nets = [FCNN(), FCN(), ...]

2. 实例化一个自定义优化器并将这些网络的所有参数传递给它

   parameters = [p for net in nets for p in net.parameters()] # 所有网络的参数列表MY_LEARNING_RATE = 5e-3optimizer = torch.optim.Adam(parameters, lr=MY_LEARNING_RATE, ...)

3. 将nets和optimizer传递给求解器： solver = Solver1D(..., nets=nets, optimizer=optimizer)

问：我得到了一个糟糕的解决方案。

与传统的数值方法（FEM、FVM 等）不同，基于神经网络的解决方案需要进行一些超级调整。该库提供了最大的灵活性来尝试任何超参数组合。

• 要使用不同的网络架构，您可以传入自定义的torch.nn.Module 。

• 要使用不同的优化器，您可以将自己的优化器传递给solver = Solver(..., optimizer=my_optim) 。

• 要使用不同的采样分布，您可以使用内置生成器或从头开始编写自己的生成器。

• 要使用不同的采样大小，您可以调整生成器或更改solver = Solver(..., n_batches_train) 。

• 要在训练期间动态更改超参数，请查看我们的回调功能。

问：有什么经验法则吗？

• 不要使用ReLU进行激活，因为它的二阶导数相同为 0。

• 以无量纲形式重新调整 PDE/ODE，最好使所有内容都在[0,1]范围内。使用像[0,1000000]这样的域很容易失败，因为a) PyTorch 将模块权重初始化得相对较小， b)大多数激活函数（如 Sigmoid、Tanh、Swish）在 0 附近最非线性。

• 如果您的 PDE/ODE 太复杂，请考虑尝试课程学习。开始在较小的域上训练网络，然后逐渐扩展，直到覆盖整个域。

贡献

欢迎大家为这个项目做出贡献。

在向该存储库做出贡献时，我们考虑以下流程：

1. 打开一个问题来讨论您计划进行的更改。

2. 仔细阅读贡献指南。

3. 如果对界面进行了更改，请在分叉存储库上进行更改并更新 README.md。

4. 打开拉取请求。