I want to plot the Poincare map for Arnold-Beltrami-Childress magnetic field for parameters $A=1, B=0.816, C=0.5773$ in Python for the Poincare section $z=0$.

Also, I am not able to understand what initial conditions should I vary for $x$, $y$, and $z$ so as to plot the entire map? I am able to get only a subset of the plot for $x=0$, $y=0$, and $z=0$ but how should I vary $x$, $y$, and $z$?

ODEs that I have used for the same:


where $\mathbf{r}=(x,y,z)$ and $ \mathbf{B}=\left( \begin{array}{c} A\sin(z)+C\cos(y)\\ B\sin(x)+A\cos(z)\\ C\sin(y)+B\cos(x) \end{array} \right) $.

Can you please provide me with the Python code?

  • 1
    $\begingroup$ Would Computational Science be a better home for this question? $\endgroup$
    – Qmechanic
    Commented Jul 7, 2018 at 16:39
  • $\begingroup$ Yes, it belongs there $\endgroup$
    – Nakul Aggarwal
    Commented Jul 7, 2018 at 22:20

1 Answer 1


One first observation is that as $x,y,z$ are angles, one should take their values mod $2\pi$. Without that, the trajectory of the initial point $[0,0,0]$ exhibits this interesting pattern

def ABC_ode(u,t):
    A, B, C = 1, 0.816, 0.5773
    x, y, z = u
    B = np.array([A*np.sin(z)+C*np.cos(y), B*np.sin(x)+A*np.cos(z), C*np.sin(y)+B*np.cos(x)])
    return B/(1e-12+sum(B**2))**0.5

def mysolver(u0, tspan): return odeint(ABC_ode, u0, tspan, atol=1e-10, rtol=1e-11)

trajectory of the solution starting in origin

So one has to reduce the coordinates by removing multiples of $2\pi$, but then also exclude the so-introduced "long jumps" by about $2\pi$ from the root computation.

To fill the cross section at $z=0$ maximally in a reasonable time, that is without computing similar trajectories, divide the square into smaller squares and remember for each small square in the grid if it was already visited. Only start trajectories from squares without visit. Here linear interpolation is used to find the intersection points, one could also use cubic interpolation with the derivative values or an additional integration step with smaller step size to get a more accurate intersection location. But at the level of the plot that is not necessary.

u0 = [0,0,0]
t = np.arange(0, 1500, 0.2)

grid = np.zeros([N,N], dtype=int)
for i in range(N):
    for j in range(N):
        if grid[i,j]>0: continue;
        x0, y0 = (2*i+1)*pi/N-pi, (2*j+1)*pi/N-pi 
        u0 = [x0,y0,0]
        px = []
        py = []

        u = mysolver(u0, t); u0 = u[-1]
        u = np.mod(u+pi,2*pi)-pi
        x,y,z = u.T

        for k in range(len(z)-1): 
            if z[k]<=0 and z[k+1]>=0 and z[k+1]-z[k]<pi:
                s = -z[k]/(z[k+1]-z[k])  # 0 = z[k] + s*(z[k+1]-z[k])
                rx, ry = (1-s)*x[k]+s*x[k+1], (1-s)*y[k]+s*y[k+1]
                m, n = int((rx+pi)*N/(2*pi)), int((ry+pi)*N/(2*pi))

        #plt.plot(px, py, '-k', lw=0.2)
        plt.plot(px, py, '.', ms=3)

enter image description here


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.