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

  • $\begingroup$ I think this is a bug, but I could be mistaken. If so, it goes back to V9.0.1 or earlier. $\endgroup$
    – Michael E2
    Aug 11 '15 at 20:30
  • $\begingroup$ WRI tech support confirmed the bug. $\endgroup$
    – Michael E2
    Aug 18 '15 at 23:58

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]];
   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"];

Block[{NDSolve`SPRKDump`CheckSeparability = cs},
 NDSolve[eqn, var, {t, 0, 10}, Method -> method, 
  WorkingPrecision -> 10, MaxStepSize -> 0.001, MaxSteps -> ∞]
  • $\begingroup$ at least now it's not throwing an error. Thanks for the reply. However, solution is not right. It must be oscillatory, but mine, using your code, is not. It exponentially diverges. $\endgroup$
    – gidman
    Aug 13 '15 at 23:10
  • 2
    $\begingroup$ @DmitryGrinko That's what I got, and it seems wrong to me, too. But the change above is only called in the initialization stage, as far as I can tell. It has nothing to do with the actual integration of the system. So I don't think the divergence is because of my code (for example, cs does not even return a value). It might be because you chose a somewhat low WorkingPrecision combined with a small step size. In any case, FWIW, my first guess is that its a numerics issue. $\endgroup$
    – Michael E2
    Aug 13 '15 at 23:28

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.