随着讨论消除了原始帖子中的一些混淆，这里是迄今为止讨论的摘要：

我想实现一个隐式 ODE 求解器，但不知道当微分方程 (DE) 具有以下形式的不连续性时该怎么办：

1. 更常见的类型： $$\begin{align}\lim_{x\rightarrow 0}\frac{\sin{x}}{x}\end{align}$$
2. 不太常见，因为它只出现在一个 DE $$\begin{align}\min(\alpha x,\beta)\end{align}$$

正如基里尔所指出的，第一种“不连续性”是可移除的，在这种情况下，您可以使用挤压定理或 L'Hôpital 规则来证明这一点。如果没有更具体的信息，L'Hôpital 规则是解决这些类型限制的更通用策略之一，因此对于问题的其他部分，这是一个很好的启发式策略。

第二种“不连续性”实际上更多的是不可微分。

我正在尝试自己实现它，因为它是一组大规模耦合、相同的 DE，并且由于我们的管道，我们希望在 Matlab 的 GPU 上运行它们。

我不知道任何在 MATLAB 中的 GPU 上实现 ODE 求解器的库。SimEngine是一个 MIT 许可的 MATLAB 工具箱，似乎为 GPU 实现了 ODE 求解器。它似乎也不受支持（编写它的公司似乎不再存在，其网站也不再存在）。

除非您真的知道自己在做什么，否则我不建议您自己实施此类软件。正如我在评论中多次说过的，有 BSD 许可的库（例如 SUNDIALS），您可以根据自己的目的进行修改。

我考虑过调整/修改 Matlab ode15，但我仍然想要准确评估雅可比，因为系统有大约 10k 状态变量，所以对雅可比进行数值近似会很昂贵。

您不一定需要一个精确的雅可比矩阵，但您确实需要一个函数，该函数使用一些分析表达式为雅可比矩阵提供足够好的近似值。（W 方法不需要精确的雅可比矩阵，在实践中，许多库实现的方法只需要足够精确的雅可比矩阵。）您是正确的，通过有限差分评估雅可比矩阵会很昂贵。

我通常考虑在不连续性之后重新启动 ODE 求解器，但我不确定在实践中如何做到这一点：我是否在不连续性周围定义一个区域，然后简单地检查求解器是否到达它，然后重新启动？$\epsilon$

您有两个通用选项：

1. 忽略不连续性，在这种情况下，您的 ODE 求解器将是一阶精确的（对于有限的 ODE；在无限维情况下，误差将为，其中是一步的大小）。$O(h^{1/2})$$h$

2. 如果您想使用阶方法，请在您的右手边或其任何导数中找到任何不连续点，直到阶。在每个不连续处，重新启动积分器。通常，不连续性用连续可微的“事件函数”的根来表示。然后通过寻根定位不连续性。此功能在 SUNDIALS 中实现，Hairer、Nørsett 和 Wanner 可能在他们的教科书中讨论如何实现它。$k$$k - 1$