寻找对导数精度持宽松态度的 ODE 积分器/求解器
我有一个（一阶）ODE 系统，计算导数相当昂贵。
然而，在给定的误差范围内，导数的计算成本要低得多，因为导数是根据收敛级数计算的，并且边界可以放置在丢弃项的最大贡献上，或者通过使用存储在 kd 树中的预先计算的范围信息/八叉树查找表。
不幸的是，我还没有找到任何可以从中受益的通用 ODE 求解器； 他们似乎都只是给你坐标并想要返回准确的结果。 （请注意，我不是 ODE 方面的专家；我熟悉 Runge-Kutta、Numerical Recipies 书中的材料、LSODE 和 Gnu Scientific Library 的求解器）。
即对于我见过的所有求解器，您提供一个 derivs 回调函数，接受 t 和 x 数组，并返回一个数组dx/dt 返回； 但理想情况下，我正在寻找一个能够提供回调 t、xs、和一组可接受的错误，并接收 dx/dt_min 和 dx/dt_max 数组返回，导数范围保证在所需的精度范围内。 （可能有许多同样有用的变化）。
任何指向考虑到此类问题而设计的求解器的指针，或解决该问题的替代方法（我不敢相信我是第一个想要这样的东西的人）将不胜感激。
I have a system of (first order) ODEs with fairly expensive to compute derivatives.
However, the derivatives can be computed considerably cheaper to within given error bounds, either because the derivatives are computed from a convergent series and bounds can be placed on the maximum contribution from dropped terms, or through use of precomputed range information stored in kd-tree/octree lookup tables.
Unfortunately, I haven't been able to find any general ODE solvers which can benefit from this; they all seem to just give you coordinates and want an exact result back. (Mind you, I'm no expert on ODEs; I'm familiar with Runge-Kutta, the material in the Numerical Recipies book, LSODE and the Gnu Scientific Library's solver).
ie for all the solvers I've seen, you provide a derivs callback function accepting a t and an array of x, and returning an array of dx/dt back; but ideally I'm looking for one which gives the callback t, xs, and an array of acceptable errors, and receives dx/dt_min and dx/dt_max arrays back, with the derivative range guaranteed to be within the required precision. (There are probably numerous equally useful variations possible).
Any pointers to solvers which are designed with this sort of thing in mind, or alternative approaches to the problem (I can't believe I'm the first person wanting something like this) would be greatly appreciated.
如果你对这篇内容有疑问，欢迎到本站社区发帖提问 参与讨论，获取更多帮助，或者扫码二维码加入 Web 技术交流群。
绑定邮箱获取回复消息
由于您还没有绑定你的真实邮箱，如果其他用户或者作者回复了您的评论，将不能在第一时间通知您！
发布评论
评论（7）
粗略地说，如果您知道 f' 的绝对误差 eps，并从 x0 积分到 x1，则来自导数误差的积分误差将 <= eps*(x1 - x0)。 还有来自 ODE 求解器的离散化错误。 考虑 eps*(x1 - x0) 对您来说有多大，并向 ODE 求解器提供使用误差 <= eps 计算出的 f' 值。
Roughly speaking, if you know f' up to absolute error eps, and integrate from x0 to x1, the error of the integral coming from the error in the derivative is going to be <= eps*(x1 - x0). There is also discretization error, coming from your ODE solver. Consider how big eps*(x1 - x0) can be for you and feed the ODE solver with f' values computed with error <= eps.
我不确定这是一个恰当的问题。
在许多算法中，例如非线性方程求解，f(x) = 0，导数 f'(x) 的估计就可以用于牛顿法等算法，因为您只需要朝着答案的“大方向”前进。
然而，在这种情况下，导数是您正在求解的 (ODE) 方程的主要部分 - 如果导数错误，您将得到错误的答案； 这就像尝试仅使用 f(x) 的近似值来求解 f(x) = 0 一样。
正如另一个答案所建议的，如果您将 ODE 设置为应用 f(x) + g(x)，其中 g(x) 是一个误差项，您应该能够将导数中的错误与输入中的错误相关联。
I'm not sure this is a well-posed question.
In many algorithms, e.g, nonlinear equation solving, f(x) = 0, an estimate of a derivative f'(x) is all that's required for use in something like Newton's method since you only need to go in the "general direction" of the answer.
However, in this case, the derivative is a primary part of the (ODE) equation you're solving - get the derivative wrong, and you'll just get the wrong answer; it's like trying to solve f(x) = 0 with only an approximation for f(x).
As another answer has suggested, if you set up your ODE as applied f(x) + g(x) where g(x) is an error term, you should be able to relate errors in your derivatives to errors in your inputs.
经过更多思考后，我发现区间算术可能是关键。 我的 derivs 函数基本上返回间隔。 使用区间算术的积分器会将 x 维持为区间。 我感兴趣的是在最终的 t 处获得 x 上足够小的误差范围。 一个明显的方法是迭代地重新集成，提高每次迭代引入最大误差的样本质量，直到我们最终得到可接受范围内的结果（尽管这听起来可能是“比疾病更糟糕的治疗方法”）整体效率）。 我怀疑自适应步长控制可以很好地适应这种方案，选择步长以保持“隐式”离散化误差与“显式误差”（即间隔范围）相当。
无论如何，谷歌搜索“ode求解器间隔算术”或只是“间隔颂”会出现大量有趣的新的相关内容（VNODE 及其参考文献）。
Having thought about this some more, it occurred to me that interval arithmetic is probably key. My
derivsfunction basically returns intervals. An integrator using interval arithmetic would maintain x's as intervals. All I'm interested in is obtaining a sufficiently small error bound on thexs at a finalt. An obvious approach would be to iteratively re-integrate, improving the quality of the sample introducing the most error each iteration until we finally get a result with acceptable bounds (although that sounds like it could be a "cure worse than the disease" with regards to overall efficiency). I suspect adaptive step size control could fit in nicely in such a scheme, with step size chosen to keep the "implicit" discretization error comparable with the "explicit error" ie the interval range).Anyway, googling "ode solver interval arithmetic" or just "interval ode" turns up a load of interesting new and relevant stuff (VNODE and its references in particular).
如果您有一个刚性系统，您将使用某种形式的隐式方法，在这种情况下，导数仅在牛顿迭代中使用。 使用近似雅可比行列式将导致牛顿迭代中严格的二次收敛，但这通常是可以接受的。 或者（大多数情况下，如果系统很大），您可以使用无雅可比行列式的牛顿-克雷洛夫方法来求解阶段，在这种情况下，您的近似雅可比行列式仅成为预处理器，并且在牛顿迭代中保留二次收敛性。
If you have a stiff system, you will be using some form of implicit method in which case the derivatives are only used within the Newton iteration. Using an approximate Jacobian will cost you strict quadratic convergence on the Newton iterations, but that is often acceptable. Alternatively (mostly if the system is large) you can use a Jacobian-free Newton-Krylov method to solve the stages, in which case your approximate Jacobian becomes merely a preconditioner and you retain quadratic convergence in the Newton iteration.
您是否考虑过使用odeset？ 它允许您为 ODE 求解器设置选项，然后将选项结构作为第四个参数传递给您调用的求解器。 您可能最感兴趣的是误差控制属性（RelTol、AbsTol、NormControl）。 不确定这是否正是您需要的帮助，但这是我能想到的最好建议，我上次使用 MATLAB ODE 函数是在几年前。
另外：对于用户定义的导数函数，您是否可以将容差硬编码到导数的计算中，或者您真的需要从求解器传递误差限制吗？
Have you looked into using odeset? It allows you to set options for an ODE solver, then you pass the options structure as the fourth argument to whichever solver you call. The error control properties (RelTol, AbsTol, NormControl) may be of most interest to you. Not sure if this is exactly the sort of help you need, but it's the best suggestion I could come up with, having last used the MATLAB ODE functions years ago.
In addition: For the user-defined derivative function, could you just hard-code tolerances into the computation of the derivatives, or do you really need error limits to be passed from the solver?
不确定我贡献了多少，但在制药建模领域，我们使用 LSODE、DVERK 和 DGPADM。 DVERK 是一个很好的快速简单阶 5/6 Runge-Kutta 求解器。 DGPADM 是一个很好的矩阵指数求解器。 如果您的 ODE 是线性的，那么矩阵指数是迄今为止最好的。 但你的问题有点不同。
顺便说一句，T 参数只是为了一般性而存在。 我从未见过依赖于 T 的实际系统。
您可能正在进入新的理论领域。 祝你好运！
补充：如果你正在进行轨道模拟，在我看来，我听说过基于圆锥曲线的特殊方法。
Not sure I'm contributing much, but in the pharma modeling world, we use LSODE, DVERK, and DGPADM. DVERK is a nice fast simple order 5/6 Runge-Kutta solver. DGPADM is a good matrix-exponent solver. If your ODEs are linear, matrix exponent is best by far. But your problem is a little different.
BTW, the T argument is only in there for generality. I've never seen an actual system that depended on T.
You may be breaking into new theoretical territory. Good luck!
Added: If you're doing orbital simulations, seems to me I heard of special methods used for that, based on conic-section curves.
检查具有线性基函数和中点求积的有限元方法。 求解以下 ODE 只需要每个元素对 f(x)、k(x) 和 b(x) 各进行一次计算：
-k(x)u''(x) + b(x)u'(x) = f(x)
答案的逐点误差与您的评估误差成正比。
如果需要更平滑的结果，可以使用二次基函数，每个元素对上述每个函数进行 2 次评估。
Check into a finite element method with linear basis functions and midpoint quadrature. Solving the following ODE requires only one evaluation each of f(x), k(x), and b(x) per element:
-k(x)u''(x) + b(x)u'(x) = f(x)
The answer will have pointwise error proportional to the error in your evaluations.
If you need smoother results, you can use quadratic basis functions with 2 evaluation of each of the above functions per element.