3
$\begingroup$

I have a system of parameterized ODEs that I would like to solve using MATLAB and its ode45 solver, and was wondering if it is possible to perform such a task in parallel. For example, I have a simple model of a rocket falling in gravity, where I would like to try different parameter values of tburnstart and tburnend, which tell the program when to turn a rocket's engines on and off:

function xdot = equations_of_motion_1D_landing(t, x, G, g, M, Isp, T, tburnstart, tburnend)

    if t < tburnstart
        xdot = [
                x(3); %xdot1: x-velocity
                x(4); %xdot2: y-velocity
                -(G*M*x(1))/((x(1)^2+x(2)^2)^(3/2)); %xdot3: x-acceleration
                -(G*M*x(2))/((x(1)^2+x(2)^2)^(3/2)); %xdot4: y-acceleration
                0 %xdot5: engine mass flow rate
                ];
    else if  t >= tburnstart && t <= tburnend
        xdot = [
                x(3); %xdot1: x-velocity
                x(4); %xdot2: y-velocity
                -(G*M*x(1))/((x(1)^2+x(2)^2)^(3/2)); %xdot3: x-acceleration
                -(G*M*x(2))/((x(1)^2+x(2)^2)^(3/2))+(T/x(5)); %xdot4: y-acceleration
                -T/(g*Isp) %xdot5: engine mass flow rate
                ]; 
    else
        xdot = [
                x(3); %xdot1: x-velocity
                x(4); %xdot2: y-velocity
                -(G*M*x(1))/((x(1)^2+x(2)^2)^(3/2)); %xdot3: x-acceleration
                -(G*M*x(2))/((x(1)^2+x(2)^2)^(3/2)); %xdot4: y-acceleration
                0 %xdot5: engine mass flow rate
               ];
    end
end

Currently I test different values for tburnstart and tburnend using the following code in my main file:

%Setting up burn start and end times
tmin = 0;
tmax = 200;
tburnstart = (tmax-tmin).*rand(pop_size,1)+tmin; %burn start times
tburnend = (tmax-tburnstart).*rand(pop_size,1) + tburnstart + 10; %burn end times (must be larger than start times) (+10 makes sure tspan has more than 1 element)
pop_init = [tburnstart, tburnend]; %creating intial population


for i = 1:length(pop_init)
    tburnstart = pop_init(i,1);
    tburnend = pop_init(i,2);
    %tspan = [0:tburnend];
    tspan = 0:0.1:tburnend;
    [t,x] = ode45(@(t,x) equations_of_motion_1D_landing(t, x, G, g, M, Isp, T, tburnstart, tburnend), tspan, xinit);
end

However, using a for loop to evaluate ode45 is extremely slow, and I was wondering if there is a way to evaluate all parameters simultaneously to save on computing time.

$\endgroup$
  • 1
    $\begingroup$ Could you be more specific on "extremely slow"? How long does a single call to ode45 take? If this is already a lot, please check the stiffness of your system. Did you try with ode23s? Next, what is your population size? $\endgroup$ – GertVdE Apr 23 '17 at 10:20
  • 1
    $\begingroup$ Also try to get the if-then-else clause out of your function. You can do this since the "decision" to switch on the engine is made outside of the function. Define two functions (one with engine on, one with engine off). Integrate "engine_off" from 0 to tburnstart. Use the final value as intiial value for the second integration from tburnstart to tburnend. $\endgroup$ – GertVdE Apr 23 '17 at 10:34
  • $\begingroup$ I am running about 10000 cases of the simulation using different parameter values on fine grids, so this can take a few minutes to complete. Thanks for the if statement tip, I'll give it a go! $\endgroup$ – InquisitiveInquirer Apr 23 '17 at 19:45
  • $\begingroup$ A for loop is the way to go. You can try using parfor, but that's not guaranteed to be faster. In addition to the above suggestions, if stiffness is not an issue, you can also try relaxing the integration tolerances or using a lower order integrator like ode23 (not to be confused with ode23s, which is for stiff systems). You might also try computing G*M/(x(1)^2+x(2)^2)^(3/2) just once. $\endgroup$ – horchler Apr 24 '17 at 18:05

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.