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:
- 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.
- 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âŚ).
