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 http://www.sciencedirect.com/science/article/pii/S0377042701004927)The code is below\[CurlyKappa] = 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}]*\[CurlyKappa]);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 -> \[Infinity]]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
Symplectic Partitioned Runge Kutta method in Mathematica
runge kutta;mathematica