ODE Solver Running Very "Slowly"

Hello,

I am using an ODE solver inside a for loop to run trajectories of particles in a magnetic field, however whenever I run the code it takes the ODE solver a very long time to solve these trajecotries. On average about 3-7 seconds in each iteration. This for loop is repeated up to a 100 times and I also want to do a double for loop so there are multiple particles, so as you can see, that time adds up quickly. My issue previously is that when inputed initial conditions into the ODE solver with the following format:

icv = [x; y; z; vx; vy; vz]; %Initial conditions defined earlier within for loop
tspan = [0 4]; %Time span
[T,S] = ode15s(@bdipuniodefun, tspan, icv); %Trajectory solver

It would be hardstuck on a single iteration. I believe this is because the ODE solver is attempting to integrate over many, discrete points within the time span which slows it down tremendously. I tried fixing this issue by defining a tiertiary time step inbetween the initial and final so that it forces larger discretization:

icv = [x; y; z; vx; vy; vz]; %Initial conditions defined earlier within for loop
tspan = [0 2 4]; %Time span
[T,S] = ode15s(@bdipuniodefun, tspan, icv); %Trajectory solver

Although this does help a lot (Computing time decreased by 100x) this is still extremely slow. Any ideas why this ODE solver is slow? Or is my formatting for the initial conditions flawed so that it makes it run slow?

PS - I need to use this ODE solver within a for loop, because there is an additional process (Collisions between the particle and other particles) that needs to be solved at each time step, therfore this solver must be within the for loop. I just defined my time step to be equivalent to the time span in "tspan". If you want to see more code, I can also post more for better context.

[madhan ravi]

please upload your script file

[Jan]

It is very strange, that the computions run faster with tspan=[0,2,4] instead of [0,4]. The only difference is that the output function is called for the provided times only, but this should not affect the runtime by a factor of 100.

As long as we cannot see the code of bdipuniodefun it is impossible to guess, why it takes so much time.

You can use the profile command to find out, where the most time is spent. Please do this and provide this important information.

[Tom Keaton]

@madhan ravi I attached the files to the question. @Jan here is akso a screenshot of the profiler's results after I ran an iteration. As I said, I attached the bdipuniodefun file and B_test file that it references within it. I also attached the main script with the double for loop. I am not sure how to interpret these results, especially the difference between "Total Time" and "Self Time". Moreover, just to be more clear, I say "100x" but this is just a wild guess as when I ran it using the first initial conditions type, it just would freeze, while on the second type, it actually would complete the calculation. This is also not at all an issue when the ode solver is used outside a for loop. It runs very well and very quickly. Adding the for loop and putting the solver into it is what causes this massive slow down.

Answer 1

[Stephan]

Hi,

it is a very bad idea what happens with B_test.m while the solver runs:

In every call of the ode function by the solver B_test is called and does the same calculation with the same result. Since B_test has no input arguments you should calculate Bx By and Bz once at the beginning. Do the symbolic computation once (in a seperate matlab session) and use the resulting function handles to calculate Bx By Bz. Do it one time only and give the calculated results as extra parameters to the ode function. They should be workspace variables that are calculated one time and then be used. using symbolic calculations in this way again and again repeated is a huge waste of time, since symbolic calculations are slow. in your case they are slow and not needed. change this and you will get much faster results.

Also do not think that you have an influence on how many time step the solver makess during the calculation. using

tspan = [0 2 4]

does not mean that you save time. it means that you get three points as result back - not that the solver calculates only 3 results.

Best regards

Stephan

[Tom Keaton]

Sorry for the late response. Why is it the case then, that when I do change the tspan to be 3 discrete points, that the solver goes from completely freezing up, to actually calculating the trajectories? (even if it takes a long time) That was how I fixed that issue previously, so I could only assume that it forced discretization at those 3 points of cacluation only instead of automatically calculating a lot of solution points between just 2.

[Jan]

With tspan=[0,2,4] the integrator should calculate some more steps to reach 2 exactly instead of using the optimal step size. You can simply get the number of evaluated steps as output to control this. Maybe if you create billions of points, the iterative creation of the output exhausts the RAM and the computer uses the slow disk as virtual memory.

I'm astonished that you need to call mupadmex 28'000 times. Do you really have to obtain such a huge number of equations? Or do you solve the same equations in each time step instead of performing this once only? Sorry, I do not have the Symbolic Toolbox, so I cannot test this by my own.

[Stephan]

the much more important point is to avoid the useless use of symbolic calculations again and again.

What the documentation says:

If tspan has more than two elements [t0,t1,t2,...,tf], then the solver returns the solution evaluated at the given points. However, the solver does not step precisely to each point specified in tspan. Instead, the solver uses its own internal steps to compute the solution, then evaluates the solution at the requested points in tspan. The solutions produced at the specified points are of the same order of accuracy as the solutions computed at each internal step.

Specifying several intermediate points has little effect on the efficiency of computation, but for large systems it can affect memory management.

[Tom Keaton]

@Jan and @stephan I can't see why it is calling mudapex 28k+ times. I isolated B_test and ran this in the profiler:

[Bx, By, Bz] = B_test();

but the results look like this:

where mudapex is being called only ~600 times or is this still a lot? I am using this to calculate the B-field in a certain volume of space so I do need this many solutions I believe.

[Stephan]

note that the solver calls your function for every timestep the solver takes. your function then calls B in every single iteration. if a single iteration means 600 calls of mupadex...

[Tom Keaton]

@Stephan Ok that is what I thought, so because I do need these solutions at this level of detail, then the solution to this speed issue here would be to do this calculation once and prior to the for loop, then just have @bdipuniodefun call the closest solution point to plug into the trajectory solver for the given time step?

[Stephan]

Yes, the results that you calculate in B do not change. They do not depend on time or other time depending variables. They are constant. Then they should be calculated once and the results should be used for every calculation that works with this constant results.

If you would work in the numeric world, maybe you would not notice that there is a problem .But you are using symbolic calculations in every iteration , which makes it time consuming.

[Tom Keaton]

I see. Is there a specific resource or key terms I can look at and that you are aware of that I can use to get a headstart on figuring out how to do this?

[Stephan]

Calculate Bx, By, Bz in a seperate session and save this results in a .mat file. Then in your calculation load this .mat-file in order to provide them for the calculations. Delete the call of b_test in your ode function, since the results are already there.

[Jan]

I'd avoid a .mat file, but create the B's once and store them as persistent objects.

function [Bx, By, Bz] = B_test()
persistent Bx_ By_ Bz_
if isempty(Bx_)
  syms x y z
  mu_0_red = 10E-7;
  m = [0,0,1.28];
  r = [x, y, z];
  B = mu_0_red *(((dot(m,r)*r*3)/norm(r)^5) - m/norm(r)^3);
  Bx_ = matlabFunction(B(1));
  By_ = matlabFunction(B(2));
  Bz_ = matlabFunction(B(3));
end
Bx = Bx_;
By = By_;
Bz = Bz_;
end

[Tom Keaton]

@Jan and @Stephan thank you very much for your answers! Regardless of the method I will attempt to use (I will try both to see if there is any computing time difference), does Matlab automatically do an approximation of the position I am looking at to the nearst Bx, By, and Bz to determine the localized field? Or is that something I will also have to define in the ode function solver?

[Bjorn Gustavsson]

The Bx_, By_ and Bz_ will be functions of x, y and z, so in your equation of motion you'll have to call those functions with the current position (most likely something like rv(1:3)). I'd do it so that I'd get one function giving the magnetic field vector out instead of three functions for the separate components.

HTH

[Tom Keaton]

Got it! Thanks!

Answer 2

[Bjorn Gustavsson]

No, no, no, noo!

Do not use the ODE-suite functions to calculate trajectories of particles in B (and E) - fields. They are utterly useless at that job, they are general-purpose ODE-integration-functions and not suitable to integrate equations of motion. For charged particles in E-n-B-fields there are constants of motion that should be conserved, and the general ode-integrating schemes does not necessarily do that, so you get all sorts of unphysical results. Try for example this little script:

m_e = 9.1094e-31;        % electron mass (kg)
q_e = 1.6022e-19;        % electron charge (C)
B = 5e-5;                % Magnetic field strength (T)
w_e = (q_e*B/m_e);       % electron gyro-frequency
T_gyro = 1/(w_e/(2*pi)); % gyro-period
v0 = 5.931e+05;          % velocity of 1 eV electron (m/s)
r_gyro = v0/w_e;         % electron gyro-radius
drdvdt = @(t,rv,B,m,q) [rv(3:4);-q*B/m*rv(4);q*B/m*rv(3)]; % 2-D equation of motion x-y-plane, B || e_z
% Calculate trajectory for 100 gyro-periods, various odeNNX
[t15s,rv15s] = ode15s(@(t,rv) drdvdt(t,rv,B,m_e,q_e),linspace(0,T_gyro*100,10001),[0;-r_gyro;v_of_E(1);0]);
[t45,rv45] = ode45(@(t,rv) drdvdt(t,rv,B,m_e,q_e),linspace(0,T_gyro*100,10001),[0;-r_gyro;v_of_E(1);0]);
[t113,rv113] = ode113(@(t,rv) drdvdt(t,rv,B,m_e,q_e),linspace(0,T_gyro*100,10001),[0;-r_gyro;v_of_E(1);0]);
[t23,rv23] = ode23(@(t,rv) drdvdt(t,rv,B,m_e,q_e),linspace(0,T_gyro*100,10001),[0;-r_gyro;v_of_E(1);0]);
[t23s,rv23s] = ode23s(@(t,rv) drdvdt(t,rv,B,m_e,q_e),linspace(0,T_gyro*100,10001),[0;-r_gyro;v_of_E(1);0]);

% Plot results:
subplot(2,2,4)
plot(t15s,[(rv15s(:,3).^2+rv15s(:,4).^2),...
           (rv113(:,3).^2+rv113(:,4).^2),...
           (rv45(:,3).^2+rv45(:,4).^2),...
           (rv23(:,3).^2+rv23(:,4).^2),...
           (rv23s(:,3).^2+rv23s(:,4).^2)])
grid on
subplot(2,2,3)
plot(t15s,[(rv15s(:,1).^2+rv15s(:,2).^2),...
           (rv113(:,1).^2+rv113(:,2).^2),...
           (rv45(:,1).^2+rv45(:,2).^2),...
           (rv23(:,1).^2+rv23(:,2).^2),...
           (rv23s(:,1).^2+rv23s(:,2).^2)].^.5)
grid on
subplot(2,2,1)
plot((rv15s([1:100,end-100:end],1).^2+rv15s([1:100,end-100:end],2).^2).^.5)
plot(rv15s([1:100,end-100:end],1),rv15s([1:100,end-100:end],2))
subplot(2,2,2)
plot(rv15s([1:100,end-100:end],3),rv15s([1:100,end-100:end],4))

For this simplest of gyro-motions perpendicular to the magnetic field the electron gyro-radius should remain constant and the kinetic energy - neither of which any of the ODE-functions manage. Instead of using the ode15s or its siblings you should implement something like the Boris-mover: Boris-mover-at-wikipedia, or a symplectic scheme.

HTH

[Tom Keaton]

Hi Bjorn! Thank you for your answer! I actually checked to see if these solvers are calculating correct trajectories in the instances of constant, uniform magnetic fields and dipole fields and the solutions aligned with theory and known results so I think it is ok for this application. Morever, the Boris-mover method is also a great and better method but I am in a little bit of a time crunch with this so I am just doing the simplest method I can with the appropriate assumptions.

[Bjorn Gustavsson]

This surprises me, since my very simple toy-example for the electron-motion in a constant uniform magnetic field clearly shows that neither the kinetic energy nor the gyro-radius are conserved during 100 gyrations - both increase by ~0.2 % per gyro-period.

I do understand the concept of aproaching deadlines, so if you get acceptable results for that then OK. Do you need it for work or are you studying something like plasma-physics?

[Tom Keaton]

Sorry for the late response Bjorn. I am studying plasma physics at my university!

[Bjorn Gustavsson]

That's what I guessed. For that you might be fine with using the general ode-functions, you will definitely be fine if you check the cnosevation of energy (magnetic moment might be a bit more vactor-algebra to check). If you're OK at programming writing a Boris-mover function would take "half a day".

Good luck with your studies.