Speeding up matlab for loop

Question

I have a system of 5 ODEs with nonlinear terms involved. I am trying to vary 3 parameters over some ranges to see what parameters would produce the necessary behaviour that I am looking for.
The issue is I have written the code with 3 for loops and it takes a very long time to get the output.
I am also storing the parameter values within the loops when it meets a parameter set that satisfies an ODE event.

This is how I have implemented it in matlab.

function [m,cVal,x,y]=parameters()

b=5000;
q=0;
r=10^4; 
s=0;
n=10^-8;

time=3000;

m=[];
cVal=[];
x=[];
y=[];

val1=0.1:0.01:5;
val2=0.1:0.2:8;
val3=10^-13:10^-14:10^-11;


for i=1:length(val1)
    for j=1:length(val2)
        for k=1:length(val3)
        options = odeset('AbsTol',1e-15,'RelTol',1e-13,'Events',@eventfunction);
        [t,y,te,ye]=ode45(@(t,y)systemFunc(t,y,[val1(i),val2(j),val3(k)]),0:time,[b,q,s,r,n],options);

         if length(te)==1
             m=[m;val1(i)];
             cVal=[cVal;val2(j)];
             x=[x;val3(k)];
             y=[y;ye(1)];

         end   
        end

    end
end

Is there any other way that I can use to speed up this process?

Profile viewer results

I have written the system of ODEs simply with the a format like

function s=systemFunc(t,y,p)
        s= zeros(2,1);
        s(1)=f*y(1)*(1-(y(1)/k))-p(1)*y(2)*y(1)/(p(2)*y(2)+y(1));
        s(2)=p(3)*y(1)-d*y(2);
end

f,d,k are constant parameters.
The equations are more complicated than what's here as its a system of 5 ODEs with lots of non linear terms interacting with each other.

Answer 1

In continuation of the other suggestions, I have 2 more suggestions for you:

• You might want to try with a different solver, ODE45 is for non-stiff problems, but from the looks of it, it might seem like your problem could be stiff (parameters have a different order of magnitude). Try for instance with the ode23s method.

• Secondly, without knowing which event you are looking for, maybe it is possible for you to use a logarithmic search rather than a linear one. e.g. the Bisection method. This will severely cut down on the number of times you have to solve the equation.

Answer 2

Tommaso is right. Preallocating will save some time.

But I would guess that there is fundamentally not a lot you can do since you are running ode45 in a loop. ode45 itself may be the bottleneck.

I would suggest you profile your code to see where the bottleneck is:

profile on 
parameters(... )
profile viewer 

I would guess that ode45 is the problem. Probably you will find that you should actually focus your time on optimizing the systemFunc code for performance. But you won't know that until you run the profiler.

EDIT

Based on the profiler output and additional code, I see some things that will help

1. It seems like the vectorization of your values is hurting you. Instead of

   @(t,y)systemFunc(t,y,[val1(i),val2(j),val3(k)])

try

@(t,y)systemFunc(t,y,val1(i),val2(j),val3(k))

where your system function is defined as

function s=systemFunc(t,y,p1,p2,p3)
        s= zeros(2,1);
        s(1)=f*y(1)*(1-(y(1)/k))-p1*y(2)*y(1)/(p2*y(2)+y(1));
        s(2)=p3*y(1)-d*y(2);
end

2. Next, note that you don't have to preallocate space in the systemFunc, just combine them in the output:

    function s=systemFunc(t,y,p1,p2,p3)
       s = [ f*y(1)*(1-(y(1)/k))-p1*y(2)*y(1)/(p2*y(2)+y(1)), 
             p3*y(1)-d*y(2) ];
    end
   

Finally, note that ode45 is internally taking about 1/3 of your runtime. There is not much you will be able to do about that. If you can live with it, I would suggest increasing your 'AbsTol' and 'RelTol' to more reasonable numbers. Those values are really small, and are making ode45 run for a really long time. If you can live with it, try increasing them to something like 1e-6 or 1e-8 and see how much the performance increases. Alternatively, depending on how smooth your function is, you might be able to do better with a different integrator (like ode23). But your mileage will vary based on how smooth your problem is.

Answer 3

I have two suggestions for you.

1. Preallocate the vectors in which you store your results and use an increasing index to populate them into each iteration.
2. Since the options you use are always the same, instantiate then outside the loop only once.

Final code:

function [m,cVal,x,y] = parameters()

    b = 5000;
    q = 0;
    r = 10^4; 
    s = 0;
    n = 10^-8;
    time = 3000;
    options = odeset('AbsTol',1e-15,'RelTol',1e-13,'Events',@eventfunction);

    val1 = 0.1:0.01:5;
    val1_len = numel(val1);
    val2 = 0.1:0.2:8;
    val2_len = numel(val2);
    val3 = 10^-13:10^-14:10^-11;
    val3_len = numel(val3);
    total_len = val1_len * val2_len * val3_len;

    m = NaN(total_len,1);
    cVal = NaN(total_len,1);
    x = NaN(total_len,1);
    y = NaN(total_len,1);
    res_offset = 1;

    for i = 1:val1_len
        for j = 1:val2_len
            for k = 1:val3_len
                [t,y,te,ye] = ode45(@(t,y)systemFunc(t,y,[val1(i),val2(j),val3(k)]),0:time,[b,q,s,r,n],options);

                if (length(te) == 1)
                    m(res_offset) = val1(i);
                    cVal(res_offset) = val2(j);
                    x(res_offset) = val3(k);
                    y(res_offset) = ye(1);
                end

                res_offset = res_offset + 1;
            end
        end
    end

end

If you only want to preserve result values that have been correctly computed, you can remove the rows containing NaNs at the bottom of your function. Indexing on one of the vectors will be enough to clear everything:

rows_ok = ~isnan(y);
m = m(rows_ok);
cVal = cVal(rows_ok);
x = x(rows_ok);
y = y(rows_ok);