In-place functions in definition of ODEProblem (DifferentialEquations)

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

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

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

It should be. Iâ€™d gladly accept a PR.

1 Like

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.

This is a good explanation of stack- vs. heap-allocated variables: https://stackoverflow.com/questions/79923/what-and-where-are-the-stack-and-heap

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

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

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.

1 Like

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