# Solve motion from Hamilton's equations

I have a system of four coordinates and four momenta (conjugates of coordinates).

I have a metric tensor $g_{ik}$ and I know the Hamiltonian:

$H(p,q,t)=\frac{1}{2}g^{ik}(q)p_ip_k$

For my current case I have the Schwartzschild metric:

$ds^2=(1-r_s/r)dt^2-\frac{1}{1-r_s/r}dr^2-r^2d\theta^2-r^2sin^2(\theta)d\phi^2$

Meaning, the Hamiltonian is:

$H=\frac{1}{2}[\frac{1}{1-r_s/r}p_t^2-(1-r_s/r)p_r^2-\frac{1}{r^2}p_\theta^2-\frac{1}{r^2sin^2(\theta)}p_\phi^2]$

As Mathematica code with $r_s = 1$,

h =  1/2 (1/(1 - 1/r) pt^2 - (1 - 1/r) pr^2 - 1/r^2 pθ^2 - 1/(r^2 Sin[θ]^2) pϕ)


The equations of motion by Hamilton's equations are:

$dt/d\tau=\frac{\partial H}{\partial p_t};$

$dr/d\tau=\frac{\partial H}{\partial p_r};$

$d\theta/d\tau=\frac{\partial H}{\partial p_\theta};$

$d\phi/d\tau=\frac{\partial H}{\partial p_\phi};$

$dp_t/d\tau=-\frac{\partial H}{\partial t};$

$dp_r/d\tau=-\frac{\partial H}{\partial r};$

$dp_\theta/d\tau=-\frac{\partial H}{\partial \theta};$

$dp_\phi/d\tau=-\frac{\partial H}{\partial \phi};$

And I want to solve these equations of motion. I've tried deriving the expressions for each partial derivatives of $H$ and then substituting them in these equations manually. However, is there any way to do this dynamically?

What I mean is: suppose, we've only got the expression of Hamiltonian:

$H=H(p_t,p_r,p_\theta,p_\phi,t,r,\theta,\phi)$

And I want to dynamically derive all the partial derivatives and construct system of 8 (in this case) first-order differential equations and then solve them with NDSolve[]. So, that next time, when I want to try something else, all I have to change is the Hamiltonian.

P.s.: I've tried lots of things, like using Function[] for Hamiltonian and throwing variables back and forth, but with no effective result.

• Did you look at NDSolve? (Or maybe DSolve?) If it's a separable Hamiltonian system (I didn't check), there are special integration methods available. See reference.wolfram.com/language/tutorial/NDSolveSplitting.html. Dec 27, 2015 at 19:40
• I don't think it's separable (at least, at first sight, it's not... there are $r$ s everywhere). Yes, I've been fighting with NDSolve[] this whole time. My problem is, I can't differentiate $H$ partially so that it's analytic. Mathematica can't seem to understand that the variables of result functions (partial derivatives) are by themselves functions of $\tau$ and I don't know how to make it understand. Dec 27, 2015 at 19:46
• I've posted an answer based on a function that I wrote earlier. To get a more specific answer regarding your exact problem with differentiation, it would be necessary for you to post the actual Mathematica code you've tried.
– Jens
Dec 27, 2015 at 19:59
• After all my attempts failed, I erased them all and substituted the partial derivatives manually. Now I regret that, though... Dec 27, 2015 at 20:05
• My codes, however, were not half as complex as yours is. I'll have to study it carefully. Thanks for sharing, I appreciate it a lot. Dec 27, 2015 at 20:06

Here is a convenience function that I use to solve for the motion under arbitrary classical Hamiltonians:

hamiltonSolve[hamilton_, varList_] := Module[
{inputVariables, initialValues, phaseSpace, t, tMax, localH},
{t, tMax} = Last[varList];
{inputVariables, initialValues} = Transpose[Most[varList]];
phaseSpace = Through[inputVariables[t]];
localH = hamilton /. Thread[inputVariables -> phaseSpace];
inputVariables /. First[
NDSolve[
D[phaseSpace, t] == Flatten[
Reverse[{-1, 1}*Partition[
D[localH, {phaseSpace}], Length[phaseSpace]/2
]]
]],
Thread[(phaseSpace /. t -> 0) == initialValues
]],
inputVariables,
{t, 0, tMax},
MaxSteps -> Infinity]
]
]

sol = hamiltonSolve[
1/2 (1/(1 - 1/r) pt^2 - (1 - 1/r) pr^2 - 1/r^2 pθ^2 - 1/(r^2 Sin[θ]^2) pϕ),
{{t1, 0.1}, {r, 2}, {θ, 0.1}, {ϕ, 0},
{pt, 1}, {pr, .1}, {pθ, 0}, {pϕ, .1}, {t, 10}}];


This provides all solutions as InterpolatingFunction objects which can be plotted. For example, here is the time as a function of proper time (I copied your Hamiltonian and set $r_s=1$):

Plot[sol[[1]][t],{t,0,10}]


The Hamiltonian is provided as the first argument of hamiltonSolve, and the second argument is a list of all the canonical variables with their initial values. As its last argument, the list additionally contains the name of the proper time variable (here I just call it t because it doesn't appear in the Hamiltonian anyway), together with the desired duration of the integration. The list of canonical variables must be even, and it is assumed that the canonical momenta are listed as the second half of the list (excluding the proper time).

So the list varList is generally of the form

{{x, x0}, {y, y0}, ... , {px, px0}, {py, py0}, ..., {t, tMax}}.

where in your case I entered the time t1 as the first coordinate, so that the corresponding momentum pt is the fifth entry in the list above.

• I wanted to use your Code for my problem but very hard to understand how you construct hamiltonSolve[...]. My Hamilton function is: -0.2 II[t] m1[t] (1 - 0.18 n0[t])^2 SS[t] + m2[t] (-0.02 II[t] + 0.2 II[t] (1 - 0.18 n0[t])^2 SS[t]) + 0.97 (-0.65 II[t] - SS[t] + n0[t] (II[t] + SS[t]))^2 and co-state vars are m1[t], m2[t] and my state variables are SS[t], II[t] and control variable is n0[t]. Initial values: SS[0]=.97, II[0]=.02, m1[0]=2, m2[0]=10. Can you run your function using these parameters? Also, I do not understand why you have {t1, .1} in your function. Dec 23, 2020 at 0:44