3

I am trying to use the DifferentialEquations.jl provided by julia, and it's working all right until I try to use it on a second order ODE. Consider for instance the second order ODE

x''(t) = x'(t) + 2* x(t), with initial conditions

x'(0) = 0, x(0) = 1

which has an analytic solution given by: x(t) = 2/3 exp(-t) + 1/3 exp(2t).

To solve it numerically, I run the following code:

using DifferentialEquations;

function f_simple(ddu, du, u, p, t)
    ddu[1] = du[1] + 2*u[1] 
end;

du0 = [0.]
u0 = [1.]
tspan = (0.0,5.0)
prob2 = SecondOrderODEProblem(f_simple, du0, u0, tspan)
sol = solve(prob2,reltol=1e-8, abstol=1e-8);

With that,

sol(3)[2] = 122.57014434362732

whereas the analytic solution yields 134.50945587649028, and so I'm a bit lost here.

Trirac
  • 31
  • 2
  • 3
    It works if you specify a standard integrator like `solve(prob2, alg=Tsit5())`. It might be good to open an issue in the github repo with your example though, since the specialized integrators seem to be quite imprecise here. – Ilya Orson May 04 '21 at 00:27
  • This problem was actually addressed by the author last year, see the comments in https://stackoverflow.com/questions/60180865/2nd-order-odes-in-julia-using-differentialequations, https://github.com/SciML/DifferentialEquations.jl/issues/523. Please add the version of the diffeq package. – Lutz Lehmann May 08 '21 at 11:11
  • Open an issue please. – Chris Rackauckas May 08 '21 at 18:56
  • 1
    I've opened an issue on GitHub. Thank you all for your help, that has helped me a bunch already! – Trirac May 12 '21 at 19:23

1 Answers1

1

According to the the documentation for DifferentialEquations.jl, Vern7() is appropriate for high-accuracy solutions to non-stiff equations:

sol = solve(prob2, Vern7(), reltol=1e-8, abstol=1e-8)
julia> println(sol(3)[2])
134.5094558872943

On my machine, this matches the analytical solution quite closely. I'm not exactly sure what the default method used is: the documentation indicates that solve has some method of choosing an appropriate solver when one isn't specified.

For more information on Vern7(), check out Jim Verner's page on Runge-Kutta algorithms.

MattG
  • 1,416
  • 2
  • 13
  • 31