In-place functions in definition of ODEProblem (DifferentialEquations)

Specific Domains Modelling & Simulations
1

I am struggling with something as basic as definition of a function for ODEProblem using DifferentialEquations package:

In Example 1 in the documentation for the package, the ODE is defined by a function on the right hand side of the ODE d/dt u(t) = f(u,p,t), that is

f(u,p,t) = 1.01*u

Then in the following Example 2 that deals with systems of ODEs, an in-place function is defined, namely

function lorenz(du,u,p,t)
 du[1] = 10.0*(u[2]-u[1])
 du[2] = u[1]*(28.0-u[3]) - u[2]
 du[3] = u[1]*u[2] - (8/3)*u[3]
end

In both cases, the ODEProblem is defined by including the name of the corresponding function among the input arguments to the ODEProblem constructor, that is, f and lorenz, respectively.

In the former (scalar) case, the function is a standard (output assigning) function while in the latter (multidimensional) case it is an in-place function, isn’t it? (By the way, shouldn’t then the latter function be named lorenz! and not just lorenz?)

Now, what if I define the problem from Example 1 using an in-place function also in the scalar case:

using DifferentialEquations
using Plots

function f!(du,u,p,t)
    du = 10.01*u
end

u0=1/2
tspan = (0.0,1.0)

prob = ODEProblem(f!,u0,tspan)
sol = solve(prob,Tsit5(),reltol=1e-8,abstol=1e-8)

plot(sol)

The above code returns an error ERROR: MethodError: no method matching similar(::Float64) after calling solve.

Thanks for any shared insight.

1 Like
2

The solver seems to assume that any in-place function will be array valued, and is trying to allocate a cache for itself to use internally while integrating. It calls similar on u0, but u0 is just a Float64, so it chokes. Changing the problem definition to

prob = ODEProblem(f!, [u0], tspan)

makes it run.

2 Likes
3

Exactly. Also, du = 1.01u isn’t actually in place since it’s not mutating anything on the heap, so that won’t work. Stack allocated variables thus need to be handled in a very different manner.

1 Like
4

It should be. I’d gladly accept a PR.

1 Like
5

Well, indeed, the solver no longer complains but the solution obtained in this way is certainly not correct (it is just constant).

Chris points out below that my f! function is actually not in-place. I confess at this point I am lost - I cannot see a difference with respect to the lorenc example (shown above). I guess this may have something to do with the difference between scalars and vectors used as function arguments.

6

This is a good explanation of stack- vs. heap-allocated variables: memory management - What and where are the stack and heap? - Stack Overflow

Basically, when du is a vector, it is a heap variable. When it is scalar, it a stack variable. The scalar du inside the function is created during the function call, and then forgotten when the function returns. It’s referring to a different location in memory than the du outside the function, so you can’t mutate it.*

du = 0.0
function foo!(u)
    du = 2.0 * u
    println("inside the function, du = $du")
end

println(du)
# 0.0
foo!(10)
# inside the function, du = 20.0
println(du)
# 0.0

*Actually you can, if you say global du = 2.0 * u in the function body. But there’s no good reason to, since allocating a single value on the stack will be as fast or faster than mutating a single-element array on the heap. Once an array gets big enough, though, it becomes faster to overwrite it in-place, rather than allocate a new block of memory for a new array.

1 Like
7

I just wanted to make sure about the in-place function definition here. Is the exclamation mark ! at the end of the function lorenz! is simply a convention or does it do something?

When I run

function lorenz(du,u,p,t)
 du[1] = 10.0*(u[2]-u[1])
 du[2] = u[1]*(28.0-u[3]) - u[2]
 du[3] = u[1]*u[2] - (8/3)*u[3]
end

and

function lorenz!(du,u,p,t)
 du[1] = 10.0*(u[2]-u[1])
 du[2] = u[1]*(28.0-u[3]) - u[2]
 du[3] = u[1]*u[2] - (8/3)*u[3]
end

I get the same results. I just wanted to make sure that this is how it should be.

Just wanted to mention that I saw the answer here in stack overflow about in-place functions, but I still wanted to make sure it is fine to use lorenz and lorenz!

8

Roi,

the exclamation mark (the bang) really does NOTHING. It is just another symbol, a part of the name. I really mean it, you can view it as a symbol “A” or “3”. The convention in Julia community (not always adhered to) is to include it in the function name if the function modifies the input arguments.

2 Likes
9

Thanks so much! I guess I just wanted to hear it from someone.

10

Yes, it’s a convention, 5th entry in in the style guide:
https://docs.julialang.org/en/v1/manual/style-guide/#bang-convention

Also confused me! (no pun intended :slight_smile: )