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!
采纳的回答
Bjorn Gustavsson
2019-6-20
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
3 个评论
Massimo Carpita
2019-6-21
Bjorn Gustavsson
2019-6-21
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
2019-6-21
更多回答（0 个）
类别
在 帮助中心 和 File Exchange 中查找有关 Programming 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
