Elliptic Integrals in loop run too slow

Hello everyone;

I'm having some trouble with a code I'm working on.

Basically the code tracks down the motion of a charged particle in different forms of an electromagnetic field (using a Runge Kutta IV order method).

In order to do so I need to evaluate the electromagnetic field components at some point in space.

I've always used the analytical expression for the components and the code showed good results, yet when it comes to a particular geometry for which I need some special functions (EllipticE/K) the code just runs too slow (about 300s for 200 steps of the particle, while I need to evaluate about 200x10000 steps).

Profiler results confirmed that this is due to using those symbolic functions, which I need to call more than 10 times at every step, yet I'm not sure about how to avoid using them.

One of my professor suggested that I could evaluate the field components at nodes of a 3D grid in space (I need to confine the particle, so I know the spatial region in which the particle will be found), store the values in 3D matrices and then evaluate the value at a certain point I need using interpolation, yet I'm not sure about how to do this.

Can anyone help me work this out or has anyone any other suggestion?

Thanks in advance!

Answer 1

[Bjorn Gustavsson]

First: scrap the Runge-Kutta! It is not a suitable scheme to solve equations of motion for charged particles in EM-fields. Instead spend 1-2 hours to implement the Boris-mover (Boris mover). With that scheme you might not get bogged down in the same way as with a adaptive RK-solver.

For the case where your fields are relatively static you could try to calculate the fields in a 3-D grid around where your particle moves, in some suitably dense grid:

[x,y,z] = mehgrid(xmin:dx:xmax,ymin:dy:ymax,zmin:dz:zmax);
[B,E] = your_EB_solver(x,y,z,t);

% Then you'd calculate the E and B-fields at position r:
B_t = interp3(x,y,z,B,r(1),r(2),r(3));
E_t = interp3(x,y,z,E,r(1),r(2),r(3));

The 3-D interpolation is not superfast either if started from scratch at every point. It might be better to use scatteredInterpolant, that function will give you a pre-calculated interpolation-function to use and reuse. Maybe you can also relax the numerical accuracy of the integration? You need a physics-related accuracy for this not a mathematical accuracy.

HTH

[Massimo Carpita]

Thank you very much!

In my case fields are static so the interpolation routine you suggested worked pretty well, yet it's still slow, I will probably need to build an interpolation function as you further suggested.

About this, I'm having some trouble with scatteredInterpolant because in its description it is said that it only accepts vectors (column) as inputs, yet I have a grid of space points and values for each B field components stored in 3D matrices:

x=linspace(xmin,xmax,100);
y=linspace(ymin,ymax,100);
z=linspace(zmin,zmax,200);

[xx yy zz]=meshgrid(x,y,z);

bx=AnalyticBX(xx,yy,zz);
by=AnalyticBY(xx,yy,zz);
bz=AnalyticBZ(xx,yy,zz);

How should I use scatterInterpolant on these data?

PS: thank you very much also for the suggestion about the scheme! I'm just a graduate student and I never heard of Boris Mover before.

[Bjorn Gustavsson]

For the scatteredInterpolant you can use the (:)-trick, something like this:

Bx_fcn = scatteredInterpolant([x(:),y(:),z(:)],bx(:),'natural');
% etc...

About the ODE-scheme - if you've taken courses in numerical methods and not heard about symplectic methods and/or Størmer-Verlet and Boris-movers then it is a sad state of scientific/academic fragmentation - when we solve equations of motion there is a rather important piece of knowledge that we have constants of motion, and using methods that somewhat preserves those is highly preferable. Good thing for you is that you learnt about this at a way younger age than I did!

[Massimo Carpita]

I managed to add the scatterInterpolant function into the code. Now it is all way faster! Thank you very much again.