我需要解决许多非常大的一阶 ODE 系统，它们描述了一些化学反应。变量的数量（在每个系统中）大约为. 我正在使用ALGLIB向量 ODE 求解器，它在后台实现Cash-Karp 自适应 ODE 求解器。
化学方程式通常被认为是僵硬的，建议使用隐式方法。不幸的是，隐式方法无法处理如此数量的变量，因为它们需要反转巨大的矩阵。这反过来又需要操作只是开始，那是在考虑舍入误差之前。上次当我“玩”一个订单较小的矩阵时（)，我需要大约 50 位数的精度才能获得合理正确的结果。当前处理器本身并不支持这种精度，并且所有这些都有效地排除了这种数量变量的隐式方法。
另一方面，自适应显式方法，如 ALGLIB 中使用的方法，计算误差并将步长减少到非常小的数字。这通常会使解决时间变得异常长（比如几个月）。
我的猜测是，当某些变量开始接近零时会发生这种情况。由于所有变量都必须是严格的非负数（它们是某些物质的浓度），一旦它们中的一些接近于零，那么相对误差就会立即变得非常大，算法就会减小步长。事实上，该算法通常会超调零，从而使一些变量为负数。在那之后，一切都会爆炸。因此，在计算导数时，我将所有小于零的变量视为精确零。这使得算法稳定，但仍然没有删除它可能会过多地减少步长的事实。
单个 ODE 系统的并行化不能很好地工作，或者，容我说，根本不能工作。这是因为 ODE 求解器反复调用一个相当快速的函数来计算导数。此类函数仅基于纯数学，并且对于此类函数（即使输入数组为大小）并行化开销只是扼杀了目的。相反，并行化是通过同时生成多个模型来实现的，其中每个模型在单个线程下运行。这可以通过产生多个线程或多个外部进程来完成（每个线程/进程运行一个“模型”）。两者都有一些优点和缺点，这与问题无关。底线是并行化单个模型会降低性能。
系统具有运动积分：系统中的总数或“原子”（所有变量的总和乘以某些权重）必须是常数。并且由于所有变量都是严格非负的，这导致了导数的上限。
我想知道如何才能更快地解决此类 ODE 系统。