DifferentialEquations can solve backward?

Using DifferentialEquations, I am solving a simple first order ODE backwards which should work, right ? But the solver fails. I can solve this using my own forward Euler scheme without any problem.

Note that the interval is (Rp, Rm) and Rp > Rm.

using DifferentialEquations

beta = 0.07
Rp = sqrt(2.0)
Rm = 0.08
Cf = 0.0036
g = 9.81
q0 = 1.2e-3
h0 = 0.003

# Radial velocity at r = Rp
u0 = -q0/(h0 * Rp)

# p is extra parameters: not used
function rhs(r, p, u)
   num = u * (  Cf * abs(u) * u^3 * r^3
              - u * g * q0 * r^2 * tan(beta)
              - q0^2 * g )
   den = q0 * r * (r * u^3 + q0 * g)
   return num/den
end

rspan = (Rp, Rm)
prob = ODEProblem(rhs, u0, rspan)
sol = solve(prob, Tsit5(), reltol=1e-8, abstol=1e-8

Did you mess up the order of arguments of the ODE right-hand side? It must be u, p, t for out-of-place problems (or du, u, p, t for in-place problems), see Ordinary Differential Equations · DifferentialEquations.jl. Here, u is the ODE solution (dependent variable), p are possible parameters, and t is the time (independent variable).

julia> using OrdinaryDiffEq # ODE part of DifferentialEquations

julia> f(u, params, time) = u
f (generic function with 1 method)

julia> ode_forward = ODEProblem(f, 1.0, (0.0, 1.0));

julia> sol_forward = solve(ode_forward, Tsit5());

julia> sol_forward.u[end], exp(1)
(2.718281708773342, 2.718281828459045)

julia> ode_backward = ODEProblem(f, 1.0, (0.0, -1.0));

julia> sol_backward = solve(ode_backward, Tsit5());

julia> sol_backward.u[end], exp(-1)
(0.3678795934275483, 0.36787944117144233)
3 Likes