# Stepsize control in ODEProblem

**URL:** <https://discourse.julialang.org/t/stepsize-control-in-odeproblem/9772>\
**Category:** Numerics\
**Tags:** diffeq\
**Created:** [March 17, 2018, 3:24pm UTC](https://discourse.julialang.org/t/stepsize-control-in-odeproblem/9772 "2018-03-17T15:24:11Z")\
**Posts on this page:** 12\
**Page:** 1

<div class="post-metadata">

**Author:** ![fedxa](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fedxa/32/488_2.png) [@fedxa](https://discourse.julialang.org/u/fedxa)\
**Post date:** [March 17, 2018, 3:24pm UTC](https://discourse.julialang.org/t/stepsize-control-in-odeproblem/9772/1 "2018-03-17T15:24:11Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [March 17, 2018, 6:24pm UTC](https://discourse.julialang.org/t/stepsize-control-in-odeproblem/9772/2 "2018-03-17T18:24:40Z")

</div>

> [@fedxa](#):
>
> 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.

`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 `i`th 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](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:

```julia
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

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

```

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

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [March 17, 2018, 6:31pm UTC](https://discourse.julialang.org/t/stepsize-control-in-odeproblem/9772/3 "2018-03-17T18:31:33Z")

</div>

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](https://github.com/JuliaDiffEq/DiffEqDocs.jl/blob/master/docs/src/extras/timestepping.md#common-setup)

---

<div class="post-metadata">

**Author:** ![fedxa](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fedxa/32/488_2.png) [@fedxa](https://discourse.julialang.org/u/fedxa)\
**Post date:** [March 17, 2018, 7:24pm UTC](https://discourse.julialang.org/t/stepsize-control-in-odeproblem/9772/4 "2018-03-17T19:24:16Z")

</div>

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:

```julia
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.

---

<div class="post-metadata">

**Author:** ![fedxa](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fedxa/32/488_2.png) [@fedxa](https://discourse.julialang.org/u/fedxa)\
**Post date:** [March 17, 2018, 7:24pm UTC](https://discourse.julialang.org/t/stepsize-control-in-odeproblem/9772/5 "2018-03-17T19:24:21Z")

</div>

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…

---

<div class="post-metadata">

**Author:** ![fedxa](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fedxa/32/488_2.png) [@fedxa](https://discourse.julialang.org/u/fedxa)\
**Post date:** [March 17, 2018, 7:24pm UTC](https://discourse.julialang.org/t/stepsize-control-in-odeproblem/9772/6 "2018-03-17T19:24:22Z")

</div>

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)

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

```

It fails with:

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

```

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [March 17, 2018, 7:42pm UTC](https://discourse.julialang.org/t/stepsize-control-in-odeproblem/9772/7 "2018-03-17T19:42:02Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![fedxa](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fedxa/32/488_2.png) [@fedxa](https://discourse.julialang.org/u/fedxa)\
**Post date:** [March 17, 2018, 7:45pm UTC](https://discourse.julialang.org/t/stepsize-control-in-odeproblem/9772/8 "2018-03-17T19:45:32Z")

</div>

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

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

```

---

<div class="post-metadata">

**Author:** ![dpsanders](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dpsanders/32/3573_2.png) [@dpsanders](https://discourse.julialang.org/u/dpsanders)\
**Post date:** [March 18, 2018, 5:48pm UTC](https://discourse.julialang.org/t/stepsize-control-in-odeproblem/9772/9 "2018-03-18T17:48:34Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![fedxa](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fedxa/32/488_2.png) [@fedxa](https://discourse.julialang.org/u/fedxa)\
**Post date:** [March 18, 2018, 6:19pm UTC](https://discourse.julialang.org/t/stepsize-control-in-odeproblem/9772/10 "2018-03-18T18:19:34Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [March 18, 2018, 8:07pm UTC](https://discourse.julialang.org/t/stepsize-control-in-odeproblem/9772/11 "2018-03-18T20:07:02Z")

</div>

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

```julia
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.

> [@fedxa](#):
>
> 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…).

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.

---

<div class="post-metadata">

**Author:** ![fedxa](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fedxa/32/488_2.png) [@fedxa](https://discourse.julialang.org/u/fedxa)\
**Post date:** [March 19, 2018, 1:31am UTC](https://discourse.julialang.org/t/stepsize-control-in-odeproblem/9772/12 "2018-03-19T01:31:29Z")

</div>

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

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

```

works perfectly, while

```julia
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 🙂 )

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