Studying ALF, the asynchronous leap frog integrator.

Ulrich Mutze 2023-06-11

The present project collect a few simulations which all are based on the template class implementation of the ALF method as defined in file cpm1/cpmalf.h.

This simulations are defined in file alftest.cpp as functions test1() to test7(). There themes are:

test1(): A quantum particle moving in a discrete one-dimensional biotope hits a smooth potential well and, in a second video, a potential wall.

test2(): Study of the 'flame ball growth equation' in one dimension. This is a well-known stiff ordinary differential equation.

test3(): Dynamics of n classical particles of the same sort. Motion is in the two dimensional continuum from stochastic initial conditions. The particles interact via rotationally invariant pair potentials and a rotationally invariant external field emanating from a point source at the center of the biotope.

test4(): Free fall of a single particle. Simple system with external field (gravity).

test5(): System of particles moving in a 2D disc to which they remain confined thjrough an external potential. The particles start from a regular grid of positions and velocities which all have the same absolute value but random directions. After some span of time, the sign of the time step value is reversed. The reversibility of the ALF integrator has the effect that all particles move back to their initial positions. The accuracy with which this happens is limited by numerical noise only in time step is reasonably small and constant. Under automatic step adjustment the return accuracy is at the limit of graphical noticibility.

void test6(): The system is similar to that in test5. However, the starting points are selected randomly. Here our interest is in the influence of friction and the regular distributed final positions (on a circle close to the (unsharpely defined) rim of the disk-shaped biotope.

void test7() Motion of the 5 heaviest bodies of the solar system in 3D.

Details on test1()

Together with the evolution of the wave function the potential is marked as a blue line and the total energy is marked as a horizontal red line. At the end of the video there is a representation of the pseudo conserved quantities norm and energy. These are named dn, de1, de2, e1, e2, where 'd' says that what is shown is the difference against the initial value where this difference is divided by the square of the time step to get values which are approximately independent of the time step. Further,'1' means real part and '2' means imaginary part. For convenient inspection of the diagrams here are links to them. Those belonging to the first video:

dn de1 de2 e1 e2

and those to the second video:

dn de1 de2 e1 e2

Notice that in the folder containing a mp4-file there are several other files which contain the main data and the main code which were used in the production of the video. In the case of the folders ./t120240528c and ./t120240528d there is the particularity that we used '#define CPM_QUAD' which means that floating point arithmetic is done with real numbers of 128 bits. This means that typical numerical noise is of the order of magnitude 1.e-32 (rather than 1.e-16 as it is for real numbers of 64 bits (type double). This confirms that de2 and e2 are 'zero up to numerical noise'. This was not the case with the old definition of energy in ./t120230607a. Surprisingly the 'QUAD-mode' works only around 20% slower than the double-mode. My expectation is that the QUAD-mode will soon become the standard mode in scientific computing.

Details on test2()

Here the most stable auto stepped method ADALF is used to generate a solution of the well known stiff ordinary differential equation x'(t) = x(t)*x(t)*(1-x(t)), x(0)=0.002. The time step can be seen to vary 'intelligently' enough to hit the correct final value 1 with great precision. With a constant timestep one can't reach the final flat phase within any reasonable execution time.

Details on test3()

The energy parts diagram shows that the pseudo-constant total energy (green horizontal line) is the sum of three contributions which are not at all nearly constant: The magenta curve shows the total kinetic energy of the particles, the black curve shows the sum of the particle pair interaction energies, and the dark blue curve, which looks like an mirror image of the kinetic energy curve is the sum of potential energies of the particles in the field of the central source. (Ignore the hardly visible yellow curve; unfortunately conversion from ppm to jpeg changes colors very badly.) The deviaton of total energy from constancy, which is not visible in the green curve due to limited graphical resolution becomes evident from the diagram relative energy error. In a program run with constant time step this diagram would look more regular. In the present run we have the following evolution of the time step. Of course, the juddery particle motion in the video reflects the variations of the auto adjusted time step.

A note concerning energy conservation in ALF-based dynamics: It came as a surprise to me to see that the expression for the observable of total energy has a quite different status for quantum particle systems and for classical particle systems. The surprise dissolves if we observe an quite obvious difference between quantum and classical mechanics: for a quantum system with known Hamiltonian knowing a wave function at a certain time allows to determine the time derivative of the state (in the Schroedinger picture) and thus a trajectory in state space. We may say: any state is a dynamical state: it determines its time evolution. In classical mechanics there are two different state-related concepts: the configurations which determines the appearance of the system at a point in time and the dynamical state (which deals not only with positions but also with velocities and determines the time evolution of the system. A different way say this is: the equation of motion is second order in differentiation with respect to time in classical mechanics and of first order for quantum mechanic. In quantum systems the definition of the dot function (see cpm1/cpmalf.h) is identical to the definition of the Hamiltonian and thus also defines an expression for the expectation value of total energy in an quantum state as represented in the ALF method. This is not the case for classical particle systems. Here we have first to convert the classical equations of motion (second order with respect to time) into a system of first order with respect to time at the price of getting twice the number of variables. My notation in code is ((x,v),(w,a)) for ((position, velocity), (internal velocity, acceleration)). The initial value for a dynamical state is given by (x,v). This determines the initial value of the ALF state as ((x,v),dot(x,v)). We write dot(x,v) as (w,a)=(dot1(x,v),dot2(x,v))=(v, acc(x,v)). Assume we have an ALF state ((x,v),(w,a)) at time t and which to compute the state ((x',v'),(w',a')) at time t'=t+dt. The ALF recipe is with h:=dt t += h; x += w*h; v += a*h; w+=(v-w)*2; a+=(acc(t,x,v)-a)*2; x += w*h; v += a*h; t += h;

In the present simulation, the particle we have been speaking about so far, will appear in two roles:

1. Its wave function, after given a uniform speed (to the right) by a kinematical boost transformation, will be taken as the initial state of our system (one free particle, no potential well in its biotope). This wave function is shown here. The edgy nature of the curves representing real and imaginary part of the wave function show that a finer discretization of position would be convenient.

2. Without boost transformation and with the potential well we have the system which will be appended to the right if any of the particles touches the right boundary of the biotope.

Both videos show how the system grows up to three particles. The log files cpmcerr...txt in the directories ./r230321a and ./r230323a show that for three particles the execution of the time evolution is achingly slow.

It should be noted that for more than 1 particle the graphics don't show the wave function but a useful derived quantity: the probability to find a particle in a place if the positions of all other particles were determined simultaneously. This gives a curve for each particle with different colors corresponding to different particles.

These computing tools were used: CPU 4.464 GHz; RAM 16 GB; OS Ubuntu; SW C++20, OpenGL, FreeGlut, Eigen3, Code::Blocks, ffmpeg (for video creation from ppm-files), ImageMagick (for conversion from ppm to jpg), C+- (self-made C++ library); Compiler GNU.