How to use ode45 efficiently on a large input vector

Hello, I am using ode45 to solve a set of coupled ODEs. I need to solve them for >2M points. This is taking a long time (~45 minutes). I was wondering if I could solve the problem faster by dividing my long input vector into many smaller vectors (i.e. of length on the order of 1000) and then "stitching" the results of the ode45 solution together. I am aware that I would need to iteratively update initial conditions as I loop over subsections. However, in principle, would reducing the input vector size and running over many of those result in a speed up?

Alternatively, do you have any other advice on how to speed up the evaluation of ode45 over many points?

[James Tursa]

Have you run the profiler? Can you post your code? In general, ode45 problems can't be solved in parallel (if I understand what you are trying to do) because each step requires the solution at the previous step. The overall solution has to be done serially.

[Dorde Gluhovic]

Here is my function and a call to it

function dC = my_ODE(t,C,rhose,rhope,rhoie,Tp,ro,t_step, m)
    % t - time vector of lenght N
    % Tp - signal vector of length N
    % rhose,rhope,rhoie,ro,t_step, m - constants
    dC(1,1) =ro*( - (1 + rhose )*C(1) - 1i*2*m*C(2)*conj(C(3)));
    dC(2,1) =ro*( - (1 + rhope )*C(2) - 1i*2*C(1)*C(3)- 1i*sqrt(2*rhope)*Tp(floor(t/t_step)+1));
    dC(3,1) =ro*( - (1 + rhoie )*C(3) - 1i*2*(1-m)*C(2)*conj(C(1)));
end

f= @(t,C) my_ODE(t,C,rhose,rhope,rhoie,Tpplus,ro,t_step, m);
[tnew,C] = ode45(f,t,y0);

As you can see, equations are intercoupled, so I understand that they have to be solved sequentally (itteratively over time). I am not sure how the internal mechanics of ode45 function works, so I was thinking that maybe it could be possible to split the input vector (Tp) into chunks, solve for them (still serially, no parfor loop would be used) and stitch solutions together.

Dorde

[Bjorn Gustavsson]

To me it seems as if your problem is in that your ode-system have discontinuities whenever you step from one element in Tp to the next, i.e. whenever floor(t/t_step) ticks up. The ode-integrating functions dont like that, they expect the RHS to be continouos and even smooth. Your best approach would be to integrate it step-wise for time-periods with constant Tp.

HTH

Answer 1

[Bjorn Gustavsson]

If the ODE-system you're trying to solve is something like a set of independend equations of motion that only couple "pairwise" ( and ) and there are no coupling between the eqs of motion for different particles then it might very well be so that the ode-integration needs to take very short time-steps during one period to get sufficient accuracy for particle i and short time-steps during some other period for particle j. When you integrate all EQs simultaneously short time-steps will be taken for all of the particles every time short time-steps is required for any particle. For this type of problem it can be faster to integrate the equations of motion particle-by-particle. However if all the ODEs couple then this is not a viable route.

HTH

[Dorde Gluhovic]

Thanks for your feedback! I thought that this could be the case.

Dorde