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$

closed as off-topic by Paul Aug 11 '15 at 1:49

  • This question does not appear to be about computational science within the scope defined in the help center.
If this question can be reworded to fit the rules in the help center, please edit the question.

  • $\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$ – Dmitry 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