Order of arguments in Differential Equations functions

I am going through the documentation of DifferentialEquation and found this in the tutorial.

The function is defined as parameterized_lorenz(du,u,p,t) but when the arguments are supplied in a different order as prob = ODEProblem(parameterized_lorenz,u0,tspan,p). Is it a typo or it does not really matter in this case as internally this is taken care of?
I tried to change the order but it throws an error.

Defining Parameterized Functions

In many cases you may want to explicitly have parameters associated with your differential equations. This can be used by things like parameter estimation routines. In this case, you use the p values via the syntax:

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


We can make our functions look nicer by doing a few tricks. For example:

function parameterized_lorenz(du,u,p,t)
  x,y,z = u
  σ,ρ,β = p
  du[1] = dx = σ*(y-x)
  du[2] = dy = x*(ρ-z) - y
  du[3] = dz = x*y - β*z
end

and then we add the parameters to the ODEProblem:

u0 = [1.0,0.0,0.0]
tspan = (0.0,1.0)
p = [10.0,28.0,8/3]
prob = ODEProblem(parameterized_lorenz,u0,tspan,p)

The order of arguments of ODEProblem does not need to be related to the order of arguments of parameterized_lorenz. It could be anything, and that particular order was chosen.

@mauro3:
If that is the case, why I get an error when I swap arguments ‘p’ and ‘tspan’ in ODEProblem. Is it not that the third arugment has to be ‘tspan’? How the function parameterized_lorenz identifies various argument when supplied as part of ODE problem.
Similarly if I swap ‘p’ and ‘t’ in parameterized_lorenz function, I get an error.
This is confusing me

p is always last in ODEProblem because it’s optional.

Why this?

I think mauro3’s point was not that the order is irrelevant, just that whoever wrote the method did not have to put them in a particular order. Order still matters.

For example, if I write

raise(x, y) = x^y

And then call raise(2, 3) I will get a different answer than if I call raise(3, 2). Someone else might have written the method as

raise(x, y) = y^x

In which case the answers will be reversed. If your methods take different types (say the second argument is supposed to be an array) and you swap the order, you’ll get an error because you maybe passed an integer where the method was expecting an array.

Edit: lol, picked the one pair of numbers that actually have the same answer :derp:

This still does not answer my question.

It’s because the methods have a specifically defined argument ordering, like all methods do. The ODEProblem one was chosen that way because p is optional, since otherwise you’d always have to provide parameters and this is consistent with the rest of Julia (and most programming languages). The derivative definition has the independent variable last because PDE definitions can have more than one independent variable, so other definitions are like f(du,u,p,x...), which is consistent with how the rest of Julia (and most programming languages) handle optional sizing.

Now you just pointed out the p is at the end of the ODEProblem but not in the f, but as you can see here, if you think about all of the conventions that exist in programming languages, there is no single rule you can choose here that is compatible with all of them. So we went with the one that is compatible with the most rules and programming language features since that would require the least amount of hacky workarounds and thus result in the least number of bugs. It’s actually very nice that this is the one concession we had to make.

@ChrisRackauckas: Thanks.
Whatever little Julia I know, I understand a function can have positional arguments or keyword arguments. The function parameterized_lorenz is a user defined function and has no keyword arguments but only positional. Positional arguments have to be supplied in the same order. Am I correct?

I agree that for ODEProblem function, p is an optional argument and so it comes at last. No issues with that.

When parameterized_lorenz function is supplied to ODEProblem function, is it not that the lorenz function expecting the arguments by position else how it is able to identify that when argument are not being supplied in the same order as were defined when specifying the function. I think something is happening inside the code to accommodate this. I might be missing some critical aspect in this.

I wonder if your confusion is about what constitutes a keyword vs positional argument. In the code you posted up top, there are no keyword argument. Keyword arguments have the form arg=value inside the function call.

Consider:

raise(x,y) = x^y

x = 3
y = 2
raise(y, x)

Will you get 8 or 9?

Maybe the confusion comes from the following: By constructing the ODEProblem nothing is passed as parameters to parameterized_lorenz. An ODEProblem is constructed by passing the function.
Afterwards solve will use that function and call it, expecting that it has the signature (du,u,p,t).
So the signature of the function, which is the first parameter of the ODEProblem constructor is unrelated to the order of the rest of the parameters you pass to the ODEProblem constructor.

Did you observe that the function named parameterized_lorenz with arguments (du,u,p,t) is a totally different function than function ODEProblem with arguments (parameterized_lorenz,u0,tspan,p)? These two functions are completely different, and have different arguments. So you cannot use the order of arguments in one of these functions to say anything about the order of arguments in the other.

Specifically, in function ODEProblem, the first argument is the name of function parameterized_lorenz. But you can not draw the conclusion that the subsequent arguments of ODEProblem (i.e., u0,tspan,p) has any direct connection with the arguments of the parameterized_lorenz function.

In general, arguments that are positional have to stay in the given order. It is only for keyword arguments that you can change argument order.

@BLI and @BeastyBlacksmith: Thanks
You are correct that my confusion was that lorenz function is being specified using some arguments but nothing is being supplied as the value of those arguments. If I infer that function is only used to specify the ODE problem and its arguments has nothing to do with the calculation of the problem, am I correct in my understanding. I stripped all the arguments of lonrenz function and it still gives the same answer. If the lorenz function is not using any of the arguments, why it is being defined that way?

Am I missing anything here?

function parameterized_lorenz()
  du=zeros(3)
#   x,y,z = u
#   σ,ρ,β = p
  du[1] = dx = σ*(y-x)
  du[2] = dy = x*(ρ-z) - y
  du[3] = dz = x*y - β*z
  return du
end
u0 = [1.0,0.0,0.0]
t = (0.0,100.0)
p = [10.0,28.0,8/3]
prob = ODEProblem(parameterized_lorenz,u0,t,p)

sol = solve(prob)

The reson is probably the multible dispatch of parametrized_lorenz
The solver calls a former definition of the function. Just restart the julia session and you will see that the example woun’t work

Based on your questions, I would guess that (i) you are very new to Julia, and (ii) you have some but not much background in programming. That is ok – I’m not not a super programmer myself, and am trying to understand some of the intricacies of Julia. So: my answers are not the final answers – I may have misunderstood things myself. But let me try to answer some of your questions… from an amateur programmer:

  1. If you define a function which describes an ordinary differential equation (example: the parameterized_lorenz function above) and want to use the differential equation solvers in the DifferentialEquations.jl package, you have to use the specified function arguments. You can not just skip arguments. As MatFi suggest, the only reason why your code worked, is because of Julia’s “multiple dispatch” feature.
  • A simple-minded understanding of “multiple dispatch” (i.e., my understanding) is that it allows you to have several mappings with the same name, but with different behavior depending on the argument types. Because you already had defined the correct version with correct arguments in parameterized_lorenz (e.g., parameterized_lorenz(du,u,p,t), when you define a second mapping with the same name but zero arguments (parameterized_lorenz()) this does not replace the first mapping. Note that in, say, MATLAB or Python, this would replace the first function with the second one, but in Julia, instead, both mappings co-exist.
  • Next, when you apply the function, Julie chooses the one which has the same arguments as needed.

Just as a stupid illustration, suppose I make three versions of a mapping my_sin:

my_sin(x) = sin(x)
my_sin() = "Good evening"
my_sin(x::Vector) = sum(x);

The first and third versions have a single argument; the first one is chosen if the single argument is a scalar, while the third is chosen if the single argument is a vector. The second version is chosen if there is zero arguments. They all co-exist:

julia> my_sin(pi/2), my_sin(), my_sin([1,2,3])
(1.0, "Good evening", 6)

Similarly in your case: you first define function parameterized_lorenz with arguments that DifferentialEquations.jl can understand, next you define the same mapping with zero arguments – which doesn’t make sense to DifferentialEquations.jl. Thus, when you define the problem via function ODEProblem, function ODEProblem chooses to use the mapping with correct arguments. In other words, you haven’t really tested it with zero arguments. To test it with zero arguments, you need to close Julia, start Julia again and only define parameterized_lorenz with zero arguments – before you call ODEProblem. And then you will get an error message, because the correct arguments (du,u,p,t) are required, and you have not defined the function with those arguments.

  1. What about “function” ODEProblem? OK – here, I’m on “thin ice”. What I guess is happening is that the guys who developed DifferentialEquations.jl created a data structure named ODEProblem, a so-called struct in Julia speak. A struct is a way to store information which is user defined.

Suppose I define a struct named MyProblem, with 4 “fields”:

struct MyProblem
    func
    x1
    xN
    N
end

In Julia, the struct name automatically is the “constructor” for the struct. So I create a “problem” as follows:

prob = MyProblem(sin,-2pi,2pi,20)

I can read the field values of the data structure prob, e.g.:

julia> prob.x1
-6.283185307179586

I can now design a plotting function which operates on my MyProblem structure:

function my_plot(problem)
    x = range(problem.x1,problem.xN,length=problem.N)
    return plot(problem.func,x)
end

…and I can apply my function my_plot on the data structure prob:

my_plot(prob)

which leads to the result:

OK – isn’t this just a complicated way of doing:

julia> plot(sin,range(-2pi,2pi,length=20))

Yes it is!!

So why would one do such a thing? Well, in my stupid case: suppose I wanted to make several versions of the my_plot function, one where I draw straight lines between the data points, another where I use quadratic interpolation, one where I use step function, etc., etc? Then I could create the data structure prob and just change an option in the my_plot function, say add a second argument: my_plot(problem, spline_order=1).

Back to your Lorenz problem. What is going on is thus (I think):
A. You define the ordinary differential equation in the parameterized_lorenz. This function is different for every model, but must have the parameters (du,u,p,t) in the given order.

B. “Function” ODEProblem is really a “constructor” which pushes the argument values into fields in a data structure. In your case, the arguments of the constructor are (parameterized_lorenz,u0,tspan,p). The first of these arguments (parameterized_lorenz) is pushed into the first field of the data structure (which you have named prob), and by the definition of the ODEProblem data structure, this should be a function. The second argument is pushed into the field of the initial value. The third argument is pushed into the field of time span for the simulation, etc. Thus, the argument order in the constructor must be correct, or else the information is pushed to the wrong field in the data structure. Again, observe that the order of arguments which are pushed to fields in the constructor has nothing to do with the order of arguments in the function which defines the model.

So function ODEProblem only stores the information in a suitable data structure, nothing more.

C. Finally, the DifferentialEquations.jl package has defined a function solve with one required positional argument; the required argument is a data structure defining the model. When this argument is a data structure of type ODEProblem, function solve knows how to pull out information from the data structure and apply the ode solver.

So why on earth go through the complexity of storing the model in a data structure? One reason is that it structures the algorithms. Another is that the solve function has a second, optional positional argument which specifies the solver. So if I write:

julia> solve(prob)

the default solver is used. If I instead write:

julia> solve(prob,Tsit5())

the Tsit5() solver is used. And I don’t have to write a complex list of arguments – the “problem” prob is the same as before; I only change the solver.

If I want to use a fourth order Runge Kutta solver, I simply write:

julia> solve(prob,RK4())

I can also use the explicit Euler solver. However, this solver only works with a fixed step length, so I need to add a keyword argument with the step length:

julia> solve(prob,Euler(),dt=1e-3)

Anyway, this was an attempt to explain what is going on when you solve ODEs. Hope it clarifies things (and that I have a correct understanding…). Attempting to explain it gave me a chance to figure out structs… (although I only scratched the surface…).

6 Likes

A side note - would you please edit the title to something a bit more specific in the interest of future searching? Maybe something like “Order of arguments in Differential Equation functions”? Having a generic title makes it harder for this conversation to be useful to others (and I think it’s a useful discussion!)

1 Like

Thanks @BLI for giving a comprehensive answer. I am aware of multiple dispatch in Julia and missed that piece when removing the arguments from lorenz function and finding that function still works.

I need one more clarification wrt your above comment. I have probably reached root of the confusion. I have tested and find that lorenz function need to have the arguments in same order as (du,u,p,t). Even I swap one of the argument it results in error, which make me to assume that internally the ODE solver code expects a certain (given) order of arguments in lorenz function so it is able to identify correct pieces.

Am I correct to assume that?

1 Like

Yes, as Chris and others (including myself) have explained in prior responses:

  • If you expect your function to work with the DifferentialEquations packages, you have to follow the (positional) argument order specified in the documentation. And for ODEs, that order is (du,u,p,t). In other words, you have zero freedom to swap order of arguments. Changing to (du,u,t,p) or (u,du,p,t), etc. will not work – and is meaningsless.

This should not be surprising. As far as I know, every programming language where you use positional arguments in the function definition, and where you expect your function to work with some given tool, will require that you adhere to the specified order of arguments. Swapping of argument order is never accepted… swapping argument order makes it impossible for the tool using the function to understand what is going on.

OK – because of multiple dispatch in Julia, you can sometimes cheat on this (but not in general in other languages without multiple dispatch)…Suppose I make the following function repeat with general arguments, and add two methods with specialized arguments:

repeat(x,y) = x^y
repeat(x::String,y::Integer) = x^y
repeat(x::Integer,y::String) = y^x

Then, by using the repeat function, I get the following behavior:

julia> repeat(2,3)
8
julia> repeat("Hi ",3)
"Hi Hi Hi "
julia> repeat(2,"Hi ")
"Hi Hi "

In this case, you see that I am allowed to swap arguments. But that only works because:

  1. The arguments are of different type, thus Julia can figure out the meaning of the arguments, and furthermore
  2. I have done the extra work to add a method for both cases of argument order.

If you look at the arguments in the ODE model, (du,u,p,t), the normal thing would be that du and u are of the same type (say, vectors), and there is no way the code can figure out the order. Etc. Furthermore, there would be no gain in trying to introduce such flexibility in the case of ODEs.

In conclusion: the order of arguments when defining ODEs is fixed.

1 Like

@BLI. This helps to remove my confusion. Many Thanks to everyone who contributed in this discussion.