Asked by Tom Keaton
on 10 Jan 2019

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.

Answer by Stephan
on 11 Jan 2019

Edited by Stephan
on 11 Jan 2019

Accepted Answer

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
on 17 Jan 2019

Bjorn Gustavsson
on 18 Jan 2019

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
on 24 Jan 2019 at 19:27

Got it! Thanks!

Sign in to comment.

Answer by Bjorn Gustavsson
on 17 Jan 2019

Edited by Bjorn Gustavsson
on 17 Jan 2019

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

Bjorn Gustavsson
on 18 Jan 2019

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
on 24 Jan 2019 at 19:27

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

Bjorn Gustavsson
on 25 Jan 2019 at 8:33

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.

Sign in to comment.

Opportunities for recent engineering grads.

Apply Today
## 3 Comments

## madhan ravi (view profile)

Direct link to this comment:https://nl.mathworks.com/matlabcentral/answers/439254-ode-solver-running-very-slowly#comment_659261

## Jan (view profile)

Direct link to this comment:https://nl.mathworks.com/matlabcentral/answers/439254-ode-solver-running-very-slowly#comment_659351

## Tom Keaton (view profile)

Direct link to this comment:https://nl.mathworks.com/matlabcentral/answers/439254-ode-solver-running-very-slowly#comment_659508

Sign in to comment.