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