# Newbie question about HamiltonianProblem

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

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

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.

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 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

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

1 Like

@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.

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).

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.

But apparently there was a typo in the documentation too Ordering of the arguments q and p in HamiltonianProblem · Issue #451 · SciML/DiffEqDocs.jl · GitHub.

1 Like

I am not there either

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

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.

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