Search code examples
juliadifferential-equations

Julia's DifferentialEquations step size control

I want to solve the double pendulum equations using DifferentialEquations in Julia. For some initial values I get the error:

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

If I use force_dtmin=true, I get:

WARNING: Instability detected. Aborting

I do not know what to do further.Here is the code:

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))
Solution

  • I recently changed this warning to instead explicitly tell the user that it's most likely a problem with the model. If you see this, then there's often two possible issues:

    1. The ODE is stiff and you're using an integrator only for non-stiff equations
    2. Your model code is incorrect.

    While (1) used to show up more often, these days the automatic algorithm will auto-detect it, so the problem is almost always (2).

    So what you can do is print out what your calculated derivatives are and see if it lines up with what you were expected. If you did this, then you'd notice that

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

    is giving you dummy values which can be infinitely small/large. The reason is because you were supposed to overwrite them! So it looks like what you really wanted to do is calculate the first two derivative terms and use them in the second set of derivative terms. To do that, you have to make sure you update the values first! One possible code looks like this:

    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

    makes:

    enter image description here