SymplecticPartitionedRungeKutta shows strange error

Question

Bug introduced in 9.0 and fixed in 11.3.0 or earlier

---

I tried to solve Hamiltonian system ($Q$ is a vector of all generalized coordinates, $P$ - of generalized momentum) $$ \frac{\mathrm{d} Q}{\mathrm{d} t}=\frac{\partial H}{\partial P} \\ \frac{\mathrm{d} P}{\mathrm{d} t}=-\frac{\partial H}{\partial Q} $$ where $$ H=\frac{1}{2} \sum_{i=1}^{n}\overrightarrow{p}_{i}^{2} + \frac{k}{2} \sum_{i=2}^{n}(\left \| \overrightarrow{q}_{i}-\overrightarrow{q}_{i-1} \right \|-rest)^2 + \frac{\kappa}{2} \sum_{i=2}^{n-1}\arccos^2\frac{(\overrightarrow{q}_{i}-\overrightarrow{q}_{i-1})\cdot (\overrightarrow{q}_{i+1}-\overrightarrow{q}_{i})}{\left \| \overrightarrow{q}_{i}-\overrightarrow{q}_{i-1} \right \|\left \| \overrightarrow{q}_{i+1}-\overrightarrow{q}_{i} \right \|} $$ with "SymplecticPartitionedRungeKutta" method numerically using Mathematica (description could be found here https://en.wikipedia.org/wiki/Symplectic_integrator http://www.sciencedirect.com/science/article/pii/S0377042701004927)

The code is below

ϰ = 20;
k = 20;
rest = Sqrt[5];
n = 3;
dim = 3;
Q = Table[Table[Subscript[q, j][i][t], {j, 1, dim}], {i, 1, n}];
P = Table[Table[Subscript[p, j][i][t], {j, 1, dim}], {i, 1, n}];
icsQ = {{-1, -1, 0}, {0, 1, 0}, {1, -1, 0}};
icsP = Table[Table[0, {j, 1, dim}], {i, 1, n}];

H = 1/2 (Sum[P[[i]].P[[i]], {i, 1, n}] + 
     Sum[(Norm[Q[[i]] - Q[[i - 1]]] - rest)^2, {i, 2, n}]*k + 
     Sum[(VectorAngle[Q[[i]] - Q[[i - 1]], 
         Q[[i + 1]] - Q[[i]]])^2, {i, 2, n - 1}]*ϰ);
Q' = Flatten[D[#, t] & /@ Q];
P' = Flatten[D[#, t] & /@ P];

Conjugate'[y_[x_][t_]] := 1;
Abs'[x_] := Sign[x];

var = Flatten[Q~Join~P];
ics = Flatten[icsQ~Join~icsP];
eqn = Table[Q'[[i]] == D[H, var[[n*dim + i]]], {i, 1, n*dim}] ~Join~
   Table[P'[[i]] == -D[H, var[[i]]], {i, 1, n*dim}] ~Join~
   Table[(var[[i]] /. t -> 0) == ics[[i]], {i, Length[var]}];

method = {"SymplecticPartitionedRungeKutta", "PositionVariables" -> var[[1 ;; dim*n]]};
sol = NDSolve[eqn, var, {t, 0, 10}, 
  Method -> method, WorkingPrecision -> 10, 
  MaxStepSize -> 0.001, MaxSteps -> ∞]

However, error pops out:

NDSolve::sprksep: "The Hamiltonian of the differential system in the method does not appear to be in separable form. Try using the method ImplicitRungeKutta with coefficients ImplicitRungeKuttaGaussCoefficients."

But it's easy to see that the Hamiltonian of the system definately has separate form H(p,q)=T(p)+V(q).

Can anyone help?

P.S. In order to prevent questions. You could notice that I defined two "strange" functions:

 Conjugate'[y_[x_][t_]] := 1;
 Abs'[x_] := Sign[x];

Idea: Mathematica doesn't take the derivative from non-analytical functions. And when you take derivative in form D[H,q[t]], terms of type D[Conjugate[q[t]],q[t]] appear. If I don't fix it, I will get an error:

NDSolve::ndnum: Encountered non-numerical value for a derivative at t == 0

Answer 1

I think there's a bug in the internal function NDSolve`SPRKDump`CheckSeparability that leads NDSolve to conclude that the system is not separable. I think you should report it and see if WRI can verify it (they would probably appreciate a link to this Q&A). It's a fair amount of work to track it down, and there is a lot of nearly unreadable stuff to wade through on the way.

To check separability, one can evaluate the P' derivatives at only the initial values of Q and, vice versa, the Q' derivatives at only the initial values of P. In both cases, if the result is a vector of numbers, then the system is separable. Unfortunately, NDSolve checks P' by plugging in the initial values of P for the variables Q. In many cases, this won't be a problem, since the number of inputs is the same and the results are thrown away anyway. But in this case, plugging the initial values of P into the variables Q causes divide-by-zero errors. Consequently the results for P' are not numeric and NDSolve says the system is not separable.

Here's a fix. NDSolve evaluates a result, but it will be up to the OP to determine whether it is reasonable result. The fix, as far as I can tell, only affects testing the system during initialization of the method; it should not affect the computation of the solution. (The fix consists of switching q0 and p0 in the two commented lines, basically a two-character fix.)

cs[ndstate_, pnf_, ppos_, qnf_, qpos_] :=     (* checks separability *)
 Module[{dir, y0, p0, pf0, q0, qf0, sd}, dir = 1;
  sd = ndstate["SolutionData"[dir]];
  y0 = ndstate["SolutionVector"[dir]];
  q0 = y0[[qpos]];
  p0 = y0[[ppos]];
  Quiet[
   NDSolve`SetSolutionDataComponent[sd, "X", p0];         (* Evaluate Q' at P0 *)
   qf0 = NDSolve`EvaluateWithSolutionData[qnf, sd];
   NDSolve`SetSolutionDataComponent[sd, "X", q0];         (* Evaluate P' at Q0 *)
   pf0 = NDSolve`EvaluateWithSolutionData[pnf, sd];];
  If[! VectorQ[qf0, NumberQ] || ! VectorQ[pf0, NumberQ], 
   NDSolve`NDSolveMessage[ndstate, "sprksep", "SymplecticPartitionedRungeKutta"];
   Throw[$Failed];];
  ]

Block[{NDSolve`SPRKDump`CheckSeparability = cs},
 NDSolve[eqn, var, {t, 0, 10}, Method -> method, 
  WorkingPrecision -> 10, MaxStepSize -> 0.001, MaxSteps -> ∞]
 ]