1
$\begingroup$

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

$\endgroup$
  • $\begingroup$ Welcome to SciComp.SE. Can you write the equations for the system that you are trying to solve? Have you tried to make a substitution of the terms that Mathematica can't differentiate? $\endgroup$ – nicoguaro Aug 9 '15 at 6:11
  • $\begingroup$ @nicoguaro 1. Done. 2. I have - see the description. $\endgroup$ – gidman Aug 9 '15 at 22:17
  • $\begingroup$ Would this be better/more likely to be answered on mathematica.stackexchange.com? $\endgroup$ – dr.blochwave Aug 10 '15 at 7:00