# Leapfrog method in Particle-in-cell

Recently, I wanted to write a 3D electromagnetic Particle-in-cell code with C++. I know that I should use leapfrog method, for example, when I calculate the position and velocity of particles，I should calculate according to the following order:

//firt move velocity back by 0.5dt
v = v - 0.5*q/m*E*dt;
//main loop
//update velocity and position
v = v + q/m*E*dt;
x = x + v*dt;

But when I want to write 3D PIC, I don't know the calculation order of $E,\,B,\,v,\,x$ (postion) and their formula. Can someone help me please?

• I'd be inclined to ask on the Computational Science SE – John Rennie Feb 24 '15 at 7:05
• I think the same method is valid - just calculate the 3 components of v instead of 1 when you get to the "calculate velocity" step; and the same for position. In other words - use vector quantities. – Floris Feb 24 '15 at 7:19
• Floris is correct, the three dimensions are treated independently. – lemon Feb 24 '15 at 9:36
• @Floris _ I am wondering if the question is more complicated - rather than just wanting position and velocity as vectors there is also the issue that $E$ and $B$ are determined from the positions of the charges and their velocities and, of course, determine the force on these particles. – tom Feb 24 '15 at 11:33
• I don't think this question is unclear. Nor does it fall afoul of the computational questions policy. – Kyle Kanos Feb 24 '15 at 13:57

For any PIC simulation, you are necessarily tying yourself to particles, particles that experience forces. Thus, we have the generic force law: \begin{align} m_i\frac{d\mathbf v_i}{dt}&=\mathbf F_i \tag{1a}\\ \frac{d\mathbf x_i}{dt}&=\mathbf v_i \tag{1b} \end{align} In the case of PIC, you are often considering the electromagnetic (Lorentz) force, so the force on particle $i$ is, $$\mathbf F_i=q_i\left(\mathbf E+\boldsymbol\beta_i\times\mathbf B\right)\tag{2}$$ with $\boldsymbol\beta=\mathbf v/c$. Since we've got particles, the electric and magnetic fields are defined as $$\mathbf E=\sum_{i\neq j}\frac{q_j}{r_{ij}^2}\hat{r}_j,\qquad\mathbf B=\sum_{j\neq i}\boldsymbol\beta_j\times\mathbf E$$

You've already mentioned the leapfrog method, which is one of the more common methods for PIC simulations (at least to my knowledge). Your question is how/when to compute $\mathbf E,\,\mathbf B$. The answer is quite simple: after you've computed $x$; after all, the electric and magnetic fields are functions of space.

Note, though, that the magnetic field term of the Lorentz force (eq (2)) contains a velocity component; this must also be leapfrogged: $$\mathbf v^{n+1/2}_i\sim\frac12\left(\mathbf v_i^{n-1/2}+\mathbf v_i^{n+1/2}\right)\times\mathbf B$$ which that $\mathbf v_i^{n+1/2}$ on the RHS makes it a little bit harder to solve (you'll likely need a linear algebra solver). There is an alternative, called the Boris method, that combines the velocity and electric field into a single term & you end up with a "rotation" without needing to involve the linear algebra solver needed; however, I am not well versed in this so I cannot say much beyond that.

Thus, the generic stepping algorithm would be something like,

while (t < tend)
computeEMFields(x, v, E, B);
newV = updateVelocity(v, E, B);
newX = moveParticles(v, newV, x);
computeDT();
t = t + dt;
x = newX
v = newV
end while
• Velocity (or position) Verlet methods are also useful here (also called Newmark methods in computational mechanics) -- these are related to the leapfrog method. Also, congrats on answering a question almost half-a-year after you left comments on it! – tpg2114 Jun 26 '15 at 1:45
• @tpg2114: Thank Community, it bumped it a bit ago & reminded me about it! – Kyle Kanos Jun 26 '15 at 1:46

The method you are using here is Euler-Cromer. It only really works for sytems with two variables (although I guess it could be modified for systems with any even number of variables - e.g. find $x, y, z$ then $vx, vy, vz$).

There is a more accurate method, which is more general, but more complicated to program. This is the 4th order Runge Kutta (RK) method. In your leapfrog method you change $x$, then $vx$, then $x$ again then $vx$ and so forth. In the RK method you first find the rates of change of $dx/dt, dy/t, dvx/dt,$ etc. and then find some new rates of change again having moved everything a half step forward in time using the old rates of change. In total four different rates of change for each independent variable are calculated and only then do we move a step forward in time and change $x$, $y$ etc. using a linear combination of all the different rates of change - more complicated to program, but more accurate and general for any number of independent variables.

• RK is "more accurate" only in a limited sense (it is not symplectic!), and thus it is likely that in the problem presented in the question RK would indeed lead to unphysical results (whereas velocity Verlet or leap-frog or even Euler-Cromer would probably work better). – alarge Jun 26 '15 at 1:14