# Newbie question about HamiltonianProblem

**URL:** <https://discourse.julialang.org/t/newbie-question-about-hamiltonianproblem/57492>\
**Category:** New to Julia\
**Tags:** diffeq\
**Created:** [March 18, 2021, 6:54pm UTC](https://discourse.julialang.org/t/newbie-question-about-hamiltonianproblem/57492 "2021-03-18T18:54:07Z")\
**Posts on this page:** 11\
**Page:** 1

<div class="post-metadata">

**Author:** ![manroygood](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/manroygood/32/23056_2.png) [@manroygood](https://discourse.julialang.org/u/manroygood)\
**Post date:** [March 18, 2021, 6:54pm UTC](https://discourse.julialang.org/t/newbie-question-about-hamiltonianproblem/57492/1 "2021-03-18T18:54:07Z")

</div>

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](https://github.com/SciML/DiffEqPhysics.jl/blob/master/examples/double_pendulum.jl)

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?

---

<div class="post-metadata">

**Author:** ![zdenek\_hurak](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/zdenek_hurak/32/53118_2.png) [@zdenek\_hurak](https://discourse.julialang.org/u/zdenek_hurak)\
**Post date:** [March 18, 2021, 10:42pm UTC](https://discourse.julialang.org/t/newbie-question-about-hamiltonianproblem/57492/2 "2021-03-18T22:42:53Z")

</div>

Have you tried having a look a the simulation output provided by the [double\_pendulum.jl](https://github.com/SciML/DiffEqPhysics.jl/blob/master/examples/double_pendulum.jl) code that you are referencing here? In particular, from the first simulation you get

```julia
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
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
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
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](https://docs.julialang.org/en/v1/manual/arrays/#Broadcasting):

```julia
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?

---

<div class="post-metadata">

**Author:** ![manroygood](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/manroygood/32/23056_2.png) [@manroygood](https://discourse.julialang.org/u/manroygood)\
**Post date:** [March 19, 2021, 2:46am UTC](https://discourse.julialang.org/t/newbie-question-about-hamiltonianproblem/57492/3 "2021-03-19T02:46:51Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![zdenek\_hurak](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/zdenek_hurak/32/53118_2.png) [@zdenek\_hurak](https://discourse.julialang.org/u/zdenek_hurak)\
**Post date:** [March 19, 2021, 1:18pm UTC](https://discourse.julialang.org/t/newbie-question-about-hamiltonianproblem/57492/4 "2021-03-19T13:18:55Z")

</div>

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

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

```

and the type of this result is

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

```

Apparently, this is a data type from the [RecursiveArrayTools.jl](https://github.com/SciML/RecursiveArrayTools.jl) package, namely the [ArrayPartition](https://github.com/SciML/RecursiveArrayTools.jl#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
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
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](https://github.com/SciML/DiffEqPhysics.jl/blob/master/examples/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.

---

<div class="post-metadata">

**Author:** ![manroygood](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/manroygood/32/23056_2.png) [@manroygood](https://discourse.julialang.org/u/manroygood)\
**Post date:** [March 19, 2021, 1:27pm UTC](https://discourse.julialang.org/t/newbie-question-about-hamiltonianproblem/57492/5 "2021-03-19T13:27:19Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![John\_Gibson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/john_gibson/32/5321_2.png) [@John\_Gibson](https://discourse.julialang.org/u/John_Gibson)\
**Post date:** [March 19, 2021, 2:44pm UTC](https://discourse.julialang.org/t/newbie-question-about-hamiltonianproblem/57492/6 "2021-03-19T14:44:46Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![manroygood](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/manroygood/32/23056_2.png) [@manroygood](https://discourse.julialang.org/u/manroygood)\
**Post date:** [March 19, 2021, 10:31pm UTC](https://discourse.julialang.org/t/newbie-question-about-hamiltonianproblem/57492/7 "2021-03-19T22:31:13Z")

</div>

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

```julia
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:

```julia
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](https://tutorials.sciml.ai/html/models/01-classical_physics.html) 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).

---

<div class="post-metadata">

**Author:** ![zdenek\_hurak](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/zdenek_hurak/32/53118_2.png) [@zdenek\_hurak](https://discourse.julialang.org/u/zdenek_hurak)\
**Post date:** [March 19, 2021, 10:37pm UTC](https://discourse.julialang.org/t/newbie-question-about-hamiltonianproblem/57492/8 "2021-03-19T22:37:33Z")

</div>

> [@manroygood](#):
>
> 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.

---

<div class="post-metadata">

**Author:** ![zdenek\_hurak](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/zdenek_hurak/32/53118_2.png) [@zdenek\_hurak](https://discourse.julialang.org/u/zdenek_hurak)\
**Post date:** [March 20, 2021, 9:54am UTC](https://discourse.julialang.org/t/newbie-question-about-hamiltonianproblem/57492/9 "2021-03-20T09:54:26Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![zdenek\_hurak](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/zdenek_hurak/32/53118_2.png) [@zdenek\_hurak](https://discourse.julialang.org/u/zdenek_hurak)\
**Post date:** [March 22, 2021, 1:09am UTC](https://discourse.julialang.org/t/newbie-question-about-hamiltonianproblem/57492/10 "2021-03-22T01:09:46Z")

</div>

> [@manroygood](#):
>
> 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
> 
> ```julia
> 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:
> 
> ```julia
> 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.

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
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
julia> f(p1,p2,q1,q2) = p1.^2 + p2.*q1 + cos.(q2)
f (generic function with 1 method)

```

```julia
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.

---

<div class="post-metadata">

**Author:** ![manroygood](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/manroygood/32/23056_2.png) [@manroygood](https://discourse.julialang.org/u/manroygood)\
**Post date:** [March 22, 2021, 3:15pm UTC](https://discourse.julialang.org/t/newbie-question-about-hamiltonianproblem/57492/11 "2021-03-22T15:15:05Z")

</div>

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.