计算科学 - Mathematica 中的辛分区 Runge Kutta 方法 - 吾爱随笔录

我试图解决哈密顿系统（是所有广义坐标的向量， - 广义动量）其中 $Q$ $P$

\frac{d Q}{d t} = \frac{\partial H}{\partial P} \frac{d P}{d t} = - \frac{\partial H}{\partial Q}

$\frac{\mathrm{d} Q}{\mathrm{d} t}=\frac{\partial H}{\partial P} \\ \frac{\mathrm{d} P}{\mathrm{d} t}=-\frac{\partial H}{\partial Q}$

H = \frac{1}{2} \sum_{i = 1}^{n} {\vec{p}}_{i}^{2} + \frac{k}{2} \sum_{i = 2}^{n} (‖ {\vec{q}}_{i} - {\vec{q}}_{i - 1} ‖ - r e s t)^{2} + \frac{κ}{2} \sum_{i = 2}^{n - 1} \arccos^{2} \frac{({\vec{q}}_{i} - {\vec{q}}_{i - 1}) \cdot ({\vec{q}}_{i + 1} - {\vec{q}}_{i})}{‖ {\vec{q}}_{i} - {\vec{q}}_{i - 1} ‖ ‖ {\vec{q}}_{i + 1} - {\vec{q}}_{i} ‖}

$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 \|}$ 使用 Mathematica 以数字方式使用“SymplecticPartitionedRungeKutta”方法（可在此处找到说明http://www.sciencedirect.com/science/article/pii/S0377042701004927）

代码如下

\[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]]

但是，错误弹出：

NDSolve::srpksep：“该方法中微分系统的哈密顿量似乎不是可分离的形式。尝试使用带有系数 ImplicitRungeKuttaGaussCoefficients 的方法 ImplicitRungeKutta。”

但是很容易看出，系统的哈密顿量肯定有单独的形式 H(p,q)=T(p)+V(q)。

任何人都可以帮忙吗？

PS为了防止提问。你会注意到我定义了两个“奇怪”的函数：

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

想法：Mathematica 不从非解析函数中求导。当你对形式求导时D[H,q[t]]，就会出现类型术语D[Conjugate[q[t]],q[t]]。如果我不修复它，我会得到一个错误：

NDSolve::ndnum: 在 t == 0 处遇到导数的非数值