Stepsize control in ODEProblem

Specific Domains Numerics
1

Hello,

I have a system of ODE, which is rather heterogeneous. That is one equation corresponds to some background property of a system, and then there are many (say, hundred) equations for individual evolution of components of the system. Thus, error estimate is quite different for the various equations - first equation is large and basically has no error for any time step, while the rest are actually important to determine the error.

So, I’d like to produce my own error estimate or adaptive step size control. What would be the suggestions? Actually, in my particular case, I don’t need something very smart – I more or less need to estimate the errors only for part of the equations I am solving.

Attempt to use the internalnorm does not help – as far as the error estimate internally uses maximum of the u[i], it anyway scales the error estimate in the same way for all equations – that is not what I need.

Any suggestions?

1 Like
2

internalnorm lets you choose any function of u that you want for the error, and it doesn’t default to the maximum of the u[i] so I have no idea where you’re getting your information from (is there an error in the docs somewhere? If so, please link). If you’re talking about the error calculation itself, it uses the maximum of uprev[i] vs u[i] to determine a scale for the ith component’s relative tolerance, but in no way is that equivalent to a maximum over all of the u.

For norm choices, please refer to the documentation here:

http://docs.juliadiffeq.org/latest/basics/common_solver_opts.html#Advanced-Adaptive-Stepsize-Control-1

For example, if you want the error norm to be l2 on the first 5 variables and maximum norm on the rest, you’d do:

internalnorm = (u)-> norm(u[1:5],2) + norm(u[6:end],Inf)`

and of course make that non-allocating with views or what not.

However, most cases are probably covered with a much simpler approach: vector tolerances. If you do something like

reltol = [1,1,1,1,1,1e-6,1e-6,1e-6]

then the first 5 components will have a much higher error tolerance.

3

I think I found what description might’ve tripped you up and edited it:

https://github.com/JuliaDiffEq/DiffEqDocs.jl/blob/master/docs/src/extras/timestepping.md#common-setup

4

I was lost mostly because the form that you suggest for internalnorm does not really work. In the current implementation the internalnorm gets evaluated both for full vector arguments, and for some scalar values… Here is (a rather silly) example:

using DifferentialEquations

function dfdt(u, p, t)
    return [ u[1], -u[2] ]
end

u0 = [1.0, 1.0]
tspan = (0.0, 10.0)
prob = ODEProblem(dfdt, u0, tspan)
sol = solve(prob, internalnorm=(u)->(u[1]^2+10*u[2]^2))

It throws an error.

The vector tolerances - interesting, I was not aware that this is possible, I’ll try.

5

Exactly, this part of the docs required a bit of guessing on what was actually going on. However, seems the real implementation is more subtle…

6

Ok, as a followup, vector reltol does not work (at least for the released version of DifferentialEquations, I have not checked for the current master version)

sol = solve(prob, reltol=[1e-6, 1e-10])

It fails with:

ERROR: MethodError: no method matching isless(::Float64, ::Array{Float64,1})
7

Oh, it’s just with the polyalgorithm I think. If you choose something like Tsit5 directly then it should work. I’ll get that fixed up tonight.

8

It fails also for Tsit5, RK4. Though the error message is different

ERROR: DimensionMismatch("Cannot left-divide transposed vector by matrix")
9

It sounds like you have a stiff system and should be using a stiff solver?

10

It is a stiff system, yes. And yes, I will check a stiff solver – however, as a starter I was interested to understand what is going on – and actually in my case I know exactly, what is going on – as far as the system is heterogeneous (several equations describe universe background expansion, and the majority of equations are Boltzman equations for fequency components), I know that I do not want to treat the errors in these parts in a similar way. And yes, if I increase precision only on the part of the equations, even the good old RK4 works perfectly and fast (checked using C program…).

Maybe, a very clever stiff solver can deal with this – but looking at what people done before, this is probably not the case.

11

I took a look at what’s going on and here’s some working scripts:

using OrdinaryDiffEq

function dfdt(u, p, t)
    [u[1],-u[2]]
end

function dfdt(du, u, p, t)
    du[1] = u[1]
    du[2] = -u[2]
end

u0 = [1.0, 1.0]
tspan = (0.0, 10.0)
prob = ODEProblem(dfdt, u0, tspan)
function internalnorm(u)
    u[1]^2+10*u[2]^2
end
internalnorm(u::Number) = norm(u)
sol = solve(prob, Tsit5(), internalnorm=internalnorm)
sol = solve(prob, Tsit5(), reltol = [1,1e-6])
sol = solve(prob, Tsit5(), reltol = [1,1e-6], internalnorm=internalnorm)
prob = ODEProblem{false}(dfdt, u0, tspan)
sol = solve(prob, Tsit5(), internalnorm=internalnorm)

# Fails
sol = solve(prob, Tsit5(), reltol = [1,1e-6], internalnorm=internalnorm)

One thing that was added recently is that you need internalnorm to also give you an appropriate norm for the element type. This was added recently because it’s required for Vector{SArray} support, but it seems to have not been documented. When that dispatch is added, everything works except when you have vector tolerances for the out-of-place form. Please feel free to open an issue there, but in most cases (especially your case with 100 equations) you should be using the in-place form for efficiency anyways so it’s really not a big deal.

You can also just solve(prob,RK4(),adaptive=false,dt=...) if necessary. RK4 has an interesting error estimator mostly for its global error control, but you can always turn these off (for the native Julia codes).

But yeah, we should update the internalnorm docs.

1 Like
12

Thank you! Works perfectly. Actually, seems that it works for the out-of-place form also, and for both Tsil5 and RK4.

By the way, while converting my function to the in-place form, I noticed a nice “feature”. I had a rescaling of the equations at the very end, and noticed that

function dfdt(du, u, p, t)
    du[1] = u[1]
    du[2] = -u[2]
    du[1:end] *= 2
end

works perfectly, while

function dfdt(du, u, p, t)
    du[1] = u[1]
    du[2] = -u[2]
    du *= 2
end

creats a new copy of du instead of modifying it in place. Seems quite a nice catch! Though this is a Julia “feature”, I think a warning in the place describing the definition of the in-place derivative functions for ODEs may be useful.

(I know, I probably should have made the multiplication while I was filling the derivatives, but that is not the point :slight_smile: )

As for the last comment – no, I do not want the constant time step! I definitely want an adaptive algorithm.