One of the things that has frustrated me about Physibounce is that when the magnetic field is on, the radius of the cyclotron motion of the particle isn't constant. Here I will write about why this problem occurred and what can be done to fix it.
In Physibounce, the magnetic fields are constant in time (other than turning on or off). When the magnetic field is turned on, the magnetic field can either point into or out of the screen. Let me call this the direction. The equations of motions for the particles in Physibounce are
where with is the magnitude of the magnetic field, is the charge of the particle, and the mass. is all forces acting on the particle that are not the magnetic field. Before now, Physibounce used an semi-implicit Euler method. Let us see why the semi-implicit Euler method is bad. For simplicity, take the case when there are no forces other than the magnetic field applied. The update equation for semi-implicit Euler is
Assuming I measure distances from the center of the circular motion of the particle, ideally the radius of the orbit does not change. Since the magnetic force is causing the circular motion, it must be that or equivalently, Thus, enforcing that the radius of orbit does not change is equivalent to enforcing that kinetic energy is conserved as (remember that magnetic fields do no work). Consider the velocity squared at the updated timestep
Thus, every timestep, a mistake of is made to the kinetic energy. Since this is proportional to the current kinetic energy, this means that the error gets worse and worse as time goes on, which is why semi-implicit Euler is bad for modeling magnetic forces. This is also why in Physibounce the radius of the orbit increases for a particle in a magnetic field.
One solution to the problem would be to use polar coordinates. Then the semi-implicit Euler update would be
which looks quite simple. In this form, it is obvious that the radius of rotation is conserved. However, if I wanted to use this method, I would have to write all of the other forces of Physibounce in polar coordinates. This would be quite a bit of work to work out, so I did not choose this method.
The solution to this problem while still using Cartesian coordinates is to use a leap-frog-like scheme, specifically one called the Boris pusher. The idea of the Boris scheme is to push the velocity by a half timestep using only the force then make a rotation due to the magnetic field, and then push the velocity another half timestep using the same force as before . As before, the position is updated with the updated velocity. I calculate the following quantities
Then I update
I will again consider the case when the magnetic field is the only force, but the more general case is considered in . In this case, the update has the nice property that
Thus, kinetic energy is conserved, and hence the radius of orbit of the particle is constant. This creates a more realistic simulation and does not make the computation of all the other forces in Physibounce any more complicated. There is a minimal amount of code required to incorporate it into a simulation already using a semi-implicit Euler method. These pluses make it ideal for solving the problem in Physibounce.
The Boris pusher is widely used in plasma physics simulations as it is quite fast (in terms of number of flops required) and represents the dynamics without too much error. It is used in N-body and PIC simulations.
I have updated Physibounce (first update in about a year!) to reflect this.
1. Qin, H. et al., 2013. Why is Boris algorithm so good? Phys. of Plas. 20-8.