Parallelized ODE45 solution with 3D spline

Hi all,

I have a set of data describing an electric potential distribution in space V(x,y,z) with 3 independent variables x, y, and z.

I want to use this function in solving a system of ODEs describing the 3D motion of a charged particle existing in this electric potential distribution.

I have no analytic expression for V(x,y,z) and the function is ill-fitted using polynomial fits, so I use a 3D spline to construct a piecewise polynomial to be used by the ODE45 solver.

Since the 3D spline introduces a large overhead time, I do it externally in a separate function then make the spline a global variable to be accessed by the ODE solver.

My problem is that now I want to solve the system of ODEs with different combinations of initial conditions and I want to parallelize the operation using parfor loops.

However, parfor loops are not allowing the use of global variables even though they are not changed within the loop itself.

Any suggestions on how I could parallelize the operation without forsaking the use of the global variables in reconstructing the function?

Here's the code mentioned:

global spx spy spz 
 sp=csapi({z_cord,y_cord,x_cord}, Data);
 spx=fnder(sp, [0,0,1]);
 spy=fnder(sp, [0,1,0]);
 spz=fnder(sp, [1,0,0]);
 
 
y0=[-12:0.01:-11];
for j=1:length(y0)
tstart=tic;
[t,sol]=ode45(@func_xy,[0 30e-6],[0;0;y0(j);0;0;0], odeset('Events', @(t,y) myEvent(t,y,tstart)));
end

function dxdt=func_xy(t,y)
global spx spy spz 
q=1.6e-19;
m=400*1.6e-27;
dxdt=[y(2);
    (-q/m)*fnval(spz,{y(1), y(3), y(5)}) ;
     y(4);
    (-q/m)*fnval(spy,{y(1), y(3), y(5)});
     y(6);
    (-q/m)*fnval(spx,{y(1),y(3), y(5)});];


end
function [value, isterminal, direction] = myEvent(t, y,tstart)
value      = double(toc(tstart)>200);
isterminal = 1;   % Stop the integration
direction  = 0;
end

Answer 1

[J. Alex Lee]

Can you just pass the variables spx spy spz to func_xy as additional arguments?

q=1.6e-19;
m=400*1.6e-27;

odeOpts = odeset('Events', @(t,y) myEvent(t,y,tstart))

% pass additional arguments by using function handle:
odefunHandle = @(t,y)func_xy(t,y , q,m,spx,spy,spz)

[t,sol]=ode45(odefunHandle,[0 30e-6],[0;0;y0(j);0;0;0], odeOpts);

%%%

function dxdt=func_xy(t,y , q,m,spx,spy,spz)
% global spx spy spz 
% q=1.6e-19;
% m=400*1.6e-27;
dxdt=[y(2);
    (-q/m)*fnval(spz,{y(1), y(3), y(5)}) ;
     y(4);
    (-q/m)*fnval(spy,{y(1), y(3), y(5)});
     y(6);
    (-q/m)*fnval(spx,{y(1),y(3), y(5)});];
end

By the way are you sure that the main cost of the spline interp is the up-front cost of csapi/fnder, and not the fnval()s?

[Hana Medhat]

Thank you for such a simple solution ! It just needed a fresh set of eyes.

The csapi function is extremely exhaustive in terms of memory & time especially for the dense matrix of data I am using.

I am sure fnval is not optimal but calling csapi with each time step is certainly not feasible.

[J. Alex Lee]

Glad that it works!

Out of curiosity, is the splines you're using any different (or better performance if same) as the base Matlab's interpolant objects (griddedInterpolant, scatteredInterpolant) when using the spline method? (disclaimer, I am not that knowledgeable about splines, except that there are lots of different types of splines).

[Hana Medhat]

I played with interp3 a little bit to see if it's a better option. It seems that it handles dense data better in terms of memory and also provides the cubic spline method among other options.

However it doesn't allow you to have the spline coefficients as a structure so in my case I have to call interp3() in every time step the ode45 solver takes which is extremely slow.

csapi() is exhaustive in terms of memory because it returns the spline coefficients as a structure, but the fnval() is then much faster than calling interp3().

[J. Alex Lee]

On closer reading, I just realized your spx,spy,spz are gradient components, so it is necessary anyway for you to have access to the spline points themselves to carry out the derivatives.

But anyway, just for reference, maybe check out

https://www.mathworks.com/help/matlab/ref/scatteredinterpolant.html#bvkm1tv-6

the "interpolant" class (scatteredInterpolant or griddedInterpolant) is supposed to solve exactly the issue of needing to evaluate interpolations multiple times. On a brief search, I couldn't find a way to efficiently calculate derivative coefficients based on the interpolant though. Seems like a missed opportunity unless it's intentional to make us buy toolboxes.