Plotting the solution of a differential equation

New to Julia
1

I have written the Julia example for solving the Lorenz equation as in the page:
https://diffeq.sciml.ai/stable/tutorials/ode_example/#Example-2:-Solving-Systems-of-Equations :

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

u0 = [1.0;0.0;0.0]
tspan = (0.0,100.0)
prob = ODEProblem(lorenz!,u0,tspan)
sol = solve(prob)

The problem appears when plotting because I obtain two different results when I was expecting the same result. Below I have attached the piece of code together with the plots.
In the code, x[1,:],x[2,:] and x[3,:] are respectively the x,y,z components of the solution points.

julia_plot
julia_plot890Ă—1570 231 KB

Thank you very much in advance.

1 Like
2

When plotting DifferentialEquations.jl solutions, data interpolation takes place within Plots. It seems to be a built in recipe…

3

Welcome with an interesting problem.

It is a little difficult to read your problem. It is recommended to insert in-line code between two single back-ticks (`). If you want to typeset more lengthy, multiline code, start the code section with three back-ticks followed by the word julia on the same line. Then new line(s) with your code. Close the code section with three back-ticks (without the word julia).

Presenting code as screen shot makes it virtually impossible to read on my screen – simply because my browser doesn’t allow for zooming the image.

4

The plot recipe is using the continuous solution sol(t) to sample uniformly spaced points and generate a nicer plot than if you used the sample points sol.t and sol[i].

5

@ChrisRackauckas, in your article it is said that the continuous sol(t) is the interpolated solution at time t. Could you shed some light about this. Is it sampled or interpolated? Thank you.

6

Chris is of course the authority here. My understanding is that the solution, assume it is named sol, is a data structure where sol.t contains the vector of the time points where the solution is computed, while sol[i] contains the solution at time sol.t[i].

The data structure also contains an interpolation function, which is accessed using the idea of a “functor”: the “functor” facility enables the creation of a function with the same name as the solution data structure (sol in this example). So sol(t) calls this interpolation function for the solution at time t.

This means that if you use the interpolation function in exactly the time points where the solution has been computed, i.e., sol(sol.t[i]), this should be identical to sol[i] (with a possible reorganization of the structure of the solution). But whereas you can only access the solution data structure sol[i] at integer values i, the interpolation function can be assessed at any floating point value (note: I have not tested it for extrapolation).

Specifically, you can create a range of time points, e.g., T = range(0,20,length=500) (just as a stupid example), and then compute the solution via the interpolation function in these time points: sol.(T), and then (using Plots), plot the result by plot(T,sol.(T)).

I’m sure my explanation is not 100% accurate, but it is relatively close I think.

So, in summary: accessing the solution by square bracket, sol[i], accesses the computed values stored in an array, while accessing the solution by parenthesis, sol(t), uses function notation and accesses the interpolation function.

2 Likes
7

… a similar idea of “functor” is used with package Polynomials.jl:

julia> using Polynomials

julia> pol = Polynomial([1,-2,0.5])
Polynomial(1.0 - 2.0*x + 0.5*x^2)

julia> length(pol)
3

julia> pol[0]     # here, I access pol as a vector with the zero order coefficient, etc.
1.0

julia> pol[1]
-2.0

julia> pol[2]
0.5

julia> pol(3.14)     # here, I access pol as a *function* 
-0.35019999999999984
2 Likes
8

Thank you very much for all the answers. I am interested in reading some documentation about this issue. I have searched in both the DifferentialEquations.jl and Plots.jl documentation without success. I would like to see how is this functor built and how is the interpolation in the plots side performed. I guess there are some default arguments, and some arguments to be modified in order to manage the interpolation.
All the best.

9

All of the code is open and fairly straightforward. Here’s the plot recipe:

https://github.com/SciML/DiffEqBase.jl/blob/master/src/solutions/solution_interface.jl#L59

By default

https://github.com/SciML/DiffEqBase.jl/blob/master/src/solutions/solution_interface.jl#L66

it plots 10,000 evenly spaced points.

https://github.com/SciML/DiffEqBase.jl/blob/master/src/solutions/solution_interface.jl#L201-L208

It just uses sol(t). When I say “evenly”, if the plot is in log-space it takes those points in log space.

Is that what you were looking for?

3 Likes
10

In the Julia on-line documentation, search for “Function-like objects”.

11

@ChrisRackauckas @BLI :
Thank you very much for your answers. Please let me know if I have understood correctly according to the following statements:
The solution sol(t) is an interpolating function that can be evaluated at any “t”. Besides, inside the object sol, there is additional information (plot_density argument, a 10,000 points grid) that can be later used by Plots. The interpolation to this dense grid is not done by Plots, but it is already built in sol(t). Is this right?
Thanks again.

12

First – if I write something that is at odds with what Chris writes, you should believe in Chris and not me – I’m just a user; Chris and his group developed the package.

I don’t know what plot density is used when plotting with the command plot(sol) and its variants. But you can always create your own sampling of the function by, say, time = range(t1,t2,length=N) and then plot(time,sol.(time)) – here, you can choose t1, t2, and N relatively freely, although I’d assume you have to restrict t1 and t2 to be within the time span used when solving the model :slight_smile: .

13

It’s done in the plot recipe.

14

Problem understood. Thanks.