1

For anyone with interests in dynamical system, i'm trying to plot the phase space of a given hamiltonian.

My code is this:

 H0[x_, y_] := (y*0.5 + 0.5*Log[(x^2 + (y - 1.)^2)/(x^2 + (y + 1.)^2)])

The hamiltonian, which i derive and get the vector field

f[{x_, y_}] := {Derivative[0, 1][H0][x, y], -1*Derivative[1, 0][H0][x, y]}

I use the Runge-Kutta 4th order method for solving the ode:

RK4[h_][z_] := Module[{u1, u2, u3, u4},
u1 = h*f[z];
u2 = h*f[z + u1/2];
u3 = h*f[z + u2/2];
u4 = h*f[z + u3];
z + (u1 + 2.*u2 + 2.*u3 + u4)/6]

Then I generate and plot the orbit:

Orb[x_, h_, T_] := NestList[RK4[h], x, Floor[T/h]]
GraphOrb[x_, h_, T_] := ListPlot[Orb[x, h, T] ]

And it works, but it takes ages to compute a single plot. I think the problem is that Mathematica recalculates everything for every nesting that occurs. I read about the Compile function, but i can't manage to make it work. Can someone help me?

P.S. Yes, I could have plotted the level sets of H0, but i'm interested in this method because i want to parametrize and the perturb H0. And to 10 seconds, for plot with the pertubation it goes to 30 seconds or more.

J. M.'s missing motivation
  • 124,525
  • 11
  • 401
  • 574
Paride Azzari
  • 11
  • 1
  • Derivative and ListPlot won't compile – vapor Jul 04 '16 at 15:07
  • 2
    If your system has a Hamiltonian, why not use a symplectic method instead of insisting on RK? NDSolve[] has symplectic methods built into it. – J. M.'s missing motivation Jul 04 '16 at 15:20
  • 1
    Removing column (:) in the definition of f and adding at the end //FullSimplify speeds up evaluation 4 -fold. – yarchik Jul 04 '16 at 16:18
  • Just a remark: using my answer here, I just do this: sol=hamiltonSolve[(y/2+Log[(x^2+(y-1.)^2)/(x^2+(y+1.)^2)]/2),{{x,1},{y,1},{t,10}}];ParametricPlot[Evaluate[Through[sol[t]]], {t, 0, 10}]. It's very fast although I let it choose the method of integration automatically. Do you really need to program your own RK solver? What are the starting values you tried, and what is the step width? – Jens Jul 04 '16 at 18:54

0 Answers0