# Numerically Integrating a Hamiltonian but getting different results when compared with alternative equivalent equations

currently trying to perform an integration of the Hamiltonian for a Schwarzschild black hole. However, I'm coming up short.

The code is as follows:

   ClearAll["Global*"]
Needs["DifferentialEquationsNDSolveProblems"];
Needs["DifferentialEquationsNDSolveUtilities"];

(*Symplectic schemes*)
q = {t[s], r[s], \[Theta][s], \[Phi][s]};
p = D[q, s];
n = Length[q];
tt = 1 - 2 m/r[s];
rr = -1/tt;
\[Theta]\[Theta] = -r[s]^2;
\[Phi]\[Phi] = -(r[s] Sin[\[Theta][s]])^2;
metric = {{tt, 0, 0, 0}, {0, rr, 0, 0}, {0, 0, \[Theta]\[Theta],
0}, {0, 0, 0, \[Phi]\[Phi]}};
inversemetric = Simplify[Inverse[metric]];
ivs = {1.3, 0, 0, 0.088}; ics = {0, 6.5, \[Pi]/2, 0};
m = 1;

hamiltonian = 1/2 (D[q, s].inversemetric.D[q, s]) == 1;
pdot = Table[-D[hamiltonian, q[[i]]], {i, 1, n}];
qdot = Table[ D[hamiltonian, p[[i]]], {i, 1, n}];

eqs1 = {{qdot, pdot}, {ics, ivs}};

invariants = hamiltonian;
time = {T, 0, 100};
solee = NDSolve[eqs1, q, time, Method -> "ExplicitRungeKutta",
StartingStepSize -> 0.1]


However, when I attempt to numerically integrate I get the following message:

NDSolve::deqn: Equation or list of equations expected instead of -True in the first argument {{{(t^\[Prime])[s]/(1+Times[<<2>>])==0,(-1+2 Power[<<2>>]) (r^\[Prime])[s]==0,-((\[Theta]^\[Prime])[s]/r[s]^2)==0,-((Csc[\[Theta][<<1>>]]^2 (\[Phi]^\[Prime])[s])/r[s]^2)==0},{-True,-(1/2 (Times[<<3>>]+Times[<<4>>]+Times[<<3>>]+Times[<<4>>])==0),-((Cot[\[Theta][<<1>>]] Csc[<<1>>]^2 (<<1>>^(<<1>>))[<<1>>]^2)/r[<<1>>]^2==0),-True}},{{0,6.5,\[Pi]/2,0},{1.3,0,0,0.088}}}.


UPDATE (initial problem solved) Following the excellent answer below, we can see that the singularity which occurs at r=2 which is reached around s=10...

Now, I have written a version of this code which performs integration of the geodesic equations of motion (as opposed to the Hamiltonian equations above) under the same initial conditions. However, there is no singularity reached for this and the results are clearly different. The results should be equivalent! I am not too sure what is the issue.

• In the code you posted, coords is undefined. – Michael E2 Jun 5 '17 at 2:48
• Have you seen this? – Michael E2 Jun 5 '17 at 2:49
• D[hamiltonian, q[[1]]] returns true instead of an equation, that's what ndsolve complains about in the message – tsuresuregusa Jun 5 '17 at 3:24
• and have you seen the output of rr or tt? it's an array with half a million elements. I think you want to define a function istead, no? – tsuresuregusa Jun 5 '17 at 3:29
• Dimensions[tt] I tried writing them as functions but it killed my kernel, let me try again in a bit. – tsuresuregusa Jun 5 '17 at 3:36

Here's my guess at how it should be done. There were a couple of changes. You need separate q and p variables in order to write the Hamiltonian and set up the corresponding equations.

(*Symplectic schemes*)
Clear[m];
q = {t[s], r[s], θ[s], ϕ[s]};
(*p=D[q,s];*)
p = {pt[s], pr[s], pθ[s], pϕ[s]};                 (* use explicit p variables *)
n = Length[q];
tt = 1 - 2 m/r[s];
rr = -1/tt;
θθ = -r[s]^2;
ϕϕ = -(r[s] Sin[θ[s]])^2;
metric = {{tt, 0, 0, 0}, {0, rr, 0, 0}, {0, 0, θθ, 0}, {0, 0, 0, ϕϕ}};
inversemetric = Simplify[Inverse[metric]];
ivs = {1.3, 0, 0, 0.088}; ics = {0, 6.5, π/2, 0};
m = 1;

(*hamiltonian=1/2 (D[q,s].inversemetric.D[q,s])-1;*)

hamiltonian = 1/2 (p.inversemetric.p) - 1;        (* define H in terms of p, q *)
pdot = Table[-D[hamiltonian, q[[i]]], {i, 1, n}];
qdot = Table[D[hamiltonian, p[[i]]], {i, 1, n}];

eqs1 = {{D[q, s] == qdot, D[p, s] == pdot},       (* adjust equations to new defs *)
{(q /. s -> 0) == ics, (p /. s -> 0) == ivs}};

invariants = hamiltonian;
time = {s, 0, 100};                               (* change T to s to match q etc. *)
solee = NDSolve[eqs1, q, time, Method -> "ExplicitRungeKutta",
StartingStepSize -> 0.1]


NDSolve::ndstf: At s == 10.620493252880072, system appears to be stiff. Methods Automatic, BDF, or StiffnessSwitching may be more appropriate.

ListLinePlot[Head /@ (q /. First@solee), PlotRange -> {-0.5, 30}, Frame -> True]


You get a slightly better result with a Gauss-Legendre collocation method of moderate order:

solee = NDSolve[eqs1, Head /@ q, time,
Method -> {"ImplicitRungeKutta", "DifferenceOrder" -> 6,
"Coefficients" -> "ImplicitRungeKuttaGaussCoefficients"}]


NDSolve::ndsz: At s == 10.834214056862363, step size is effectively zero; singularity or stiff system suspected.

The "SymplecticPartitionedRungeKutta"` method could not be used because the system did not appear to be separable. I did not look further into it.

• Wow. Very nice. Although I am confused. I have an alternate code where I integrate the geodesic equation with the same initial conditions using an explicit runge kutta, however, I get different results? They should be equivalent. For one, there is no singularity observed at all.... Hmmmm, I will need to keep checking to see what is going on. – Rumplestillskin Jun 5 '17 at 5:05
• So, for one simulation which uses a code I wrote adopted from here the orbit never approaches the singularity located at r=2 and looks okay. For the above method the singularity is reached? I'm not entirely sure why this is happening. – Rumplestillskin Jun 5 '17 at 6:33
• @Rumplestillskin I got equivalent results with both the explicit and implicit RK methods (as shown in my answer). My plot looks like the plot at the top of page 9 of the Tollerud pdf, singularity and all. Maybe it's like that or maybe there's a difference in the matrix or the ICs. My first guess is that the difference is not due to the numerical method but is in the math. – Michael E2 Jun 5 '17 at 14:24
• from what I believe the output should be equivalent to the output at the bottom of pg. 7 and the top of page 8 of the Tollerud pdf. What we have is a particle with angular velocity 0.088 starting at r=6.5 so it should make some sort of orbit shape! But I am not sure why it is not doing that. Although your answer is fantastic and much appreciated. – Rumplestillskin Jun 6 '17 at 1:08