Newbie question about HamiltonianProblem

New to Julia
1

I’m trying to teach myself Julia by using it to solve some differential equations i need solved. I’ve been looking at this example

I’d like to plot the Hamiltonian (or any other function of the variables) as a function of time. The Hamiltonian function is written to take data at a given instant in time, so I could evaluate the function H for each instant in time in a for loop, but I’m very used to the matlab style of vectorized coding, which I also think would be stylistically better. How would I accomplish this?

1 Like
2

Have you tried having a look a the simulation output provided by the double_pendulum.jl code that you are referencing here? In particular, from the first simulation you get

julia> sol1
retcode: Success
Interpolation: automatic order switching interpolation
t: 134-element Array{Float64,1}:
  0.0
  0.005
  0.05499999999999999
  ⋮
 24.71989277613792
 24.86868170266417
 25.0
u: 134-element Array{Array{Float64,1},1}:
 [1.0, 1.0, 1.0, 1.0]
 [0.9173873470403654, 0.9174516957353678, 1.0023967042847444, 1.0000000563149365]
 [0.08644609908203948, 0.09196004700351201, 1.0149015079311146, 1.0000567812684389]
 ⋮
 [-6.054231992074114, -12.056668232377719, 0.027646098435418, 0.10612941900049216]
 [-6.1495227859460195, -11.69036787508306, -6.080301318376474e-5, -0.32640770884302933]
 [-6.256313663905104, -10.239601242478551, -0.1324936832683153, -0.629101072316051]

Obviously, the simulation output comes as some more complicated data type than just a an array (in fact, a matrix) that you are used to from Matlab. But getting the array of values is straigtforward using the dot. Note that the vector of (state) variables is denoted u here:

julia> sol1.u
134-element Array{Array{Float64,1},1}:
 [1.0, 1.0, 1.0, 1.0]
 [0.9173873470403654, 0.9174516957353678, 1.0023967042847444, 1.0000000563149365]
 [0.08644609908203948, 0.09196004700351201, 1.0149015079311146, 1.0000567812684389]
 [-1.9803904064051525, -1.9817102851965904, 0.9544687664602022, 1.0006037605753557]
 [-4.5565663813837265, -4.95682194749687, 0.672002327426011, 0.9832650510821658]
 [-5.9032809983798655, -7.790178097454042, 0.3242060850308398, 0.8759059209509462]
 [-5.943613484795805, -10.775076914267512, 0.022043544720056207, 0.5667541455845946]
 [-5.7130796782664195, -12.013371802990793, -0.0426856254161595, 0.14784486944733183]
 [-5.648958300778203, -11.889774751071544, -0.007217148908584846, -0.25005242797974475]
 ⋮
 [5.587060546627331, 7.806700961000203, 0.37552602998220136, 0.8960073447463703]
 [3.8453217414924703, 4.448148406573896, 0.7144888073156342, 1.0376447446638644]
 [1.286466475379462, 1.376402400305727, 0.9166785732656021, 1.0694682774712758]
 [-2.491704854444083, -2.7950277723838823, 0.8532718060615504, 1.0576089418582855]
 [-5.274180306238754, -6.898084330762543, 0.47881885731568813, 0.9453633766685099]
 [-6.015148096215678, -10.687330502172246, 0.11381677838517706, 0.5827117250218716]
 [-6.054231992074114, -12.056668232377719, 0.027646098435418, 0.10612941900049216]
 [-6.1495227859460195, -11.69036787508306, -6.080301318376474e-5, -0.32640770884302933]
 [-6.256313663905104, -10.239601242478551, -0.1324936832683153, -0.629101072316051]

Note that it is actually an array of arrays, the outer dimension corresponding to time, therefore, the initial state vector is obtained with

julia> sol1.u[1]
4-element Array{Float64,1}:
 1.0
 1.0
 1.0
 1.0

Now, obtaining a function of this sequence of state vectors applied time-wise is easy. Perhaps you may want to get the following function

julia> f(u) = cos(u[1]) + u[2]*u[3] + 5*u[4]
f (generic function with 1 method)

Applying this to the simulation outcomes is straightforward using broadcasting:

julia> f.(sol1.u)
134-element Array{Float64,1}:
  6.54030230586814
  6.527547555276698
  6.089880158973252
  2.7133012143705875
  1.4301365857910806
  2.7826065872357417
  3.5391472089129405
  2.0938666057229924
 -0.358921886704315
  ⋮
  8.178989078710751
  7.603941442605557
  6.889574248219444
  2.106974617987561
  1.9565870812663815
  2.6614536749009683
  1.1712317404892454
 -0.6402472780756148
 -0.7891838986468183

Is this what you wanted?

3 Likes
3

The structure of sol2.u is a little bit different and I don’t understand how to access the parts. For example sol2.u[1] = ([1.0, 1.0], [1.0, 1.0]). These are the momentum and position vectors q and p needed for the Hamiltonian function, but I don’t understand how to understand the components to q and p.

4

I see. Indeed, it seems nontrivial in this case because the solution array is no longer just a simple array of arrays as in the previous case:

julia> sol2.u
501-element Array{ArrayPartition{Float64,Tuple{Array{Float64,1},Array{Float64,1}}},1}:
 ([1.0, 1.0], [1.0, 1.0])
 ([0.16925382765132868, 0.1734695426199092], [1.0155798207154305, 1.002006582279395])
 ([-0.6637297688382037, -0.6532950958415893], [1.0093732723730813, 1.002533724538233])
 ([-1.4875152970292118, -1.4795030569000462], [0.9815477972529919, 1.0014977825985978])
 ([-2.29002188286572, -2.3042184453623173], [0.9328492236209759, 0.9985418344908934])
 ⋮
 ([2.490470760013373, 8.041375539705852], [-0.16416506522085506, -1.402865198861694])
 ([2.535699705758003, 9.093276033199764], [-0.13257354719297346, -1.2883903428980794])
 ([2.5984231952697496, 10.067554278809938], [-0.11335379793647005, -1.1581032116507806])
 ([2.73213939006668, 10.903382717982783], [-0.11058657477542508, -1.0114343635160448])
 ([2.9914689050490897, 11.537448258303435], [-0.1277674836486944, -0.8477504986026044])

In particular, at a given time – say, at the initial time – we have

julia> sol2.u[1]
([1.0, 1.0], [1.0, 1.0])

and the type of this result is

julia> typeof(sol2.u[1])
ArrayPartition{Float64,Tuple{Array{Float64,1},Array{Float64,1}}}

Apparently, this is a data type from the RecursiveArrayTools.jl package, namely the ArrayPartition type. I have just learn it and I am now exploring this along with you :slight_smile: According to the documentation it allows us accessing the “subvectors” by indexing the x attribute (or whatever the accepted terminology is here). Namely,

julia> sol2.u[1].x[1]
2-element Array{Float64,1}:
 1.0
 1.0

julia> sol2.u[1].x[2]
2-element Array{Float64,1}:
 1.0
 1.0

I guess that the former is the q vector and the latter is the p vector in Hamiltonian modeling.

You can also access the desired vector components of the solution sequence conveniently as

julia> sol2[1:2,:]
2×501 Array{Float64,2}:
 1.0  0.169254  -0.66373   -1.48752  -2.29002  -3.05784  …  2.49047  2.5357    2.59842   2.73214   2.99147
 1.0  0.17347   -0.653295  -1.4795   -2.30422  -3.12625     8.04138  9.09328  10.0676   10.9034   11.5374

But now this is essentially a matrix (rows correspond to the individual variables, columns correspond to time).

As a matter of fact, I should have mentioned this fact for the previous case (sol1) too. And this is also exploited in the referenced double_pendulum.jl code within the make_pretty_gif function where they needed this exactly for the purpose of computing the x and y positions from the angles. Perhaps you can just mimic their code.

2 Likes
5

Thanks. That’s very helpful. I think I’m on my way to getting it.

1 Like
6

@zdenek_hurak: thanks, this is very helpful.

What’s the performance impact of ODE solvers returning an m x n Array{Array{Float64,1},1} rather than an m x n Array{Float64,2}? Doesn’t the former have two-level pointer chasing for access and need m system calls for m allocations of m Arrays of length n?

I tried to do some benchmarking but I can’t figure out how to force the full allocation of an m x n Array{Array{Float64,1},1} in one call.

7

I’m still not there yet. I have a few questions and one observation:

  1. After I run double_pendulum.jl, then I can run
p=sol2[1:2,:];
q=sol2[3:4,:];
H0=H(p,q,params)

and this returns the initial value of the Hamiltonian. However I can’t broadcast it because doing so I get the following error message:

Hvec=H.(p,q,params)
ERROR: DimensionMismatch("arrays could not be broadcast to a common size; got a dimension with lengths 2 and 5")
Stacktrace:
 [1] _bcs1 at ./broadcast.jl:501 [inlined]
 [2] _bcs at ./broadcast.jl:495 [inlined]
 [3] broadcast_shape at ./broadcast.jl:489 [inlined]
 [4] combine_axes at ./broadcast.jl:484 [inlined]
 [5] combine_axes at ./broadcast.jl:483 [inlined]
 [6] instantiate at ./broadcast.jl:266 [inlined]
 [7] materialize(::Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{2},Nothing,typeof(H),Tuple{Array{Float64,2},Array{Float64,2},Array{Float64,1}}}) at ./broadcast.jl:837
 [8] top-level scope at REPL[36]:1

So how do I evaluate the Hamiltonian H on the array valued p and q and the params array simultaneously.

  1. In this example the solution is given at a uniformly-spaced array of time steps that were chosen sufficiently close together for the graphs to look resolved, but in many of the examples, a resolved plot is created by doing interpolation based on the ODESolution. How do I evaluate the Hamiltonian function directly on the q and p components of the ODESolution. I see some examples of similar things in this example in this line

energy = map(x->E(x...), sol.u)

but I’m not sure how to do an equivalent thing for a HamiltonianProblem, nor do I know how the ... works in this map.

  1. In the double pendulum example, the Hamiltonian function is entered as H(p, q, params) with p before q but the ODE is defined in this line

prob = HamiltonianProblem(H, q0, p0, times, params)

has q0 before p0. I’m never going to remember when to use which. I my papers, I always write H(q,p).

8

The former was just the author’s arbitrary choice (a bit confusing, I agree with you). The latter (that is, H(q,p)) is more in agreement with the general conventions.

9

But apparently there was a typo in the documentation too https://github.com/SciML/DiffEqDocs.jl/issues/451.

1 Like
10

I am not there either :slight_smile:

Indeed, with sol2[1:2,:] you obtain a matrix and its interpretation as a sequence (or array) of vectors is lost. I am not sure if anything but a for loop (or list comprehension) can be used to evaluate some scalar function of the individual vector elements of the solution sequence (or array).

And I guess that the reason why the solution of HamiltonianProblem is returned in this vector-of-tuples-of-vectors format was exactly to facilitate solution of your problem. But I do not know how :slight_smile:

A solution that I can come up with is to introduce dedicated variables for each scalar variable:

julia> p1 = sol2[1,:]
501-element Array{Float64,1}:
  1.0
  0.16925382765132854
 -0.6637297688382036
 -1.4875152970292118
 -2.2900218828657195
 -3.0578405129575437
  ⋮
  2.4219339977447945
  2.4904707600133404
  2.535699705757957
  2.5984231952696906
  2.732139390066608
  2.991468905049008

julia> p2 = sol2[2,:]
501-element Array{Float64,1}:
  1.0
  0.1734695426199092
 -0.6532950958415893
 -1.4795030569000462
 -2.304218445362318
 -3.1262469177392993
  ⋮
  6.959103420372306
  8.041375539705818
  9.09327603319974
 10.067554278809917
 10.903382717982769
 11.537448258303423

julia> q1 = sol2[3,:]
501-element Array{Float64,1}:
  1.0
  1.0155798207154303
  1.009373272373081
  0.9815477972529917
  0.9328492236209757
  0.8647494469914832
  ⋮
 -0.2039548396853722
 -0.16416506522084431
 -0.13257354719296213
 -0.11335379793645835
 -0.11058657477541316
 -0.12776748364868243

julia> q2 = sol2[4,:]
501-element Array{Float64,1}:
  1.0
  1.002006582279395
  1.002533724538233
  1.0014977825985978
  0.9985418344908934
  0.9929106161122214
  ⋮
 -1.5019899786777668
 -1.402865198861699
 -1.288390342898084
 -1.158103211650785
 -1.0114343635160494
 -0.8477504986026092

and then you can apply broadcasting with any scalar function of the individual components. Say

julia> f(p1,p2,q1,q2) = p1.^2 + p2.*q1 + cos.(q2)
f (generic function with 1 method)

julia> f(p1,p2,q1,q2)
501-element Array{Float64,1}:
 2.5403023058681398
 0.7434317636738219
 0.3192871154893914
 1.2995401520634888
 3.636240571186884
 7.193222498581859
 ⋮
 4.515173536139257
 5.049474606198497
 5.502912268819041
 6.011685379726895
 6.78946340847042
 8.136446959124365

PS: I believe your question could have received more attention with a more descriptive title. Perhaps something like “Format of simulation outputs from HamiltonianProblem using DifferentialEquations” or something like that.

11

Thanks @zdenek_hurak.

Your solution will work well if the number of degrees of freedom in your code is fixed. I’m working on a project with N interacting vortices, which is not quite appropriate for NBodySimulation.jl since it’s not a Newtonian system, and I’d like to be able to write codes that don’t require modfication as I vary N.

Regarding your PS. I think I’ll write a second post with a more descriptive title, referencing this thread and succinctly re-asking the questions you and I still have. It sounds as if the issues I have brought up are non-trivial and worthy of wider attention.

1 Like