Extracting time series from DifferentialEquations solution


#1

I’m exploring DifferentialEquations,jl (and teaching it to high school students tomorrow :slight_smile:). First demo is simple projectile motion, a 4d ODE wih state vector u = [x, y, vx, vy]. I want to extract x and y from the solution and plot them together. But the code doesn’t behave as the documentation suggests (http://docs.juliadiffeq.org/latest/basics/solution.html).

Minimal example

using DifferentialEquations

# frictionless 2d projectile motion
#  u = [x, y, vx, vy], f(u) = du/dt = [vx, vy,  0, -9.8]
function f(t, u)
    return [u[3], u[4], 0.0, -9.8]
end

x0 = 0.0   # initial horizontal position
y0 = 1.0   # initial height
v0 = 10.0  # initial speed
θ  = 0.5   # inital angle

u0 = [x0, y0, v0*cos(θ), v0*sin(θ)] 
tspan = (0.0, 5.0) 

prob = ODEProblem(f, u0, tspan)
soln = solve(prob)

Running this produces the sensible answer

julia> soln.u
7-element Array{Array{Float64,1},1}:
 [0.0, 1.0, 8.77583, 4.79426]           
 [0.00173129, 1.00095, 8.77583, 4.79232]
 [0.0190442, 1.01038, 8.77583, 4.77299] 
 [0.192173, 1.10263, 8.77583, 4.57966]  
 [1.92346, 1.8154, 8.77583, 2.64632]    
 [19.2363, -12.0343, 8.77583, -16.687]  
 [43.8791, -97.5287, 8.77583, -44.2057] 

The DiffEQ documentation says soln.u[1,:] should extract a timeseries of the first component u1 = x (the first column of the above), but instead it extracts the state vector [u1, u2, u3, u4] at the first time step.

julia> soln.u[1,:]
1-element Array{Array{Float64,1},1}:
 [0.0, 1.0, 8.77583, 4.79426]

No permutation of colons, transposes, rows, cols produces what I need, e.g.

julia> soln.u[:,1]
7-element Array{Array{Float64,1},1}:
 [0.0, 1.0, 8.77583, 4.79426]           
 [0.00173129, 1.00095, 8.77583, 4.79232]
 [0.0190442, 1.01038, 8.77583, 4.77299] 
 [0.192173, 1.10263, 8.77583, 4.57966]  
 [1.92346, 1.8154, 8.77583, 2.64632]    
 [19.2363, -12.0343, 8.77583, -16.687]  
 [43.8791, -97.5287, 8.77583, -44.2057] 

How can I convert Array{Array{Float64,1},1} into ````Array{Float64,2}``` and then perform normal colon indexing of rows and cols?

running julia-0.6.0 on Linux, downloaded binary tarball from julialang.org


#2

Answer: use hcat and splat

julia> v = hcat(soln.u[:]...)'
7×4 Array{Float64,2}:
  0.0           1.0      8.77583    4.79426
  0.00173129    1.00095  8.77583    4.79232
  0.0190442     1.01038  8.77583    4.77299
  0.192173      1.10263  8.77583    4.57966
  1.92346       1.8154   8.77583    2.64632
 19.2363      -12.0343   8.77583  -16.687  
 43.8791      -97.5287   8.77583  -44.2057 

julia> v[:,1]
7-element Array{Float64,1}:
  0.0       
  0.00173129
  0.0190442 
  0.192173  
  1.92346   
 19.2363    
 43.8791    

#3

I don’t think it says that anywhere in the documentation? Can you point to where it does? All of the parts of the docs that I check correctly say it’s soln[1,:], i.e. don’t put the .u. Let me know the spot where you found this. When that error is corrected it works fine:

julia> soln[1,:]
7-element Array{Float64,1}:
  0.0
  0.00173129
  0.0190442
  0.192173
  1.92346
 19.2363
 43.8791

Again, there’s no .u here. Let me know where the documentation has this and I’ll fix it, but I cannot find any spots which use a .u.

julia> soln[:,1]
4-element Array{Float64,1}:
 0.0
 1.0
 8.77583
 4.79426

No, don’t use hcat and splat, that will give really bad performance for the conversion. Instead, the docs mention:

At anytime, a true Array can be created using convert(Array,sol).

julia> convert(Array,soln)
4×7 Array{Float64,2}:
 0.0      0.00173129  0.0190442  0.192173  1.92346   19.2363    43.8791
 1.0      1.00095     1.01038    1.10263   1.8154   -12.0343   -97.5287
 8.77583  8.77583     8.77583    8.77583   8.77583    8.77583    8.77583
 4.79426  4.79232     4.77299    4.57966   2.64632  -16.687    -44.2057

#4

You’re right --I’m looking all over the docs and can’t find what I thought I read. I must have seen soln.t[i] and extrapolated my expectations for u into what was actually written. Sorry for the noise, and thanks for the correction.


#5

Thanks for the update. I know I had an error like that in the docs a few months back so it scared me that it was more prevalent. Would you mind updating the solution answer so that way when users search in the future it’ll point to the correct answer + correct conversion method? That’s the best way to make sure this won’t happen again :slight_smile: .


#6

Done, and thanks again.