Second order ODE in Julia giving wrong results

368 Views Asked by At

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.

1

There are 1 best solutions below

1
On

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.