Julia 的 DifferentialEquations 步长控制

Question

我想在 Julia 中使用 DifferentialEquations 求解 double pendulum 方程。对于某些初始值，我得到错误：

WARNING: dt <= dtmin. Aborting. If you would like to force continuation with
 dt=dtmin, set force_dtmin=true

如果我使用 force_dtmin=true，我得到：

WARNING: Instability detected. Aborting

不知如何是好further.Here是代码：

using DifferentialEquations
using Plots

m = 1
l = 0.3
g = pi*pi
function dbpen(du,u,pram,t)
  th1 = u[1]
  th2 = u[2]
  thdot1 = du[1]
  thdot2 = du[2]
  p1 =  u[3]
  p2 =  u[4]  
 du[1] = (6/(m*l^2))*(2*p1-3*p2*cos(th1-th2))/(16-9*(cos(th1-th2))^2)
 du[2] = (6/(m*l^2))*(8*p2-3*p1*cos(th1-th2))/(16-9*(cos(th1-th2))^2)
 du[3] = (-0.5*m*l^2)*(thdot1*thdot2*sin(th1-th2)+(3*g/l)*sin(th1))
 du[4] = (-0.5*m*l^2)*(-thdot1*thdot2*sin(th1-th2)+(g/l)*sin(th2))
end

u0 = [0.051;0.0;0.0;0.0]
tspan = (0.0,100.0)
prob = ODEProblem(dbpen,u0,tspan)
sol = solve(prob)

plot(sol,vars=(0,1))

Answer 1

您的 ODE 似乎很僵硬，默认情况下需要极小的 dt。您可以切换到刚性 ODE 求解器或给出如下提示：

sol = solve(prob,alg_hints=[:stiff])

参考：包文档中的 ODE 示例

Answer 2

我最近更改了此警告，改为明确告诉用户这很可能是模型有问题。如果你看到这个，那么通常有两个可能的问题：

ODE 是刚性的，您只对非刚性方程使用积分器
您的型号代码不正确。

虽然 (1) 过去经常出现，但现在自动算法会自动检测它，所以问题几乎总是 (2)。

所以你可以做的是打印出你计算的导数是什么，看看它是否符合你的预期。如果你这样做，那么你会注意到

thdot1 = du[1]
thdot2 = du[2]

为您提供可以无限 small/large 的虚拟值。原因是因为你应该覆盖它们！所以看起来你真正想做的是计算前两个导数项并在第二组导数项中使用它们。为此，您必须确保首先更新值！一种可能的代码如下所示：

function dbpen(du,u,pram,t)
  th1 = u[1]
  th2 = u[2]
  p1 =  u[3]
  p2 =  u[4]  
  du[1] = (6/(m*l^2))*(2*p1-3*p2*cos(th1-th2))/(16-9*(cos(th1-th2))^2)
  du[2] = (6/(m*l^2))*(8*p2-3*p1*cos(th1-th2))/(16-9*(cos(th1-th2))^2)
  thdot1 = du[1]
  thdot2 = du[2]
  du[3] = (-0.5*m*l^2)*(thdot1*thdot2*sin(th1-th2)+(3*g/l)*sin(th1))
  du[4] = (-0.5*m*l^2)*(-thdot1*thdot2*sin(th1-th2)+(g/l)*sin(th2))
end

品牌：

Julia 的 DifferentialEquations 步长控制

Julia's DifferentialEquations step size control

differential-equations

julia