我有一个由扩散反应物和中间体组成的系统模型。对于每个变量，，最终的代表方程看起来像标准的反应扩散形式：
在哪里是表示反应项的非线性函数。我正在使用有限差分法 (FDM) 用 Neumann 边界条件离散空间导数：
使用中心差分公式：
ode15s（基于 RK4 的求解器）及时对方程进行数值积分。
我在 PDE 上的工作不多，FDM 上的在线资源（甚至是维基百科）也描述了时间导数的同时离散化。
据我了解，这些案例中显示的离散化描述了简单的欧拉积分（显式或隐式）。甚至 Crank-Nicolson 方案也是这样显示的。
我的方法是否正确（离散空间导数并使用优化的 ode 求解器求解）？会不会有数值稳定性问题？
系统生物学家在这里研究反应扩散系统。您开始使用的方法可能是最常见的。事实上，我会说这是尝试解决此类问题的第一件事。如果只让为三对角矩阵 (1,-2,1)，则可以将系统写为：
其中是所有反应物的向量，是其向量化形式。然后可以将其放入任何 ODE 求解器中。这是 FancyPants 以简洁的形式指出的线法 (MOL) 方法。对于许多反应扩散问题，这已经足够了。如果您的不稳定性来自刚性反应方程（即非空间模型是刚性的），那么使用这种形式的刚性 ODE 求解器通常就足够了。不仅如此，由于这个方程的高度矢量化形式，您可以轻松地利用 GPGPU/Xeon Phi/多节点计算来“强力”解决方案（在 MATLAB 中，这需要很少/甚至不需要额外的努力） .
但是，如果您的扩散常数很高（或者您的扩散是一个高度可变的方程），您可能无法使系统足够稳定以使 MOL 方法起作用。如果是这种情况，您可以采取两种方式。一种常见的方法是 Crank-Nicholson 方法。一个快速的直觉是，它只是空间中的中点方法和时间中的中点方法。几乎任何计算 PDE 文本都会概述该方法。
但是，这将导致一个很大的隐式系统，您必须使用像牛顿法这样的非线性求解器来求解。真正的问题是你的隐式方程都是耦合的，这意味着它们不能并行求解，这可能会对大型系统造成伤害。因此，该领域最先进的方法之一是隐式积分因子方法。例如，如果您查看链接论文中的公式 26，如果您要将先前值的函数计算为某个大常数，然后您会看到隐式方程解耦为它在空间中的每个点都是独立的。这意味着您可以在空间中的每个点运行非线性求解器（意味着更小的雅可比行列式以获得更快的解决方案（并且您可以轻松地用笔/纸解决它而不是数值逼近它），并且您可以在空间中的每个点求解并行/GPU/等）并通过曲柄尼科尔森方法获得大量实际加速。
总结：首先尝试使用 ODE45 的线条方法，如果这不起作用，请尝试使用刚性求解器。否则去 Crank-Nicholson （如果你已经解决了那个代码），但如果你要解决很多这类问题，试试隐式积分因子方法。
注意：如果您正在查看周期性问题，您可能还想尝试Spectral 方法。另外，我假设您使用的是一维系统。如果不是，请查看算子拆分方法（ADI 或积分因子方法）或有限元方法。
您所做的称为线法，它是离散化演化 PDE 的最常见方法。Python中提供了如何稳定地解决反应扩散问题（由我编写）的示例。
正如你所说，你只能离散空间导数，从而获得一组 DAE（微分代数方程），或者以某种方式离散空间和时间以获得一组 AE（代数方程）。
正如评论中提到的，MOL（线法）是一种离散化方法，它假设仅离散空间域并将时间保留为计数。主要优点是求解器根据某些错误检查规则的执行自动选择时间步长。
如果您要实现一组 AE，从而离散化空间和时间域，那么由您来选择时间步长和/或实施能够根据某些错误检查规则接受/拒绝时间步长大小的算法。
此外，如果您要实现后一种公式，在某些情况下，数值方案可能会不稳定（在文献中有许多例子显示了这种特殊行为，同时用不同的数值方法以固定步长求解同一组方程）。最后，在使用非线性偏微分方程（如您正在使用的偏微分方程）时，对于 AE 实现，您还需要围绕工作点线性化非线性项，以实现某种牛顿-拉夫森方法来迭代求解非线性。无论如何，有很多方法可以做到这一点：Newman 教授开发了一种称为 BAND 的有效算法，该算法也可用于 Fortran 例程。这方面将由 DAE 的求解器自动管理，例如 ode15i、ode23、SUNDIALS 套件、DASSL 等。