Оптимизация итегратора и Compile [Mathematica]

[sithif]

Помогите разобраться как оптимизизовать интегратор (один шаг) для гамильтоновых систем вида $H=H_1(p)+H_2(q)$
Сам процесс состоит из двух частей:
1) Рекурсивно строится приближенное решение четного порядка $2k+2$, используя метод Yoshida, добавил только сложение коммутирующих частей. В результате получается "отображение" для одного шага интегрирования в виде (YOSHADAONESTEPINTEGRATOR):

{{2.56479 px+x, px, 2.56479 py +y, py}, {...}, ...

2) Итегрирование для заданых нач. условий (PUSHONESTEP)

Нужно максимально оптимизировать оба шага, при этом первый шаг не обязательно должен быть быстрым, так как он считается один раз, а второй -- вызывается несколько раз.

По первому шагу: думаю можно уменьшить длину "отображения", если сделать, например, попарную композицию соседних элементов. Полную композицию всех элементов едва ли можно сделать, это долго.

Для второго шага хочется сделать Compile с CompilationTarget -> "C", но у меня не получается нормально передать "отображение" в Compile. Это "отображение" висит в памяти в виде листа {{},{},...,{}}. Кажется, что можно конвертировать каждый элемент этого листа в Function или в CompileFunction. Получится что-то вроде {f1,f2,...,fn}, тогда в основном Compile можно использовать fn@...@f2@f1[...] без вызова MainEvaluate. Но может что проще можно сделать?

Можно было бы написать one-step метод для NDSolve, так как там есть auto-compile, и не париться, но мне еще нужна касательная динамика и как ее довавить в NDSolve мне не понятно.

Также думал, что проше сделать, например, FortranForm на "отображение" и дальше работать с Fortran (или C), но мне нужна параллельность по разным начальным условиям, и надстраивать MPI не хочется. При этом Compile может сделать Listable, что очень заманчиво :)

Моя стратегия оптимизации пока такая:
1) Добавить, где получится, символьную параллельность в символьное построение "отображения"
2) Конвертировать каждый элемент листа "отображения" в Funtion и сделать частичную композицию.
3) Конвертировать результат композиций в Compile
4) Собрать все части в основном Compile
5) profit

Опыта с Compile у меня почти нет, буду рад Вашим советам:)

Пример:

(* INTEGRATOR ORDER = 2 KORDER + 2 *)
KORDER = 10;
SLICES = 50;
ITERATIONS = 100;
NLENGTH = 100;
NSTRENGTH = 0.0005467499999999999;
(* GENERATE ONE-STEP MAP *)
NHAMILTONIANLIST = {1/2 (px^2 + py^2), -NSTRENGTH (x^2 - y^2)/(x^2 + y^2)^2};
NMAP = YOSHADAONESTEPINTEGRATOR[NHAMILTONIANLIST, NLENGTH/SLICES, KORDER];  // AbsoluteTiming
(* INTEGRATE ONE-STEP *)
INI = {1.621810314392223, 0.009258483424276098, 0.6939644172202891, -0.02163568943615997};
PUSHONESTEP[NMAP, INI] // AbsoluteTiming

Source:

(Оффтоп)

(* GLOBALS *)
(* PHASE-SPACE COORDINATES *)
GCOORDINATES={x,px,y,py};
(* PHASE-SPACE DIMENSION *)
GDIMENSION=Length@GCOORDINATES;
(* SYMPLECTIC 4D-IDENTITY MATRIX *)
GSIMPLECTICIDENTITY={{0,1,0,0},{-1,0,0,0},{0,0,0,1},{0,0,-1,0}};

(* YOSHADA ONE-STEP INTEGRATOR *)
(* YOSHIDA COEFFICIENTS *)
YOSHIDAX[0]=1/2;
YOSHIDAX[k_Integer]:= 1/(2-2^(1/(2k+1))) /; k>0;
YOSHIDAY[0]=1;
YOSHIDAY[k_Integer]:= -2^(1/(2k+1))/(2-2^(1/(2k+1)))/; k>0;
(* YOSHADA COEFFICIENTS PRECISION *)
YOSHADAPRECISION=$MachinePrecision;
(* ONE-STEP MAP GENERATOR *)
YOSHADAONESTEPINTEGRATOR[PARTS_List,STEP_,ORDER_Integer]:=Module[
{
GRADIENTFIELD,
SOLUTIONS,
PARAMETER,
STEPMAP,
OUTPUT
},
GRADIENTFIELD=Table[Simplify@(D[PARTS[[i]],GCOORDINATES[[j]]]),{i,1,Length@PARTS},{j,1,GDIMENSION}];
SOLUTIONS=Table[GCOORDINATES+PARAMETER GSIMPLECTICIDENTITY.GRADIENTFIELD[[i]],{i,1,Length@GRADIENTFIELD}];
STEPMAP[0]= N@Join[
   SOLUTIONS[[1;;-2]]/. PARAMETER -> PARAMETER YOSHIDAX[0],
   {SOLUTIONS[[-1]] /. PARAMETER-> PARAMETER YOSHIDAY[0]},
   Reverse@SOLUTIONS[[1;;-2]] /. PARAMETER-> PARAMETER YOSHIDAX[0]
];
If[ORDER>0,
   {
   STEPMAP[1]=N@Join[
      STEPMAP[0][[1;;-2]] /. PARAMETER-> PARAMETER YOSHIDAX[1],
      {STEPMAP[0][[-1]] /. PARAMETER-> PARAMETER (YOSHIDAX[1]+YOSHIDAY[1])},
      STEPMAP[0][[2;;-2]] /. PARAMETER-> PARAMETER YOSHIDAY[1],
      {STEPMAP[0][[-1]] /. PARAMETER-> PARAMETER (YOSHIDAX[1]+YOSHIDAY[1])},
      Reverse@STEPMAP[0][[1;;-2]] /. PARAMETER -> PARAMETER YOSHIDAX[1]
   ];
   Do[
   {
      STEPMAP[k]=N@Join[
         STEPMAP[k-1][[1;;-2]] /. PARAMETER -> PARAMETER YOSHIDAX[k],
         {STEPMAP[k-1][[-1]] /. PARAMETER -> PARAMETER (YOSHIDAX[k]+YOSHIDAY[k])},
         STEPMAP[k-1][[2;;-2]] /. PARAMETER -> PARAMETER YOSHIDAY[k],
         {STEPMAP[k-1][[-1]] /. PARAMETER -> PARAMETER (YOSHIDAX[k]+YOSHIDAY[k])},
         Reverse@STEPMAP[k-1][[1;;-2]] /. PARAMETER -> PARAMETER YOSHIDAX[k]
      ];
   },{k,2,ORDER}];
   OUTPUT=STEPMAP[ORDER];
   },
   {
   OUTPUT=STEPMAP[ORDER];
   }
];
OUTPUT /. PARAMETER-> STEP
];

(* ONE STEP MAP PUSHER *)
PUSHONESTEP[MAP_List,RAY_List]:= Module[
   {OUTPUT},
   OUTPUT=RAY;
   Do[{OUTPUT=MAP[[i]] /. Thread[GCOORDINATES -> OUTPUT];}, {i,1,Length@MAP}];
   OUTPUT
];