Julia:如何调试微分方程求解问题

这篇文章是 Chris Rackauckas 的帖子的翻译和总结,但是并不按照原文完全翻译,有个人的取舍,可以的话建议观看原文。原文地址是:PSA: How to help yourself debug differential equation solving issues

1/ 如何自己 debug 微分方程求解的问题

Debug 微分方程求解问题基本上都是在做同样的一件事,所以这篇文章是一篇总结。如果想要查看有关特殊问题的更多信息,可以去看 DifferentialEquations.jl 文档的 FAQ 一节

2/ 如何调试微分方程求解器发散?

dt <= dtmin. Aborting. There is either an error in your model specification or the true solution is unstable.

NaN dt detected. Likely a NaN value in the state, parameters, or derivative value caused this outcome.

Instability detected. Aborting

如果发现自己遇到了微分方程求解器出现了问题,第一件事情就是降低精度,或者说降低容忍度(tolerance)。很多情况下,降低容忍度可以提高稳定性。可以试试让 abstol=1e-10,reltol=1e-10 看看是不是容忍度高的问题。

如果这个方法没有用,试试使用更稳定的求解器,例如一些对于刚性方程(stiff equations)的求解器:TRBDF2(),KenCarp4() 或者 QNDF().

如果都没有用,问题可能出现在模型上。如果怀疑 Julia 求解器有问题,有一个办法去证明:使用非 Julia 的求解器。只要简单地把 solve(prob,Tsit5()) 改为 solve(prob,CVODE_BDF()) 就会使用经典 Sundials C/C++ 库。

  • Sundials.jl,C/C++ SUNDIALS 库的封装,包括 CVODE_Adams, CVODE_BDF, IDAARKODE.
  • ODEInterfaceDiffEq.jl,经典 Hairer Fortran 代码的封装,例如 dorpi5, dop853, radau, rodas 等等.
  • LSODE.jl,经典 lsoda 算法的封装.
  • MATLABDiffEq.jl,MATLAB ODE 求解器 ode45, ode15s 等的封装.
  • SciPyDiffEq.jl,SciPy 的 odeint (LSODA) 和其它方法(LSODE 等)的封装.
  • deSolveDiffEq.jl,R 语言库常用方法的封装.

注意:这些求解器包没有默认安装,在使用之前需要先安装包,例如在使用 Sundials 求解器之前通过 ]add Sundials; using Sundials 来安装先

如果你的模型在所有主流求解器上都失败了,包括从 C/C++ 和 Fortran 调用的所有主流求解器,那么问题不在于求解器,而在于你的模型。用所有语言创建的每个求解器都是不正确的,而不是你几个小时前编写的代码的可能性非常小。

以下是常见问题列表:

  • 仔细逐项地检查自己的模型。看看模型里面是否有哪一些项会无限增大的,哪一项的导数会变得非常大,为什么变大了,变大是不是正常的。
  • 仔细检查模型的假设。记住导数并不一定随着 u 变为零而变为负数。u' = -sqrt(u) 在有限时间内达到零,仅仅只是正在建模的系统有一个属性(为正数),但是并不意味着模型实际上在求解的时候也具有这个属性。可以去查一查导致与该属性相反的项,看看导数的值是不是正确。
  • 仔细检查是不是违反了 ODE 假设。ODE 右侧的 f 函数应该始终提供相同的结果,即 u' = f(u,p,t) 需要唯一定义,否则一定无法求解。
    • 如果 f 函数出现了随机性,求解器的自适应性会认为 ODE 正在以高的错误率求解(因为导数不断变化),为了让随机性降低到零,这样就会达到 dtmin。如果确实需要随机性,用 SDE 或者 RODE 求解器);
    • 如果 f 函数会修改 u,那么用不同的步长调用 f 会是不确定的,这样也会导致求解失败。如果确实需要这么做,可以使用 callback;
    • 如果 f 函数缓存上一步的值,意味着如果改变了 dt,那么就是在改变 f,这样 u' 也不再被定义了。自适应性 ODE 求解器不一定固定地往前求解,有可能会先尝试一个大的步长,然后再选择小的步长;

3/ 性能表现的问题

以下的一些网址是在讲如何以最好的表现去求解微分方程:

Solving Stiff Equations

Optimizing DiffEq Code

Optimizing Serial Code

Optimizing Serial Code in Julia 1: Memory Models, Mutation, and Vectorization

如果你已经阅读了这些教程,但仍然有性能问题,或者有清晰的问题要问,可以去给 ChrisRackauckas 提问,可以选择在 GitHub 上的 DifferentialEquations.jl 上提交 issue。但是在提问之前请查看这些教程,因为其中涵盖了人们需要的大部分内容!

智能之路 文章被收录于专栏

包括机器学习、神经网络、深度学习、强化学习各种方面的文章

全部评论

相关推荐

点赞 评论 收藏
分享
评论
点赞
收藏
分享

创作者周榜

更多
牛客网
牛客网在线编程
牛客网题解
牛客企业服务