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

简介: 这篇文章是 Chris Rackauckas 的帖子的翻译和总结,但是并不按照原文完全翻译,有个人的取舍。
这篇文章是 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。但是在提问之前请查看这些教程,因为其中涵盖了人们需要的大部分内容!

目录
相关文章
|
网络架构
子网划分中subnet-id可以取全0和全1吗?子网计算实战
子网划分划分中的全0 和全 1在不同模式下处理情况不同。分为 classful 和classless,如果你的路由器工作在classful环境下,全0 和全1网段是不能使用的,而classless的掩码任何时候都和IP地址成对地出现。所以说要看题目给的具体情况,
808 0
|
机器学习/深度学习 数据可视化 PyTorch
【PyTorch】TensorBoard基本使用
【PyTorch】TensorBoard基本使用
1031 0
|
前端开发 Java 测试技术
基于SSM的大学生电信诈骗宣传系统设计与实现
基于SSM的大学生电信诈骗宣传系统设计与实现
581 0
|
数据可视化 数据挖掘 Python
【数据分析与可视化】Matplotlib中动态rc参数设置详解与实战(图文解释 附源码)
【数据分析与可视化】Matplotlib中动态rc参数设置详解与实战(图文解释 附源码)
829 0
|
4月前
|
消息中间件 数据管理 Serverless
阿里云消息队列 Apache RocketMQ 创新论文入选顶会 ACM FSE 2025
阿里云消息团队基于 Apache RocketMQ 构建 Serverless 消息系统,适配多种主流消息协议(如 RabbitMQ、MQTT 和 Kafka),成功解决了传统中间件在可伸缩性、成本及元数据管理等方面的难题,并据此实现 ApsaraMQ 全系列产品 Serverless 化,助力企业提效降本。
|
6月前
|
消息中间件 存储 Kafka
分布式消息中间件设计与实现
本文深入探讨了消息中间件的核心功能实现与高并发、高可用设计。在生产者设计中,涵盖消息构造、序列化、路由策略及可靠性保障(如ACK机制)。消费者部分分析了拉取/推送模式、分区分配与消息确认机制。同时,Broker作为核心组件,负责消息路由、存储和投递,并通过索引技术实现快速检索。 高并发设计方面,重点讨论了文件存储(顺序写入、分段存储)、日志结构存储及负载均衡策略(如哈希分区、轮询分区)。为确保高可用性,文章详细解析了主从复制、故障转移机制以及同城/异地多活容灾方案。
|
设计模式 前端开发 C语言
【设计模式】 观察者模式介绍及C代码实现
观察者模式(Observer Pattern)是一种常用的设计模式,它定义了一种一对多的依赖关系,让多个观察者对象同时监听某一个主题对象,当主题对象发生变化时,它的所有观察者都会收到通知并更新自己的状态。观察者模式又称为发布-订阅模式。Subject(主题):被观察的对象,它将所有观察者对象的引用保存在一个集合中,并提供了添加和删除观察者对象的方法。Observer(观察者):观察者接口,定义了更新自己的状态的方法,以便主题在状态发生变化时通知观察者。ConcreteSubject(具体主题)
687 0
【设计模式】 观察者模式介绍及C代码实现
|
缓存 关系型数据库 MySQL
为啥MySQL官方不推荐使用uuid或者雪花id作为主键
为啥MySQL官方不推荐使用uuid或者雪花id作为主键
331 1
|
算法 vr&ar UED
硬核解决Sora的物理bug!美国四所顶尖高校联合发布:给视频生成器装个物理引擎
【5月更文挑战第16天】美国四所顶级高校联合推出PhysDreamer,将物理引擎集成到视频生成模型,以实现更真实的3D对象动态交互。该技术利用动态先验知识估计物体物理属性,生成逼真的动态视频。实验显示PhysDreamer在动态逼真度上超越现有方法,但在计算成本和处理复杂物理交互方面仍有局限。研究团队对未来持乐观态度,期待改善效率并扩展应用范围。这一创新将推动虚拟体验技术的发展,增强VR/AR的沉浸感和多领域应用。[论文链接](https://arxiv.org/pdf/2404.13026.pdf)
248 2
|
人工智能 编解码 搜索推荐
AI绘画入门:从小白到入门,轻松玩转AI作画
随着AI技术的不断发展,AI绘画已经不再是遥不可及的梦想,它正逐渐走入大众视野,成为了一种新兴的艺术创作形式。即使没有绘画基础,你也可以通过AI工具轻松创作出精美的作品。本文将带你从小白入门,学习AI绘画的基础知识和操作技巧,让你快速体验AI绘画的乐趣。
1001 0