0
$\begingroup$

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:

$$\frac{d\mathbf{r}}{dt}=\frac{\bf{B}}{|\mathbf{B}|}$$

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?

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

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)

N=54
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]
                px.append(rx); 
                py.append(ry);
                m, n = int((rx+pi)*N/(2*pi)), int((ry+pi)*N/(2*pi))
                grid[m,n]=1

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

enter image description here

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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