MATLAB: 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 200×10000 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!

Best Answer

  • 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.